Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

OLS versus IV em pequenas amostras: Evidências Usando Simulações de Monte Carlo

0 votos
36 visitas
perguntada Out 28 em Economia por Thiago Trafane (36 pontos)  
editado Out 28 por Thiago Trafane

Trabalho final do curso de Métodos numéricos e modelos computacionais em
economia.

Compartilhe

2 Respostas

0 votos
respondida Out 28 por Thiago Trafane (36 pontos)  
editado Out 28 por Thiago Trafane

1. Introdução

Na presença de regressores endógenos, o estimador de OLS é viesado e inconsistente. Nesse contexto, a estratégia econométrica usual é usar variáveis instrumentais (IV), aplicando métodos como 2SLS ou GMM, que são consistentes.

Um grande desafio é que a superioridade de estimadores de IV frente ao OLS é apenas assintótica. Contudo, nada garante que a amostra seja grande o suficiente em problemas do mundo real. E, no caso de amostras pequenas, não se sabe qual estimador apresenta desempenho superior. Intuitivamente, imagina-se que o OLS apresentaria maior viés e menor variância frente aos estimadores de IV, de modo que a performance relativa em termos de Mean Squared Error (MSE) seria inconclusiva.

O objetivo do trabalho é investigar o comportamento desses estimadores em amostras pequenas. Para isso, eu calculo uma proxy do MSE do OLS e do estimador de IV em amostras pequenas usando simulações de Monte Carlo. Eu considerei um problema exatamente identificado, isto é, em que o número de instrumentos é igual ao número de regressores. Nesse caso, os estimadores de 2SLS e GMM são exatamente equivalentes e são normalmente chamados simplesmente de estimador de variáveis instrumentais.

Mais especificamente, o trabalho visa (i) avaliar o desempenho dos estimadores quando variamos o grau de endogeneidade e a força do instrumento e (ii) testar uma proposta de método de escolha do estimador que pode ser aplicada em problemas práticos, que consiste em escolher o estimador de menor MSE assintótico estimado.

Conforme esperado, os resultados apontam que o estimador de IV tem desempenho melhor em comparação ao OLS quanto maior for endogeneidade e mais forte for o instrumento, tudo o mais permanecendo constante. Além disso, as simulações mostram resultados promissores para o método de escolha proposto, especialmente em amostras pequenas.

2. Metodologia

2.1. Setup

A equação básica do problema é
\[y = \beta z + \varepsilon\]
em que \(\beta=5\), \(E(z)=0\), \(\varepsilon \sim N(0,1)\) e \(corr(\varepsilon,z) \neq 0\).

Sem perda de generalidade dado o Teorema da Projeção Ortogonal, podemos escrever o regressor endógeno como
\[z = \bar{x} + \theta_z \varepsilon\]
em que \(\varepsilon \perp \bar{x}\), \(\theta_z \neq 0\) e \(\bar{x} \sim N(0,1)\).

É fácil mostrar que \(\bar{x}\) é a melhor IV possível, isto é, a variável exógena mais correlacionada com \(z\). Para ver isso, tome uma variável \(\tilde{x} \perp \varepsilon\) qualquer. Então,
\[Cov(\tilde{x},z) = Cov(\tilde{x},\bar{x}) + \theta_z Cov(\tilde{x},\varepsilon) = Cov(\tilde{x},\bar{x})\]

\[corr(\tilde{x},z) = \left( \frac{\sigma_{\tilde{x}} \sigma_{\bar{x}}}{\sigma_{\tilde{x}} \sigma_z} \right) corr(\tilde{x},\bar{x}) = \left( \frac{\sigma_{\bar{x}}}{\sigma_z} \right) corr(\tilde{x},\bar{x}) \]

que é maximizado quando \(\tilde{x} = \bar{x}\). Assim, \(\bar{x}\) é a IV ideal ou ótima.

Na prática, podemos não observar \(\bar{x}\), mas sim um instrumento \(x\), que é apenas uma medida imperfeita de \(\bar{x}\):
\[x = \bar{x} + \theta_x \omega\]

em que \(\omega \perp \bar{x}\), \(\omega \perp \varepsilon\) e \(\omega \sim N(0,1)\).

2.2. Usando correlações para identificar parâmetros

Ao invés de fazer as simulações atribuindo diretamente valores para os parâmetros \(\theta_z\) e \(\theta_x\), acredito que seja mais intuitivo definir correlações, obtendo, então, tais parâmetros.

