Considere que você deseja estimar os parâmetros da equação (1) abaixo que pertence ao sistema de equações simultâneas:
\(y_1=\alpha_1 + \beta_1 y_2 + \gamma_1 z_1 + u_1\; (1)\)
\(y_2= \alpha_2 + \beta_2 y_1 + \gamma_2 z_2 + u_2\; (2),\)
onde \(y_1\) e \(y_2\) são endógenas e \(z_1\) e \(z_2\) são exógenas.
Podemos reescrever esse modelo da seguinte forma (isolando os valores de \(y_1\) e \(y_2\)):
\(y_1 (1-\beta_1 \beta_2)=\alpha_1 + \beta_1 \alpha_2 + \beta_1 \gamma_2 z_2 + \gamma_1 z_1 + \beta_1 u_2 + u_1\; (3) \)
\(y_2 (1-\beta_1 \beta_2)=\alpha_2 + \beta_2 \alpha_1 + \beta_2 \gamma_1 z_1 + \gamma_2 z_2 + \beta_2 u_1 + u_2 \; (4)\)
Vamos proceder fazendo simulações Monte Carlo. Com o conhecimento dos valores dos \(\alpha\)s, dos \(\beta\)s e dos \(\gamma\)s e de amostras das variáveis \(z\)s e \(u\)s geramos várias amostras de variáveis \(y\)s usando as equações (3) e (4) acima.
No código abaixo, implementamos dois estimadores:
a) Métodos de mínimos quadrados (que é viesado) supondo que estamos estimando a Equação (1) independente do conhecimento da segunda.
b) Estimador de variáveis instrumentais a partir do método dos momentos e supondo que
\(E[u_1]=0\; (5)\)
\(E[z_1 u_1]=0\; (6)\)
\(E[z_2 u_1]=0\; (7).\)
Note que para a Equação (7) ser válida (ou seja, dar informação adicional as equações (5) e (6)), nós precisamos que o coeficiente \(\delta_2\) na Equação (8) abaixo seja significante
\[y_2=\delta_0 + \delta_1 z_1 + \delta_2 z_2 + \epsilon.\]
Veja abaixo também que esse coeficiente ser significante implica que \(z_2\) consegue recuperar a casualidade de \(y_2\) para \(y_1\).
A figura abaixo apresenta as estimativas do parâmetro \(\beta_1\) na Equação (1), onde a linha vermelha apresenta o valor correto desse parâmetro. Note que o estimador (a) [figura da esquerda] é viesado e o estimador (b) [figura da direita]não é:
O código usado para implementar essa figura está apresentado a seguir:
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
plt.figure(num=None, figsize=(8, 4), dpi=80, facecolor='w', edgecolor='k')
def geracaoDados(n):
alpha1=0.5
alpha2=0.7
beta1=-1
beta2=1
gamma1=-0.5
gamma2=2
z1=np.random.normal(2,1,n)
z2=np.random.normal(2,1,n)
u1=np.random.normal(0,1,n)
u2=np.random.normal(0,1,n)
y1=((alpha1+beta1*alpha2)*np.ones(n)+beta1*gamma2*z2+gamma1*z1+beta1*u2+u1)/(1-beta1*beta2)
y2=((alpha2+beta2*alpha1)*np.ones(n)+beta2*gamma1*z1+gamma2*z2+beta2*u1+u2)/(1-beta1*beta2)
return y1,y2,z1,z2
def estimacaoVarInstrumental(y1,y2,z1,z2):
sumZ1=np.sum(z1)
sumZ2=np.sum(z2)
sumY1=np.sum(y1)
sumY2=np.sum(y2)
sumY1Z1=np.dot(y1,z1)
sumY2Z1=np.dot(y2,z1)
sumY1Z2=np.dot(y1,z2)
sumY2Z2=np.dot(y2,z2)
sumZ1Z1=np.dot(z1,z1)
sumZ1Z2=np.dot(z1,z2)
A=np.array([[n,sumY2,sumZ1],[sumZ1,sumY2Z1,sumZ1Z1],[sumZ2,sumY2Z2,sumZ1Z2]])
b=np.array([sumY1,sumY1Z1,sumY1Z2])
beta=np.linalg.solve(A, b)
return beta
def estimacaoOls(x,y):
hatBeta = np.linalg.lstsq(x,y)[0]
return hatBeta
if __name__ == '__main__':
numeroAmostras=1000
n=50 # numero de observacoes
# Estimacao OLS
vetorOLS=np.empty([numeroAmostras])
for i in range(0,numeroAmostras):
[y1,y2,z1,z2]=geracaoDados(n)
x=np.vstack([np.ones(n),y2,z1]).T
hatBetaOls=estimacaoOls(x,y1)
vetorOLS[i]=hatBetaOls[1]
axes=plt.subplot(121)
plt.setp(axes, xticks=[-1.3,-1.1,-0.9,-0.7,-0.5,-0.3,-0.1,0.1],yticks=[20,40,60,80,100,120,140])
plt.hold(True)
plt.plot([-1, -1],[0,100],'r-')
plt.hist(vetorOLS,bins=30)
plt.title('OLS: Amostra de tamanho 1000')
# Estimacao Variavel Instrumental
vetorInst=np.empty([numeroAmostras])
for i in range(0,numeroAmostras):
[y1,y2,z1,z2]=geracaoDados(n)
hatBetaInst=estimacaoVarInstrumental(y1,y2,z1,z2)
vetorInst[i]=hatBetaInst[1]
axes=plt.subplot(122)
plt.hold(True)
plt.plot([-1, -1],[0,100],'r-')
plt.hist(vetorInst,bins=30)
plt.title('I: Amostra de tamanho 1000')