
import numpy as np
from math import gamma
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import pandas as pd

import os
os.chdir("C:\\Users\\user\\Downloads\\Thesis_paper_1\\calibration")

np.set_printoptions(suppress=True) # supress scientific notation

import matplotlib.ticker as mtick
import locale
locale.setlocale(locale.LC_NUMERIC, "de_DE")
plt.rcParams['axes.formatter.use_locale'] = True


## parameters



def get_data(year):
    w='data_' + str(year)
    alfa = np.array(pd.read_excel(f'{w}.xlsx').iloc[:,5])
    phi = np.array(pd.read_excel(f'{w}.xlsx').iloc[:,4])
    return alfa, phi



def pars(year):
    global alfa, rho, beta, theta, eta, psi, kappa, lamb, gamma1, G, phi
    rho = 0.19
    beta = 0.69
    theta = 3.44
    eta = 0.25
    psi = eta**eta*(1-eta)**(1-eta)
    kappa = np.divide(1, (1 - eta))

    lamb = 1 - np.divide(kappa, theta)
    gamma1 = np.multiply(eta, kappa)
    
    G = gamma(1 - np.divide(kappa, np.multiply(theta, (1-rho)) ) )    
    
    alfa, phi = get_data(year)


# taus

st=3
def taus2(st):        
    tau = np.random.uniform(low=-0.989, high=5e4, size=(st)) 
#    tau[0] = 0
    x1 = np.array(tau)    
    return x1


# Time spent at school

def sf( ):
    global s
    s = np.power( (1 + np.divide( (1-eta), ( np.multiply(beta, phi) ) ) ), -1 )
    #s = s.reshape(st, 1)
    return s


## N_i = H_i / w_i^gamma1


def N_if(x1):
    s = sf()
    T = np.multiply(np.power((1 + x1), (-1*theta*lamb-gamma1)), np.power(psi, theta*lamb) )    
    S = np.power(s, (phi*(theta*lamb + 1 + gamma1 ) ) )
    S2 = np.multiply(np.power( (1 - s), (theta*lamb/kappa*beta)  ), np.power(eta, gamma1) )
    N_i = np.multiply(G, np.multiply(T, np.multiply(S, S2)))
    return N_i


#  linear system 
pars(92)

A = np.repeat(alfa, st).reshape(st, st).T  # matriz de parâmetros    
B = np.repeat(alfa, st).reshape(st, st).T        # Matriz de alfas 
for i in range(st):
    for j in range(st):
        if i==j:
            A[i, j] = 1 - (A[i, j] - 1)*(theta*lamb + gamma1)
            B[i, j] = B[i, j] - 1
        if i != j:
            A[i, j] = -A[i, j]*(theta*lamb + gamma1) 

    
def lin_sys(x1):
    global N_i
    N_i = np.log( N_if(x1) )          # Logaritmo de N_i
    C = np.log(alfa)                  # Logaritmo dos alfas
    D = np.dot(B, N_i) + C            # produto matricial      
    lin = np.linalg.solve(A, D)       # Resolvendo o sistema  A w = D
    res = np.exp(lin)
    return res


# w_til


def w_til_f(x1):
    s = sf()
    A = np.multiply(np.power((1 + x1), -eta), w )
    S = np.multiply(np.power(s, phi), np.power((1 - s), ((1-eta)/beta) ) )    
    w_til = np.multiply(psi, np.multiply(A, S) )    
    return w_til


# p_i

def p_i_f(x1):
    w_til = w_til_f(x1)
    p_i =  np.power(w_til, theta)/ np.sum( np.power(w_til, theta) )   
    return p_i


# wages - model 

def W_i_f(x1):
    w_til = w_til_f(x1)
    sum_wtil = np.power( np.sum( np.power(w_til, theta) ) , (kappa/theta) ) 
    S = np.power((1-s), (-1/beta))
    W_i = np.multiply(kappa, np.multiply(G, np.multiply(S, sum_wtil) ) )    
    return W_i


# PNAD - DATA


def get_p_i(year):
    global p_s2
    w='data_' + str(year)
    p_i = np.array(pd.read_excel(f'{w}.xlsx').iloc[:,2])
    p_s2 = p_i #[1:st]
    return p_i


def get_W_i(year):
    global W_id
    w='data_' + str(year)
    W_id = np.array(pd.read_excel(f'{w}.xlsx').iloc[:,1])
    return W_id


get_W_i(15)

array([4.16695872, 5.69981764, 5.62506513])


# Objective function



def obj2(x1):    
    global w
    x1 = x1.reshape((st)) 
#    x1[0] = 0 
    w = lin_sys(x1) 
    W_i = W_i_f(x1) 
    p_i = p_i_f(x1)   
    #p_i = p_i[1:st]    
    #D = np.sum(np.power( np.divide(( W_id - W_i ), W_id )  , 2 ) ) 
    D = np.sum(np.power( np.divide(( p_s2 - p_i ), p_s2 )  , 2 ) )