Assim, para obter \(\theta_z\), eu uso a \(corr(\varepsilon,z)\). Do que foi apresentado na seção 2.1, é fácil ver que \(\sigma_z^2 = \sigma_{\bar{x}}^2 + \theta_z^2 \sigma_{\varepsilon}^2\) e \(Cov(\varepsilon,z) = \theta_z \sigma_{\varepsilon}^2\). Logo,
\[corr(\varepsilon,z) \equiv \frac{Cov(\varepsilon,z)}{\sigma_{\varepsilon} \sigma_z} = \frac{\theta_z \sigma_{\varepsilon}^2}{\sigma_{\varepsilon} \sqrt{\sigma_{\bar{x}}^2 + \theta_z^2 \sigma_{\varepsilon}^2}} = \frac{1}{\sqrt{\frac{\sigma_{\bar{x}}^2}{\theta_z^2 \sigma_{\varepsilon}^2} + 1}} \]

\[ \theta_z^{-2} = \frac{1-corr(\varepsilon,z)^2}{ (\sigma_{\bar{x}}^2/\sigma_{\varepsilon}^2) corr(\varepsilon,z)^2}\]

\[\theta_z = \pm \sqrt{\frac{(\sigma_{\bar{x}}/\sigma_{\varepsilon}) corr(\varepsilon,z)}{1-corr(\varepsilon,z)^2}} = \frac{(\sigma_{\bar{x}}/\sigma_{\varepsilon}) corr(\varepsilon,z)}{\sqrt{1-corr(\varepsilon,z)^2}} = \frac{ corr(\varepsilon,z)}{\sqrt{1-corr(\varepsilon,z)^2}}\]

em que na última passagem eu usei \(\sigma_{\bar{x}}=\sigma_{\varepsilon}=1\). Na penúltima passagem eu usei o fato do sinal da correlação ser dado pelo sinal da covariância, o que implica que \(\theta_z\) deve ser positivo (negativo) quando a correlação for positiva (negativa).

No caso de \(\theta_x\), existem duas opções de correlação: \(corr(x,\bar{x})\) ou \(corr(x,z)\). No primeiro caso, usamos o que foi apresentado na seção 2.1 para obter \(\sigma_x^2 = \sigma_{\bar{x}}^2 + \theta_x^2 \sigma_{\omega}^2\) e \(Cov(x,\bar{x}) = \sigma_{\bar{x}}^2\). Então,
\[corr(x,\bar{x}) \equiv \frac{Cov(x,\bar{x})}{\sigma_x \sigma_{\bar{x}}} = \frac{\sigma_{\bar{x}}^2}{\sigma_{\bar{x}} \sqrt{\sigma_{\bar{x}}^2 + \theta_x^2 \sigma_{\omega}^2}} = \frac{1}{\sqrt{\frac{\theta_x^2\sigma_{\omega}^2}{ \sigma_{\bar{x}}^2} + 1}} \]

\[ \theta_x = \pm \sqrt{\frac{1-corr(x,\bar{x})^2}{ (\sigma_{\omega}^2/\sigma_{\bar{x}}^2) corr(x,\bar{x})^2}} = \frac{\sqrt{1-corr(x,\bar{x})^2}}{ (\sigma_{\omega}/\sigma_{\bar{x}}) corr(x,\bar{x})} = \frac{\sqrt{1-corr(x,\bar{x})^2}}{corr(x,\bar{x})}\]

em que na última passagem eu usei \(\sigma_{\omega}=\sigma_{\bar{x}}=1\). Note que só podemos considerar correlações positivas, já que a covariância é sempre positiva nesse caso. Assim, o sinal de \(\theta_x\) é irrelevante, o que justifica a penúltima passagem.

No segundo caso, usamos novamente o que foi apresentado na seção 2.1 para obter \(Cov(x,z) = Cov(x,\bar{x}) + \theta_z Cov(x,\varepsilon) = Cov(x,\bar{x}) = \sigma_{\bar{x}}^2\). Então, usando \(\sigma_x^2\) e \(\sigma_z^2\) obtidos anteriormente, temos que
\[corr(x,z) \equiv \frac{Cov(x,z)}{\sigma_x \sigma_z} = \frac{\sigma_{\bar{x}}^2}{\sqrt{\sigma_{\bar{x}}^2 + \theta_x^2 \sigma_{\omega}^2} \sqrt{\sigma_{\bar{x}}^2 + \theta_z^2 \sigma_{\varepsilon}^2}} = \frac{1}{\sqrt{\frac{\theta_x^2 \sigma_{\omega}^2}{\sigma_{\bar{x}}^2} + 1} \sqrt{\frac{\theta_z^2 \sigma_{\varepsilon}^2}{\sigma_{\bar{x}}^2}+1}} \]

