#!/usr/bin/env python
# coding: utf-8

# #  Tutorial: CNLS_VRS model

# - Author   : Sheng Dai (sheng.dai@aalto.fi)
# - Objective: Estimates a variable returns to scale (VRS) production function model and overall economic efficiency
# - Data     : data genrate pprocess (DGP) 
# - Date     : April 16, 2020

# In[1]:


import CNLS_VRS
import numpy as np
import pandas as pd
from cvxopt import matrix


# In[2]:


## Generate data
## Number of DMUs
m = 50     

## set seed
np.random.seed(1)

## x1~U[1,10]
x1 = np.random.uniform(1, 10, m)

## composite error term
#  v ~  N(0, 0.708^2)
#  u ~ |N(0, 1.174^2)|
y = 3 + x1**0.5+np.log(x1)+np.random.normal(0,0.708,m)- abs(np.random.normal(0,1.174,m))

x = np.asmatrix(x1).T


# In[3]:


## calculate the VRS model; 
##           call function cnls_vrs()
res_cnls = CNLS_VRS.cnls_vrs(x, y, 'MOSEK')

## results: yhat and beta
yhat  = res_cnls[:m]
beta  = res_cnls[m:2*m]


# In[4]:


## calculate the residuals (eps)
eps   = matrix(0.0, (m,1))
for i in range(m):
    eps[i]   = y[i] - yhat[i]  


# In[5]:


## calculate the constant term (alpha)
alpha = matrix(0.0, (m,1))
for i in range(m):
 alpha[i] = y[i] - eps[i] - beta[i] * x1[i]


# In[6]:


alpha = pd.DataFrame(alpha)
beta  = pd.DataFrame(beta)
eps   = pd.DataFrame(eps)

results = pd.concat([alpha, beta, eps], axis=1) 

var_name = ['Alpha', 'Beta', 'Residual']

results.columns = var_name


# In[7]:


## print estimates
print(results)