#    D = np.log(D)
    return D


get_p_i(15)
obj2(taus2(3))

1.2617985663077622


# bounds

Bd = ((-0.99, 5e4), )*st 
Bd = np.array(Bd, dtype=object)


# callback


cc = 0
def callback(x):
    global cc
    cc += 1
    fobj = obj2(x)
    print(f'\033[1;033mObjetivo: {np.around(fobj, 5)}, iter: {cc}')


pars(15)
get_W_i(15)
get_p_i(15)
st = 3
cc=0
sol= minimize(obj2, taus2(3),  method='L-BFGS-B', callback=callback,
              bounds=Bd,
              tol=1e-9, options={'maxiter':2e4})


# solver 1
pars(92)
get_W_i(92)
get_p_i(92)
st = 3
cc=0
sol= minimize(obj2, sol.x,  method='Nelder-Mead', callback=callback,
              tol=1e-9, options={'maxiter':2e4})


#z1 = pd.DataFrame({'dist_92': d_92, 'dist_95': d_95, 'dist_2':d_2, 'dist_15':d_15})

#ocp = np.array(['AGR', 'IND', 'SERV'])
#z1.index = ocp

#z1.to_excel('dist.xlsx')


# occupations

#ocp = np.array(['AGR', 'IND_EXTR', 'IND_TRANSF', 'IND_CONST',
#                'SERV_COMER', 'SERV_TRANSP_ARM',  'SERV_FIN_IMOB', 'ADM_SAUDE_EDUC'])

#names = np.array(pd.read_excel('IPEA_PIBS.xls', sheet_name='nomes', header=None)).reshape(1, 8)



ocp = np.array(['AGR', 'IND', 'SERV'])
names = np.array(['Agropecuária', 'Indústria', 'Serviços'])

for i in range(st):
    print(f'{names[i]:-<20}> {ocp[i]}')

Agropecuária--------> AGR
Indústria-----------> IND
Serviços------------> SERV


# get distortions


def get_dist(year):
    global w
    pars(year)
    get_p_i(year)    
    l ='dist_' + str(year)
    dist = np.array(pd.read_excel('dist.xlsx').loc[:,l])
    w = lin_sys(dist)
    return dist


# Distortions of 1992, 1995, 2002 and 2015

pd.read_excel('dist.xlsx').iloc[:,0:5]


def g_pif(x):
    p_i = p_i_f(get_dist(x))
    p_s = get_p_i(x)       
    fig,ax = plt.subplots(figsize=(10,8))
    plt.scatter(p_i, p_s)
    for tt, txt in enumerate(names):
        plt.annotate(txt, (p_i[tt], p_s[tt]), size=18)  

    plt.plot([0.0, max(p_s)+0.1], [0.0, max(p_s)+0.1], 'k-', lw=2, label='linha de 45°', color='blue')
    plt.grid(True)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.xlabel(r"$p_i$ - Modelo", fontsize=20)
    plt.ylabel(r"$p_i$ - PNAD", fontsize=20)
    plt.legend(loc="lower right", prop={'size': 20})
    plt.tight_layout()  
    #plt.savefig("C:/Users/user/Downloads/Thesis_paper_1/paper_tex/fig5.eps", format='eps', dpi=2000)


g_pif(15)


def Wi_plot(year):
    N_i = N_if(get_dist(year))
    H_i = np.multiply(N_i, np.power(w, gamma1) ) 
    fig,ax = plt.subplots(figsize=(10,8))
    plt.scatter(H_i, W_i_f(get_dist(year)) )
    for tt, txt in enumerate(ocp):
        plt.annotate(txt, (H_i[tt], W_i_f(get_dist(year))[tt] ), size=15)  
    plt.grid(True)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.xlabel(r"$H_i$ - Modelo", fontsize=20)
    plt.ylabel(r"$W_i$ - Modelo", fontsize=20)
    plt.tight_layout()  
    
    
Wi_plot(15)


# human capital

def hc_plt(year):
    N_i = N_if(get_dist(year))
    H_i = np.multiply(N_i, np.power(w, gamma1) ) 
    fig,ax = plt.subplots(figsize=(10,8))
    plt.scatter(H_i, w, s=50)
    for tt, txt in enumerate(ocp):
        plt.annotate(txt, (H_i[tt], w[tt]), size=15)  
    plt.grid(True)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.xlabel(r"$H_i$ - Model", fontsize=20)
    plt.ylabel(r"$w_i$ - Model", fontsize=20)
    plt.tight_layout() 

hc_plt(15)