\[\theta_x = \pm \sqrt{\frac{\frac{1}{(\sigma_{\varepsilon}^2/\sigma_{\bar{x}}^2)\theta_z^2+1} - corr(x,z)^2}{(\sigma_{\omega}^2/\sigma_{\bar{x}}^2) corr(x,z)^2}} = \frac{\sqrt{\frac{1}{(\sigma_{\varepsilon}^2/\sigma_{\bar{x}}^2)\theta_z^2+1} - corr(x,z)^2}}{(\sigma_{\omega}/\sigma_{\bar{x}}) corr(x,z)}\]

em que na última passagem eu usei o fato do sinal de \(\theta_x\) ser irrelevante, por raciocínio análogo ao feito para o caso anterior.

Da expressão de \(corr(x,z)\) para \(\theta_x = 0\), o que implica que \(x=\bar{x}\), temos que \(corr(\bar{x},z)^2 = \frac{1}{(\sigma_{\varepsilon}^2/\sigma_{\bar{x}}^2)\theta_z^2+1} \), que é exatamente o termo que aparece no numerador da expressão de \(\theta_x\). Esse resultado tem duas implicações interessantes. Em primeiro lugar, vemos que \(\theta_x\) só está bem definido para \(corr(x,z) \leq corr(\bar{x},z)\), o que é natural dado que \(\bar{x}\) é o instrumento ótimo. Em segundo lugar, usando a expressão da \(corr(\varepsilon,z)\), é fácil verificar que \(corr(\bar{x},z)^2 = 1 - corr(\varepsilon,z)^2\). Ou seja, quanto maior a endogeneidade, mais fraco é o instrumento ótimo, o que mostra que endogeneidade e força do instrumento não são totalmente independentes.

Usando esses resultados na expressão de \(\theta_x\), obtemos
\[\theta_x = \frac{\sqrt{1 - corr(\varepsilon,z)^2 - corr(x,z)^2}}{(\sigma_{\omega}/\sigma_{\bar{x}}) corr(x,z)} = \frac{\sqrt{1 - corr(\varepsilon,z)^2 - corr(x,z)^2}}{corr(x,z)}\]

em que na última passagem eu usei \(\sigma_{\omega}=\sigma_{\bar{x}}=1\).

2.3. Estimativas OLS e IV e o método de escolha

Sejam \(Z\), \(X\) e \(Y\) os vetores contendo as \(n\) observações de \(z\), \(x\) e \(y\), respectivamente. Então, sob a hipótese de homocedasticidade, temos os estimadores OLS e IV usuais do parâmetro e da variância assintótica:
\[\hat{\beta}_{OLS} = (Z’Z)^{-1}Z’Y\]

\[\widehat{AVar}(\hat{\beta}_{OLS}) = (Z’Z)^{-1} s_{OLS}^2\]

\[s_{OLS}^2 = \frac{\sum_{i=1}^n \hat{\varepsilon}_{OLS,i}^2}{n-1}\]

\[ \hat{\varepsilon}_{OLS,i} = y_i -\hat{\beta}_{OLS,i}z_i, \forall i\]
\[\hat{\beta}_{IV} = (X’Z)^{-1}X’Y\]

\[\widehat{AVar}(\hat{\beta}_{IV}) = (X’Z)^{-1} (X’X) (Z’X)^{-1} s_{IV}^2\]

\[s_{IV}^2 = \frac{\sum_{i=1}^n \hat{\varepsilon}_{IV,i}^2}{n-1}\]

\[ \hat{\varepsilon}_{IV,i} = y_i -\hat{\beta}_{IV,i}z_i, \forall i\]

O método de escolha do estimador que eu proponho aqui consiste em escolher o estimador que tem o menor MSE assintótico estimado. É fácil mostrar que o \(MSE(\hat{\beta}) = Var(\hat{\beta}) + Bias(\hat{\beta})^2\). Então, para estimar o MSE assintótico, resta obter o viés assintótico. O estimador de IV é não viesado assintoticamente e, logo, o MSE assintótico será igual à variância assintótica. No caso do OLS, existe um viés assintótico. Para calculá-lo, pré-multiplique a equação básica do problema por \(z\) e aplique a esperança:
\[E(zy) = E(zz’) \beta + E(z \varepsilon)\]

\[ \beta = E(zz’)^{-1}E(zy) - E(zz’)^{-1}E(z \varepsilon)\]

O termo \(E(zz’)^{-1}E(zy)\) é o parâmetro da regressão populacional de \(y\) em \(z\), de modo que é possível mostrar que \(plim \hat{\beta}_{OLS} = E(zz’)^{-1}E(zy)\). Logo, o termo \(E(zz’)^{-1}E(z\varepsilon)\), que é o parâmetro da regressão populacional de \(\varepsilon\) em \(z\), é quem determina o viés:
\[ ABias(\hat{\beta}_{OLS}) \equiv plim \hat{\beta}_{OLS} - \beta = E(zz’)^{-1}E(z \varepsilon)\]

Logo, o grande desafio em estimar o viés assintótico é o fato de \(\varepsilon\) não ser observado. A estratégia que eu segui foi estimar essa regressão usando \(\hat{\varepsilon}_{IV}\), se apoiando na consistência do estimador de IV. Assim, sendo \(\hat{\epsilon}_{IV}\) o vetor contendo as \(n\) observações de \(\hat{\varepsilon}_{IV}\),
\[ \widehat{ABias}(\hat{\beta}_{OLS}) = (Z’Z)^{-1}Z’\hat{\epsilon}_{IV}\]

2.4. A medida de MSE

Os estimadores de OLS e IV são assintoticamente normais. Isso implica que os momentos e, logo, o MSE estão bem definidos para amostras grandes. Contudo, para amostras pequenas, tais momentos podem não estar bem definidos, tornando inadequado o cálculo do MSE. Para contornar esse problema, eu optei por calcular uma proxy para o MSE que é passível de ser calculada para qualquer distribuição por usar apenas percentis. Formalmente,
\[ MSE_n(\hat{\beta}) = Var_n(\hat{\beta}) + Bias_n(\hat{\beta})^2\]

\[ Var_n(\hat{\beta}) = \left[ \frac{Q_3(\hat{\beta})-Q_1(\hat{\beta})}{Q_3^{Z}-Q_1^{Z}} \right]^2 = \left[ \frac{Q_3(\hat{\beta})-Q_1(\hat{\beta})}{2Q_3^{Z}} \right]^2\]

\[ Bias_n(\hat{\beta}) = Q_2(\hat{\beta})-\beta \]

em que \(Q_i(\hat{\beta})\) e \(Q_i^Z\) são o i-ésimo quartil de \(\hat{\beta}\) e da normal padrão, respectivamente. Desse modo, \(Q_2(\hat{\beta})\) é a mediana de \(\hat{\beta}\).

Essa medida do MSE tem duas características interessantes. Em primeiro lugar, ela tem um apelo intuitivo, ao usar medidas tradicionais de tendência central e de dispersão quando se trabalha com percentis: mediana e desvio interquartílico, respectivamente. Em segundo lugar, essa medida é exatamente igual ao MSE tradicional caso o estimador seja normalmente distribuído. Logo, essa proxy do MSE é assintoticamente equivalente ao MSE usual no problema aqui considerado.

2.5. Algoritmo de simulação

O algoritmo de simulação consiste nos seguintes passos:

  1. Defina a \(corr(\varepsilon,z) \neq 0\) e calcule \(\theta_z\) como descrito na seção 2.2.
  2. Defina a \(corr(x,\bar{x})>0\) ou a \(corr(x,z)>0\) e calcule \(\theta_x\) como descrito na seção 2.2.
  3. Para cada tamanho de amostra \(n \in \{2,...,80\}\):
    (a) Repita os seguintes passos 100.000 vezes:
    i. Extraia \(n\) observações de \(\bar{x}\), \(\varepsilon\) e \(\omega\).
    ii. Calcule \(y\), \(z\) e \(x\).
    iii. Estime o parâmetro, a variância assintótica, o viés assintótico e o MSE assintótico por OLS e IV, como descrito na seção 2.3.
    (b) Calcule as medidas de variância, viés e MSE dos estimadores apresentadas na seção 2.4 usando os parâmetros estimados em cada iteração. Obtenha também as medianas das estimativas assintóticas dessas estatísticas feitas em cada iteração.

2.6. Experimentos

Eu considerei dois experimentos, que buscam avaliar os estimadores em diferentes circunstâncias. Tais experimentos são utilizados para atender os dois objetivos desse trabalho. O primeiro objetivo é avaliar o desempenho dos estimadores quando variamos o grau de endogeneidade e a força do instrumento. Em particular, eu quero investigar o comportamento do tamanho de amostra mínimo a partir do qual o estimador de IV tem MSE menor que o OLS, que denominei de \(n^*\). O segundo objetivo é testar a proposta de método de escolha do estimador em diferentes contextos.

No primeiro experimento eu mantenho constante a força do instrumento fixando \(corr(x,z)=0,5\) e aumento a intensidade da endogeneidade usando \(corr(\varepsilon,z)= 0,2,\text{ }0,5,\text{ }0,8\). Espera-se que \(n^*\) diminua com o aumento da endogeneidade, dada a inconsistência do OLS nesse ambiente.