def hc_plt2(year):
    N_i = N_if(get_dist(year))
    H_i = np.multiply(N_i, np.power(w, gamma1) )
    
    N_i2 = N_if(np.repeat(get_dist(year)[2], st))
    H_i2 = np.multiply(N_i2, np.power(w, gamma1) ) 
    fig,ax = plt.subplots(figsize=(10,8))
    plt.scatter(H_i, H_i2, s=10)
    for tt, txt in enumerate(names):
        plt.annotate(txt, (H_i[tt], H_i2[tt]), size=18)  
    plt.plot([0.0, max(H_i)], [0.0, max(H_i)], 'k-', lw=2, label='linha de 45°', color='blue')
    plt.grid(True)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.xlabel(r"$H_i$ - Antes da mudança", fontsize=20)
    plt.ylabel(r"$H_i$ - Depois da mudança", fontsize=20)
    plt.legend(loc="lower right", prop={'size': 20})
    plt.tight_layout() 
    #plt.savefig("C:/Users/user/Downloads/Thesis_paper_1/paper_tex/fig7.eps", format='eps', dpi=2000)

hc_plt2(15)


# p_i

def pi_plot(year):
    fig,ax = plt.subplots(figsize=(10,8))
    plt.scatter( p_i_f(get_dist(year)), W_i_f(get_dist(year)), s=20)
    for tt, txt in enumerate(ocp):
        plt.annotate(txt, (p_i_f(get_dist(year))[tt], W_i_f(get_dist(year))[tt]), size=15)  
    plt.grid(True)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.xlabel(r"$p_i$ - Model", fontsize=20)
    plt.ylabel(r"$W_i$ - Model", fontsize=20)
    plt.tight_layout()


pi_plot(15)


# output

t_y = np.arange(0, 1000, 0.1)
Y_tau = []

for i  in t_y:
    res = get_dist(15) + i    
    N_i = N_if(res)
    H_i = np.multiply(N_i, np.power(w, gamma1) ) 
        
    Y_tau.append( np.prod(np.power(H_i, alfa )) )
    
def pct(a):
   return ((a - a[0])/a[0] )*100



fig, ax = plt.subplots(figsize=(10,8))
ax.plot(t_y, pct(Y_tau))
plt.grid(True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel(r"$\tau_i$", fontsize=20)
plt.ylabel(r"Variação percentual em $Y$ ", fontsize=20)
plt.tight_layout() 
#plt.savefig("C:/Users/user/Downloads/Thesis_paper_1/paper_tex/fig6.eps", format='eps', dpi=2000)


# alphas
get_dist(15)

fig,ax = plt.subplots(figsize=(10,8))
plt.scatter( alfa, W_i_f(get_dist(15)), s=20)
for tt, txt in enumerate(ocp):
    plt.annotate(txt, (alfa[tt], W_i_f(get_dist(15))[tt] ), size=15)  
plt.grid(True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel(r"$\alpha_i$ - Model", fontsize=20)
plt.ylabel(r"$W_i$ - Model", fontsize=20)
#plt.tight_layout()

Text(0, 0.5, '$W_i$ - Model')


# spend in education

def geduc_f(x1):
    P = np.power(np.divide(1, p_i_f(x1)) , kappa/eta )
    N = np.power(np.divide(np.multiply(w, np.power(sf(), phi)), (1 + x1) ), kappa )
    G = np.multiply(np.power(eta, kappa), gamma1 )
    educ = np.multiply(np.multiply(P, N), G)
    return educ


geduc_f(taus2(3))

array([0.00044847, 0.00000182, 0.04363348])


t_y = np.arange(0, 100, 0.1)
h_tau = []
spend = []

for i  in t_y:
    res = get_dist(15) - i    
    w = lin_sys(res)
    N_i = N_if(res)
    H_i = np.sum(np.multiply(N_i, np.power(w, gamma1) ) )
    educ = np.sum(geduc_f(res) )
    
    
    spend.append(educ)
    h_tau.append( H_i )
    
def pct(a):
   return ((a - a[0])/a[0] )*100



fig, ax = plt.subplots(figsize=(10,8))
ax.plot(t_y, pct(h_tau), label='H')
ax.plot(t_y, pct(spend), label='e')
plt.grid(True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel(r"$\tau_i$", fontsize=20)
plt.ylabel(r"Variação percentual em $H$ ", fontsize=20)
plt.tight_layout() 
plt.legend(loc="lower right", prop={'size': 20})

#plt.savefig("C:/Users/user/Downloads/Thesis_paper_1/paper_tex/fig6.eps", format='eps', dpi=2000)

<matplotlib.legend.Legend at 0x242563191f0>


%notebook "C:\Users\user\Downloads\Thesis_paper_1\calibration\file1.ipynb"

	Unnamed: 0	dist_92	dist_95	dist_2	dist_15
0	AGR	64353.362239	63953.101872	52486.467244	46081.801140
1	IND	13422.416752	13059.938405	16113.739739	17391.597780
2	SERV	9479.312406	10685.841401	10083.784299	10650.878049

Results¶

Fit Quality¶

Human capital¶

Proportion of people in each occupation¶

Increases in distortions and effects on GDP¶

Alphas and w¶