No segundo experimento eu mantenho constante o grau de endogeneidade fixando \(corr(\varepsilon,z)=0,2\) e aumento a força do instrumento usando \(corr(x,\bar{x})= 0,6,\text{ }0,7,\text{ }0,8\). Nesse caso, é esperado que \(n^*\) se reduza à medida que o instrumento se torne mais forte, pois é sabido que instrumentos fortes tendem a melhorar o desempenho do estimador de IV.

3. Desempenho dos estimadores

As Figuras 1 e 2 apresentam o desempenho dos estimadores para cada tamanho de amostra nos experimentos 1 e 2, respectivamente. Em cada figura, temos três grupos de gráficos, representando as três simulações efetuadas em cada experimento. A estrutura de cada grupo de gráficos é sempre a mesma: o primeiro gráfico mostra o MSE, o segundo gráfico a variância e o terceiro gráfico o viés ao quadrado. As linhas sólidas representam os valores calculados com base nas 100.000 estimativas do parâmetro para cada tamanho de amostra usando as medidas da seção 2.4, enquanto que os marcadores representam a mediana dos valores assintóticos da seção 2.3 estimados em cada uma das iterações.

De modo geral, os resultados mostram que os valores assintóticos fornecem estimativas bastante precisas para amostras grandes. Para amostras pequenas, o comportamento é bastante aceitável, notadamente para a variância. Além disso, o comportamento de \(n^*\), que é identificado pelo ponto de intersecção das linhas vermelhas e azul nos gráficos do MSE, é o esperado. Afinal, \(n^*\) diminui se a endogeneidade aumenta para um dado nível de força do instrumento (Figura 1) e se o o instrumento se torna mais forte para um dado nível de endogeneidade (Figura 2). Assim, o estimador de IV terá um desempenho relativamente melhor em amostras pequenas quanto maior for a endogeneidade e maior for a força do instrumento, tudo o mais permanecendo constante.

A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
Figura 1: OLS versus IV - aumento marginal do grau de endogeneidade

A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
Figura 2: OLS versus IV - aumento marginal da força do instrumento

4. Desempenho do método de escolha

As Figuras 3 e 4 apresentam o desempenho do método de escolha dos estimadores para cada tamanho de amostra nos experimento 1 e 2, respectivamente. Em cada figura, temos três grupos de gráficos, que têm sempre a mesma estrutura. O primeiro gráfico mostra a porcentagem de acerto do método, isto é, a proporção de vezes em que o método com menor MSE calculado foi de fato o escolhido. O segundo gráfico mostra a proporção de vezes em que o método escolhido teve desempenho pior do que a estratégia econométrica usual de sempre usar o estimador de IV. Ou seja, mostra a proporção de vezes em que o OLS foi escolhido, mas o estimador de IV tinha menor MSE. O terceiro gráfico mostra a proporção de vezes em que o método escolhido teve desempenho melhor do que a estratégia econométrica usual. Isto é, mostra a proporção de vezes em que o OLS foi escolhido e ele era, de fato, o melhor estimador em termos de MSE.

Por fim, o quarto gráfico apresenta o MSE calculado de cada estimador juntamente com o valor esperado do MSE com base no método de escolha, que é a média ponderada do MSE usando como pesos as proporções de escolha do método em cada tamanho de amostra. Além da média, no quarto gráfico eu também apresento uma medida alternativa de tendência central, que é o MSE da moda do método, isto é, o MSE do método mais vezes recomendado em cada tamanho de amostra.

Dos primeiros gráficos de cada grupo, vemos que a porcentagem de acerto do método é aceitável, com os erros se concentrando em períodos em que o MSE calculado de ambos os estimadores são muito próximos. Portanto, o método erra com maior frequência quando esse erro é pouco custoso, mas tende a se comportar bem quando a escolha do estimador é de fato relevante em termos de MSE.

Além disso, a porcentagem de piora frente à estratégia econométrica usual de sempre usar IV na presença de regressores endógenos é, em geral, baixa (segundos gráficos de cada grupo). E nos casos em que essa porcentagem é relevante, vemos que os MSE's de ambos os estimadores são essencialmente iguais. Novamente, os erros do método ocorrem mais frequentemente quando o custo de tais erros é pequeno.

Não apenas esse método não nos leva a escolhas piores do que a estratégia usual, como muitas vezes nos leva a escolhas melhores. Como se vê nos terceiros gráficos de cada grupo, a porcentagem de melhora é substancial, especialmente em amostras pequenas. Isso sugere que o método proposto pode ser bastante útil para problemas reais.

Por fim, nos últimos gráficos de cada grupo, vemos que tanto o valor esperado do MSE quanto a moda do MSE com base no método de escolha são bastante próximos do MSE do estimador mais recomendado em cada tamanho de amostra. Ou seja, o método usualmente indica o uso do melhor estimador.

A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
Figura 3: Método de escolha - aumento marginal do grau de endogeneidade

A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.
Figura 4: Método de escolha - aumento marginal da força do instrumento

0 votos
respondida Out 28 por Thiago Trafane (36 pontos)  

(Continuação)

5. Conclusão

Como esperado, as simulações indicam que o estimador de IV tem desempenho melhor em comparação ao OLS quanto maior for endogeneidade e mais forte for o instrumento, tudo o mais permanecendo constante. Além disso, as simulações mostram resultados promissores para o método de escolha do estimador baseado no MSE assintótico estimado, especialmente em pequenas amostras.

Como qualquer resultado obtido a partir de simulações de Monte Carlo, a principal dúvida é a respeito da generalidade das conclusões. Serão essas conclusões gerais ou fruto dos experimentos e métodos aqui considerados? Essa é uma questão crucial, que pode ser explorada em futuras investigações.

Apêndice - Código do programa em Python

# 0.1) Importando pacotes
import numpy as np
import matplotlib.pyplot as plt
import time
from scipy.stats import norm

# 0.2) Iniciando contagem do tempo de execução
start_time = time.time()

# 1) Definindo funções
# 1.1) Função que faz os plots      
def plot(corr_z_e,corr_x,type_corr_x,n_list,mse_ols,mse_ols_asymp,mse_iv,mse_iv_asymp,var_ols,var_ols_asymp,var_iv,var_iv_asymp,bias2_ols,bias2_ols_asymp,bias2_iv,bias2_iv_asymp,acerto_met,mse_mean,mse_mode,piora_iv,melhora_iv):

        # a) Definindo "base" do título do gráfico e do nome do arquivo
        title_fig = 'corr(e,z) = '+str(corr_z_e)+' e corr(x,'+type_corr_x+') = '+str(corr_x)      
        name_fig = str(corr_z_e).replace('.','')+'_'+str(corr_x).replace('.','')+'_'+type_corr_x     

        # b) Gráfico 1: OLS versus IV
        # Construindo a figura e definindo o título
        fig1,axs1 = plt.subplots(3,sharex=True)        
        axs1[0].set_title(title_fig)          

        # Construindo os gráficos
        y_list = ['mse','var','bias2'] #Nomes das variáveis a serem plotadas
        ylabel_list = ['MSE','VAR','BIAS^2'] #Nomes dos eixos, que identificará as variáveis de cada gráfico

        for i in range(len(y_list)):        
            axs1[i].plot(n_list,locals()[y_list[i]+'_ols_asymp'],'.',color='cornflowerblue',alpha=0.5,markersize=15,label='OLS assint.')
            axs1[i].plot(n_list,locals()[y_list[i]+'_iv_asymp'],'.',color='lightcoral',alpha=0.5,markersize=15,label='IV assint.')
            axs1[i].plot(n_list,locals()[y_list[i]+'_ols'],'-',color='darkblue',linewidth=2,label='OLS')
            axs1[i].plot(n_list,locals()[y_list[i]+'_iv'],'-',color='darkred',linewidth=2,label='IV')   
            axs1[i].set_ylim([0,0.5])    
            axs1[i].set_ylabel(ylabel_list[i])     

        # Alinhando ylabes, adicionando a legenda e salvando a figura        
        axs1[i].legend(loc="upper center", bbox_to_anchor=(0.5, -0.25),ncol=4,frameon=False)   
        fig1.align_ylabels(axs1)
        fig1.savefig('g1_'+name_fig+'.pdf',format='pdf')            
        fig1.savefig('g1_'+name_fig+'.png',format='png')  

        # c) Gráfico 2: método de escolha estimador
        # Construindo a figura e definindo o título
        fig2,axs2 = plt.subplots(4,sharex=True)        
        axs2[0].set_title(title_fig)          

        # Construindo gráfico da porcentagem de acerto
        axs2[0].plot(n_list,acerto_met,'k-',linewidth=2)
        axs2[0].set_ylim([0,100])           
        axs2[0].set_ylabel('% ACERTO')  

        # Construindo gráfico da porcentagem de piora frente IV
        axs2[1].plot(n_list,piora_iv,'k-',linewidth=2)
        axs2[1].set_ylim([0,100])           
        axs2[1].set_ylabel('% INF. IV')  

        # Construindo gráfico da porcentagem de melhora frente IV
        axs2[2].plot(n_list,melhora_iv,'k-',linewidth=2)
        axs2[2].set_ylim([0,100])           
        axs2[2].set_ylabel('% SUP. IV')  

        # Construindo gráfico do MSE    
        axs2[3].plot(n_list,mse_mode,'.',color='gray',alpha=0.5,markersize=15,label='Moda mét.')   
        axs2[3].plot(n_list,mse_ols,'-',color='darkblue',linewidth=2,label='OLS')
        axs2[3].plot(n_list,mse_iv,'-',color='darkred',linewidth=2,label='IV')   
        axs2[3].plot(n_list,mse_mean,':',color='black',linewidth=2,label='Média mét.')          
        axs2[3].set_ylim([0,0.5])           
        axs2[3].set_ylabel('MSE')   
        axs2[3].legend(loc="upper center", bbox_to_anchor=(0.5, -0.3),ncol=4,frameon=False)  

        # Alinhando ylabes e salvando a figura
        fig2.align_ylabels(axs2)           
        fig2.savefig('g2_'+name_fig+'.pdf',format='pdf')            
        fig2.savefig('g2_'+name_fig+'.png',format='png')   

# 1.2) Função que calcula o Erro Quadrático Médio (mse) de cada estimador e mostra os plots
def MSE(n_list,n_simul,beta,corr_z_e,corr_x,type_corr_x,sd_xbar=1,sd_e=1,sd_w=1):

        # a) Obtendo valor do percentil da distribuição normal
        perc = 0.75
        percentil_normal = norm.ppf(perc)

        # b) Calculando theta_z e theta_x
        theta_z = (corr_z_e*sd_xbar/sd_e)/((1-corr_z_e**2)**0.5)   

        if type_corr_x == 'xbar':
            cte_theta_x = 1 
        else:
            cte_theta_x = 1-corr_z_e**2 #Nesse caso, essa constante denota a correlação entre z e xbar           
            if cte_theta_x-corr_x**2 < 0:
                print('Erro: corr(z,x) > corr(z,xbar)')
                return #Interrompendo execução da função        
        theta_x = ((cte_theta_x-corr_x**2)**0.5)/(corr_x*sd_w/sd_xbar)

        # c) Criando as listas a serem preenchidas
        for met in ['ols','iv']:    
            for v in ['var','bias2','mse']:            
                locals()[v+'_'+met+'_list'] = [] 
                locals()[v+'_'+met+'_asymp_list'] = []
        acerto_met_list = []
        mse_mean_list = []
        mse_mode_list = []
        piora_iv_list = []
        melhora_iv_list = []

        # d) Calculando as variáveis de interesse        
        for n in n_list:
            beta_ols_nlist = []
            var_ols_asymp_nlist = []
            bias2_ols_asymp_nlist = []
            mse_ols_asymp_nlist = []
            esc_ols_nlist = []

            beta_iv_nlist = []              
            var_iv_asymp_nlist = []

            # Para cada tamanho de amostra, a regressão é estimada por OLS e IV n_simul vezes
            for i in range(n_simul):            

                # Gerando variáveis aleatórias                
                e = np.resize(np.random.normal(0,sd_e,n),(n,1))            
                xbar = np.resize(np.random.normal(0,sd_xbar,n),(n,1))
                w = np.resize(np.random.normal(0,sd_w,n),(n,1)) 

                # Calculando variáveis a partir das variáveis aleatórias
                z = xbar + theta_z*e
                x = xbar + theta_x*w  
                y = beta*z + e     

                # Comparando correlações teóricas com empíricas (opcional) - rodar com n_list=[10000000] e n_simul=1
                #corr_z_e_emp = np.corrcoef(z[:,0],e[:,0])[0,1]
                #corr_x_emp = np.corrcoef(x[:,0],xbar[:,0])[0,1] if type_corr_x == 'xbar' else np.corrcoef(x[:,0],z[:,0])[0,1]                             
                #print((corr_z_e,round(corr_z_e_emp,2)),(corr_x,round(corr_x_emp,2)))

                # Estimando as regressões
                Szz = z.T@z/n
                Sxx = x.T@x/n
                Sxz = x.T@z/n       
                Szy = z.T@y/n
                Sxy = x.T@y/n

                Szz_inv = np.linalg.inv(Szz)
                Sxz_inv = np.linalg.inv(Sxz)

                beta_ols = Szz_inv@Szy     
                e_ols = y-beta_ols*z
                s2_ols_asymp = sum(e_ols**2)/(n-1)                
                var_ols_asymp = Szz_inv*s2_ols_asymp/n

                beta_iv = Sxz_inv@Sxy  
                e_iv = y-beta_iv*z                             
                s2_iv_asymp = sum(e_iv**2)/(n-1)                
                var_iv_asymp = Sxz_inv@Sxx@Sxz_inv.T*s2_iv_asymp/n     

                bias_ols_asymp = Szz_inv@(z.T@e_iv/n)       
                bias2_ols_asymp = bias_ols_asymp[0,0]**2
                mse_ols_asymp = var_ols_asymp[0,0] + bias2_ols_asymp  

                esc_ols = 1 if mse_ols_asymp <= var_iv_asymp[0,0] else 0                       

                # Adicionando estimativas às listas
                beta_ols_nlist.append(beta_ols[0,0])
                var_ols_asymp_nlist.append(var_ols_asymp[0,0])
                bias2_ols_asymp_nlist.append(bias2_ols_asymp)
                mse_ols_asymp_nlist.append(mse_ols_asymp)
                esc_ols_nlist.append(esc_ols)

                beta_iv_nlist.append(beta_iv[0,0])                
                var_iv_asymp_nlist.append(var_iv_asymp[0,0])

            # Depois que as regressões foram estimadas n_simul vezes, podemos calcular as estatísticas dos estimadores
            for met in ['ols','iv']:
                beta_list = locals()['beta_'+met+'_nlist']     
                var = ((np.quantile(beta_list,perc)-np.quantile(beta_list,1-perc))/(2*percentil_normal))**2
                bias2 = (np.median(beta_list)-beta)**2
                mse = var + bias2                   

                var_asymp = np.median(locals()['var_'+met+'_asymp_nlist'])           
                bias2_asymp = np.median(bias2_ols_asymp_nlist) if met == 'ols' else 0
                mse_asymp = np.median(mse_ols_asymp_nlist) if met == 'ols' else var_asymp                                                    

                # Adicionando estatísticas dos estimadores para o tamanho da amostra atual às listas
                for v in ['var','bias2','mse']:
                    locals()[v+'_'+met+'_list'].append(locals()[v])    
                    locals()[v+'_'+met+'_asymp_list'].append(locals()[v+'_asymp'])    

            # Calculando estatísticas do método de escolha do estimador
            frac_ols = sum(esc_ols_nlist)/n_simul           
            acerto_met = frac_ols*100 if locals()['mse_ols_list'][-1] <= locals()['mse_iv_list'][-1] else (1-frac_ols)*100    
            acerto_met_list.append(acerto_met)   

            mse_mean = frac_ols*locals()['mse_ols_list'][-1] + (1-frac_ols)*locals()['mse_iv_list'][-1]
            mse_mean_list.append(mse_mean)       

            mse_mode = locals()['mse_ols_list'][-1] if frac_ols>=0.5 else locals()['mse_iv_list'][-1]
            mse_mode_list.append(mse_mode)

            piora_iv = frac_ols*100 if locals()['mse_ols_list'][-1] > locals()['mse_iv_list'][-1] else 0
            piora_iv_list.append(piora_iv)    

            melhora_iv = frac_ols*100 if locals()['mse_ols_list'][-1] < locals()['mse_iv_list'][-1] else 0
            melhora_iv_list.append(melhora_iv) 

        # e) Construindo os gráficos usando a função plot            
        plot(corr_z_e,corr_x,type_corr_x,n_list,locals()['mse_ols_list'],locals()['mse_ols_asymp_list'],locals()['mse_iv_list'],locals()['mse_iv_asymp_list'],locals()['var_ols_list'],locals()['var_ols_asymp_list'],locals()['var_iv_list'],locals()['var_iv_asymp_list'],locals()['bias2_ols_list'],locals()['bias2_ols_asymp_list'],locals()['bias2_iv_list'],locals()['bias2_iv_asymp_list'],acerto_met_list,mse_mean_list,mse_mode_list,piora_iv_list,melhora_iv_list)

# 2) Gerando os resultados      
if __name__ == '__main__':    

    # a) Definindo parâmetros gerais das simulações
    n_list = list(range(2,80,1))
    n_simul = 100000
    beta = 5

    # b) Simulação 1: aumento marginal da intensidade da endogeneidade 
    corr_z_e_list = [0.2,0.5,0.8]
    corr_x_z = 0.5
    for corr_z_e in corr_z_e_list:        
        MSE(n_list,n_simul,beta,corr_z_e,corr_x_z,type_corr_x='z')

    # c) Simulação 2: aumento marginal da força do instrumento 
    corr_z_e = 0.2
    corr_x_xbar_list = [0.6,0.7,0.8]    
    for corr_x_xbar in corr_x_xbar_list:        
        MSE(n_list,n_simul,beta,corr_z_e,corr_x_xbar,type_corr_x='xbar')

# 3) Finalizando contagem e exibindo tempo de execução
print("--- %s seconds ---" % (time.time() - start_time))
...