• Para terminar o curso voltemos ao modelo de regressão linear. Porém, desta vez vamos analisar o modelo como uma aproximação para o processo gerador de dados. Com essa abordagem vamos quantificar a incerteza de nossas estimativas.

• O modelo de regressão linear com \(p\) variáveis explicativas (que podem ser chamadas de independentes ou preditoras) é definido por: \[Y_i = \alpha + \beta_1 X_{i1}+\beta_2X_{i2}+\cdots+\beta_pX_{ip} + \varepsilon_i\]
• No modelo temos \(p+1\) coeficientes que devemos estimar.

• De acordo com o modelo a variável de resultado (ou dependente) é gerada como uma função linear das variáveis explicativas e do termo de erro. Tanto a variável de resultado quanto as variáveis explicativas são observadas, apenas o erro não é observado.
• Por isso é fundamental fazer hipóteses a respeito da distribuição do erro.

• A hipótese de exogeneidade para o modelo de regressão linear é definida como: \[E(\varepsilon_i|X_1, X_2, \dots, X_n)=E(\varepsilon_i)=0\]
• Esta hipótese implica que a parte não observada do resultado, contida no termo de erro, não é relacionada com as variáveis explicativas \(X_{ij}\) para \(i=1,2,\dots,n\) e \(j=1,2,\dots,p\).

• De acordo com a hipótese de exogeneidade a epserança do erro condicional às variáveis explicativas é igual a esperança não condicional do erro.
• Desde que o modelo tenha uma constante, a esperança não condicional do erro pode sempre ser suposta igual a zero sem maiores problemas.

• A esperança condicional de uma variável aleatória \(Y\) dada outra variável aleatória \(X\) é denotada por \(E(Y|X)\) e dada por: \[E(Y|X)=\begin{cases}\sum_y y\times f(y|X) \qquad \ \mbox{se }Y\mbox{é discreta}\\ \int y\times f(y|X)dy\qquad \mbox{se }Y\mbox{é contínua}\end{cases}\] Onde \(f(y|X)\) é função de probabilidade condicional (função de densidade condicional) de uma variável aleatória \(Y\) discreta (contínua) dado \(x\).

• Com a hipótese de exogenidade temos que: \[E(Y_i|X_1,X_2,\dots,X_p)= \alpha + \beta_1 X_{i1}+\beta_2 X_{i2} + \cdots + \beta_p X_{ip}\]
• Para observar esse fato lembre que \(E(\beta_j X_{ij}|X_1, X_2, \dots, X_p)=\beta_j X_{ij}\)

• Em estudos aleatórios controlados, a hipótese de exogeneidade é respeitada porque a escolha das unidades tratadas é feita de forma aleatória. Em termos de modelo de regressão linear isso significa que a variável de tratamento, \(X\), é independente de todas as características, observadas ou não, que existem antes do tratamento que estão contidas em \(\varepsilon\).

• Considere o caso de mulheres e políticas públicas que vimos na Unidade 4, a variável explicativa de interesse, \(X\), era se existia reserva para mulheres no Gram Panchayat (GP) da vila. Como a distribuição das reservas foi aleatória, a vila ter a reserva não está relacionado a outras variáveis que podiam influenciar a variável de resultado como tamanho da população ou nível de renda.

• Em estudos com dados observados a validade da hipótese de exogeneidade é bem mais complicada. No exemplo dos GPs se não existisse o programa de reserva de vagas e usássemos a presença de mulheres no GP a exogeneidade poderia ficar comprometida.
• Imagine, por exemplo, que vilas com maior nível de educação tenham uma visão mais liberal e elejam mais mulheres e dêem mais valor a água tratada. Como saber se o maior investimento em equipamentos para taratmento da água foi por conta da presença de mulheres no GP ou pelas características dos eleitores?

• No estudos com dados observados, confundidores não observados (e.g. nível de educação) podem estar contidos no termo de erro, se também estiverem correlacionados com as variáveis explicativas (e.g. gênero dos políticos eleitos) a hipótese de exogeneidade fica comprometida.

• Existem várias técnicas para contornar o problema da não exogeneidade em estudos com dados observados. Uma delas, que discutimos na Unidade 2, é encontar unidades muito parecidas com as tratadas e que não receberam o tratamento.
• Foi a estratégia usada no estudo do salário mínimo em NJ, os autores do estudo usaram restaurantes das mesmas cadeias localizados na vizinha PA como controle para avaliar os efeito do aumento do salário mínimo em NJ.

• Em um modelo de regressão linear, se todos os confundidores estiverem no modelo (e a hipótese de relação linear se aplicar) o coeficiente estimado para a variável de tratamento é um estimador não viesado do efeito médio do tratamento.
• No exemplo do salário mínimo, suponha que os únicos confundidores entre os restaurantes fast-food em NJ e PA sejam a cadeia de fast-food, o salário e a proporção de empregados de tempo integral antes da mudança do salário mínimo em NJ.
• Vamos estimar o modelo com o emprego tempo integral como variável de resultado, estar em NJ como tratamento e as três variáveis acima como confundidores.

library(tidyverse)
minwage <- read.csv("minwage.csv")

minwage <- minwage %>%
  mutate(fullPropBefore = fullBefore/(fullBefore + partBefore),
         fullPropAfter = fullAfter/(fullAfter + partAfter),
         NJ = if_else(location == "PA", 0, 1))

• Para que o R estime os coeficientes de cada cadeia de fast-food vamos retirar a constante do modelo.

fit_minwage <- lm(fullPropAfter ~ -1 + NJ + fullPropBefore + wageBefore + chain, 
                  data = minwage)

fit_minwage

## 
## Call:
## lm(formula = fullPropAfter ~ -1 + NJ + fullPropBefore + wageBefore + 
##     chain, data = minwage)
## 
## Coefficients:
##              NJ   fullPropBefore       wageBefore  chainburgerking  
##         0.05422          0.16879          0.08133         -0.11563  
##        chainkfc        chainroys      chainwendys  
##        -0.15080         -0.20639         -0.22013

• Se deixarmos a constante, o coeficiente dacadeia de referência, Burger King, é retirado e os coeficientes de cada cadeia passam a representar a diferença para o coeficiente do Burger King.

fit_minwage1 <- lm(fullPropAfter ~ NJ + fullPropBefore + wageBefore + chain, 
                   data = minwage)

fit_minwage1

## 
## Call:
## lm(formula = fullPropAfter ~ NJ + fullPropBefore + wageBefore + 
##     chain, data = minwage)
## 
## Coefficients:
##    (Intercept)              NJ  fullPropBefore      wageBefore        chainkfc  
##       -0.11563         0.05422         0.16879         0.08133        -0.03517  
##      chainroys     chainwendys  
##       -0.09076        -0.10451

• Os modelos estimados com e sem a constante são os mesmos, só muda a forma de apresentar os resultados. Repare, por exemplo, que o coeficiente da variável de interesse, \(NJ\), é o mesmo.
• Para ilustrar que os modelos são os mesmos vamos comparar as previsões de cada modelo.

library(modelr)

pred_1 <- minwage %>%
  slice(1) %>% 
  add_predictions(fit_minwage) %>%
  select(pred) %>%
  mutate(model = "fit_minwage")

pred_compare <- minwage %>%
  slice(1) %>%
  add_predictions(fit_minwage1) %>%
  select(pred) %>%
  mutate(model = "fit_minwage1") %>%
  bind_rows(pred_1)

pred_compare

##        pred        model
## 1 0.2709367 fit_minwage1
## 2 0.2709367  fit_minwage

• O uso do modelo linear para fazer inferência necessita da hipótese de exogeneidade. Essa hipótese será violada quando existem confundidores não observados. Para tornar a hipótese de exogeneidade mais plausível, pesquisadores medem variáveis confundidoras e as incluem como variáveis explicativas do modelo de regressão linear.

• Quão precisos são os coeficientes estimados no modelo de regressão linear? Com a hipótese que o modelo linear descreve o verdadeiro processo gerador de dados, podemos quantificar a incerteza associada aos coeficientes estimados.
• Por simplicidade, a análise será feita com a hipótese de apenas uma variável explicativa, mas os resultados valem para modelos com mais de uma variável explicativa.

• Modelo: \[Y_i = \alpha + \beta X_i + \varepsilon_i\]

• Coeficientes estimados \[\begin{align*}\hat{\alpha} &= \bar{Y}-\hat{\beta}\bar{X} \\ \hat{\beta} &= \frac{\sum_{i=1}^n (Y_i-\bar{Y})(X_i - \bar{X})}{\sum_{i=1}^n (X_i - \bar{X})^2}\end{align*}\]

• Com a hipótese de exogeneidade \(\hat{\alpha}\) e \(\hat{\beta}\) são estimadores não viesados de \(\alpha\) e \(\beta\), ou seja, \(E(\hat{\alpha})=\alpha\) e \(E(\hat{\beta})=\beta\).

• Para mostrar esse resultado considere que: \[\begin{align*}\bar{Y} &= \alpha + \beta \bar{X} + \bar{\varepsilon} \Rightarrow \\ \hat{\alpha} &= \alpha + \beta \bar{X}+\bar{\varepsilon}-\hat{\beta}\bar{X}\Rightarrow \\ \hat{\alpha}&=\alpha +(\beta - \hat{\beta})\bar{X}+\bar{\varepsilon}\end{align*}\]
• O erro de estimação, \(\hat{\alpha}-\alpha\) é dado por \((\hat{\beta}-\beta)\bar{X}+\bar{\varepsilon}\)

• Para encontrar \(\hat{\beta}-\beta\):

\[\begin{align*}\hat{\beta}&=\frac{\sum_{i=1}^n (\beta X_i+\varepsilon_i-\beta\bar{X}-\bar{\varepsilon})(X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}\Rightarrow\\ \hat{\beta}&=\frac{\sum_{i=1}^n[\beta(X_i-\bar{X})+(\varepsilon_i-\bar{\varepsilon})](X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}\Rightarrow\\ \hat{\beta}&=\frac{\sum_{i=1}^n \beta (X_i-\bar{X})^2}{\sum_{i=1}^n(X_i - \bar{X})^2}+\frac{\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})(X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}\Rightarrow \\ \hat{\beta} &= \beta + \underbrace{\frac{\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})(X_i-\bar{X})}{\sum_{i=1}^n(X_i - \bar{X})^2}}_{\mbox{erro de estimação}}\end{align*}\]

• O numerador do erro de estimação acima pode ser escrito como:

\[\begin{align*}\sum_{i=1}^n(\varepsilon_i-\bar{\varepsilon})(X_i-\bar{X}) &= \sum_{i=1}^n\varepsilon_i(X_i - \bar{X})+\sum_{i=1}^n\bar{\varepsilon} (X_i - \bar{X})\\ &= \sum_{i=1}^n\varepsilon_i(X_i - \bar{X}) + \bar{\varepsilon}\left(\sum_{i=1}^nX_i - n\bar{X} \right) \\ &= \sum_{i=1}^n\varepsilon_i(X_i - \bar{X})\end{align*}\]

• Desta forma: \[\hat{\beta}-\beta = \frac{\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})}{\sum_{i=1}^n(X_1-\bar{X})^2}\]

• Para seguir com o argumento vamos calcular a esperança condicional de \(\hat{\beta}-\beta\): \[E(\hat{\beta}-\beta|X)=\frac{1}{\sum_{i=1}^n(X_i - \bar{X})^2}\sum_{i=1}^nE(\varepsilon_i|X)(X_i-\bar{X})\]
• Pela hipótese de exogeneidade, \(E(\varepsilon_i|X) = 0\), logo: \[\begin{align*}E(\hat{\beta}-\beta|X) &= 0\\E(\hat{\alpha}-\alpha|X) &= E(\hat{\beta}-\beta|X)\bar{X}+E(\bar{\varepsilon}|X)=0\end{align*}\]
• Pois: \[E(\bar{\varepsilon}|X)=\frac{1}{n}\sum_{i=1}^nE(\varepsilon_i|X)=0\]

• A lei das expectativas iteradas diz que para quaisquer duas variáveis aleatórias \(X\) e \(Y\) vale que: \[E(Y)=E(E(Y|X))\]

• Pela lei das expectativas iteradas:

\[\begin{align*} E(\hat{\alpha})&=E(E(\hat{\alpha}|X))=E(\alpha)=\alpha\\E(\hat{\beta})&=E(E(\hat{\beta}|X))=E(\beta)=\beta\end{align*}\]

• Já temos os valores esperados para os estimadores, precisamos agora do erro padrão.

\[V(\hat{\beta}|X)=V(\hat{\beta}-\beta|X)=V\left(\frac{\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})}{\sum_{i=1}^n(X_i-\bar{X})^2}\bigg|X\right)\\=\frac{1}{\left[\sum_{i=1}^n(X_i - \bar{X})^2\right]^2}V\left(\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})\bigg|X\right)\]

• Para simplificar a expressão acima usamos a hipótese de homocedasticidade, qual seja, condicional em \(X\) os erros \(\varepsilon_i\) são independentes e que a variância do erro não depende de \(X\).

• A hipótese de homocedasticidade dos erros consiste em:
  • \(\varepsilon_i\) é independente de \(\varepsilon_j\) condicional em \(X\) para todo \(i \neq j\);
  • A variância do erro não depende de \(X\), \(V(\varepsilon_i|X)=(\varepsilon_i)\)

• Com a hipótese de homocedasticidade, temos que:

\[V\left(\sum_{i=1}^n\varepsilon_i(X_i - \bar{X})\bigg|X\right)=\sum_{i=1}^n V(\varepsilon_i|X)(X_i-\bar{X})^2\\=V(\varepsilon_i)\sum_{i=1}^n(X_i-\bar{X})^2\] - Portanto: \[V(\hat{\beta}|X)=\frac{V(\varepsilon_i)}{\sum_{i=1}^n(X_i - \bar{X})^2}\]

• A lei da variância total diz que para quaisquer variáveis aleatórias \(X\) e \(Y\) vale que: \[V(Y)=V[E(Y|X)]+E[V(Y|X)]\]

• Aplicando a lei da variância total: \[V(\hat{\beta})=V[E(\hat{\beta}|X)]+E[V(\hat{\beta}|X)]=V(\beta)+E\left(\frac{V(\varepsilon_i)}{\sum_{i=1}^n(X_i - \bar{X})^2}\right)\\=V(\varepsilon_i)E\left(\frac{1}{\sum_{i=1}^n(X_i - \bar{X})^2}\right)\]

• Repare que o variância não condicional de \(\hat{\beta}\) é igual ao valor esperado da variância condicional de \(\hat{\beta}\), por conta disso um bom estimador da variância condicional também é um bom estimador da variância não condicional.

• Com esse resultado podemos calcular o erro padrão de \(\hat{\beta}\) com a raiz do estimador de variância. \[\mbox{erro padrão de }\hat{\beta}=\sqrt{\widehat{V(\hat{\beta})}}=\sqrt{\frac{\frac{1}{n}\sum_{i=1}^n \hat{\varepsilon}_i^2}{\sum_{i=1}^n(X_i - \bar{X})^2}}\]
• Para obter esse resultado usamos que média amostral dos resíduos é zero, de forma que: \(\frac{1}{n}\sum_{i=1}^n(\hat{\varepsilon}_i-\bar{\hat{\varepsilon}})^2=\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon}_i^2\).

• Sem a hipótese de homocedasticidade a fórmula acima não é correta para estimar o erro padrão do estimador.
• Não é objetivo desse curso, mas o R tem pacotes com funções que calculam o erro padrão robusto para heterocedasticidade.

• Com o estimador para o erro padrão podemos construir intervalos de confiança para \(\hat{\beta}\). Usando o teorema do limite central temos que: \[\mbox{z-score de }\hat{\beta}=\frac{\hat{\beta}-\beta}{\mbox{erro padrão de }\hat{\beta}}\sim N(0,1)\]

• O intervalo com nível de confiança \((1-\alpha)\times 100\%\) é dado por:

\[\mbox{CI}(\alpha)=\left[\hat{\beta}-z_{\frac{\alpha}{2}}\times\mbox{erro padrão de }\hat{\beta}, \hat{\beta}+z_{\frac{\alpha}{2}}\times\mbox{erro padrão de }\hat{\beta}\right]\]

• Também podemos fazer testes de hipóteses a respeito dos coeficientes, por exemplo, podemos testar se o coeficiente estimado é igual a \(\beta^\star\), na maioria das vezes queremos testar se o coeficiente é igual a zero.

• Com as hipótese que estamos utilizando a estatística de teste para a hipótese nula \(\beta = \beta^\star\) é dada por \(z^\star = \frac{\hat{\beta}-\beta^\star}{\mbox{erro padrão de }\hat{\beta}}\). A distribuição de \(z^\star\) sob a hipótese nula é uma normal padrão.
• Desta forma, podemos calcular o p-valor usando a função de distribuição acumulada da normal padrão.

• Assim como na seção anterior, é comum usar a distribuição \(t\) para fazer testes de hipóteses sobre os coeficientes estimados.
• Com as hipóteses que o termo de erro tem distribuição normal e é homocedastico a estatística \(z^\star\) tem distribuição \(t\) com \(n-2\) graus de liberdade. No caso é \(n-2\) porque estimamos dois parâmetros, \(\alpha\) e \(\beta\).

• Como primeiro exemplo vamos usar o caso das mulheres nos GPs na Índia.

women <- read.csv("women.csv")

fit_women <- lm(water ~ reserved, data= women)

• A função summary() retorna os valores estimados e as estatísticas relevantes.
• A função tidy() do pacote broom retornam essas informações em um tibble.

summary(fit_women)

## 
## Call:
## lm(formula = water ~ reserved, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.991 -14.738  -7.865   2.262 316.009 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   14.738      2.286   6.446 4.22e-10 ***
## reserved       9.252      3.948   2.344   0.0197 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.45 on 320 degrees of freedom
## Multiple R-squared:  0.01688,    Adjusted R-squared:  0.0138 
## F-statistic: 5.493 on 1 and 320 DF,  p-value: 0.0197

library(broom)
tidy(fit_women)

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    14.7       2.29      6.45 4.22e-10
## 2 reserved        9.25      3.95      2.34 1.97e- 2

• A função tidy() permite calcular os intervalos de confiança

tidy(fit_women, conf.int=TRUE)

## # A tibble: 2 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)    14.7       2.29      6.45 4.22e-10    10.2       19.2
## 2 reserved        9.25      3.95      2.34 1.97e- 2     1.49      17.0

• Intervalo de confiança a 99%

tidy(fit_women, conf.int=TRUE, conf.level=0.99)

## # A tibble: 2 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)    14.7       2.29      6.45 4.22e-10    8.81       20.7
## 2 reserved        9.25      3.95      2.34 1.97e- 2   -0.977      19.5

• Repare que esse intervalo de confiança contém o zero.

• Como foi dito, os resultado obtidos com uma variável explicativa valem para modelos com mais variáveis explicativas. Não é objetivo desse curso demonstrar esses resultados em modelos com muitas variáveis (para isso tem a sequência de econometria), mas vamos usar o estudo do salário mínimo como exemplo de regressão múltipla no R.

tidy(fit_minwage, conf.int=TRUE)

## # A tibble: 7 × 7
##   term            estimate std.error statistic p.value conf.low conf.high
##   <chr>              <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 NJ                0.0542    0.0332     1.63  0.103   -0.0111      0.120
## 2 fullPropBefore    0.169     0.0566     2.98  0.00307  0.0574      0.280
## 3 wageBefore        0.0813    0.0389     2.09  0.0374   0.00478     0.158
## 4 chainburgerking  -0.116     0.179     -0.646 0.518   -0.467       0.236
## 5 chainkfc         -0.151     0.183     -0.824 0.411   -0.511       0.209
## 6 chainroys        -0.206     0.187     -1.11  0.270   -0.574       0.161
## 7 chainwendys      -0.220     0.188     -1.17  0.243   -0.591       0.150

• O resultado da regressão é que não podemos rejeitar a hipótese nula que o coeficiente de \(NJ\) é zero, ou seja, o restaurante estar em NJ, onde houve aumento do mínimo, não parece ser relevante para explciar a proporção de empregados em tempo integral.

Referência: A contribution to the empirics of economic growth; Gregory Mankiw, David Romer e David Weil, QJE, 1992

• “This paper takes Robert Solow seriously.”
• Os autores argumentam que os dados suportam a relação da taxa de poupança e da taxa de crescimento da população com nível de renda que é obtida no Modelo de Solow, o modelo também explica parte considerável da variação entre as rendas per capitas dos diversos países.
• Uma versão ampliada do modelo para incluir capital humano “corrige” encontra melhor suporte nos dados.

• A introdução de capital humano têm dois efeitos relevantes:
• Para uma dada taxa de acumulação de capital humano, maiores taxas de poupança ou menores taxas de crescimento populacional levam a um maior nível de renda e, portanto, a um maior nível de capital humano; desta forma, a acumulação de capital físico e taxa de crescimento da população têm maior impacto na renda quando a acumulação de capital humano é considerada.
• Acumulação de capital humano pode estar correlacionada com taxa de poupança e crescimento da população; a omissão da acumulação de capital humano colocam um viés nas estimativas dos coeficientes da taxa de poupança e da taxa de crescimento da população.

• Equação estimada: \[ \ln\left( \frac{Y}{L}\right) = a + \frac{\alpha}{1-\alpha} \ln(s) - \frac{\alpha}{1-\alpha} \ln(g + n + \delta) + \varepsilon \]

• Hipótese: as taxas de poupança e de crescimento da população são independentes de fatores específicos a países que deslocam a função de produção, ou seja, \(s\) e \(n\) são independentes de \(\varepsilon\).

• Como os fatores são remunerados por seus produtos marginais, \(\alpha\) é fração da renda destinada ao capital.
• Estimativas sugerem que \(\alpha = 1/3\), nesse caso a previsãpo do modelo é que o coeficiente de \(\ln(s)\) é 0,5 e o coeficiente de \(\ln(n+g+\delta)\) é -0,5.

• Três amostras:

• A primeira inclui todos os países da PWT com dados disponíveis excluindo os países que tem a produção de petróleo como principal indústria, é formada por 98 países.

• A segunda exclui os países com dados classificados como {} na PWT, também foram excluídos os países com menos de um milhão de habitantes em 1960, é formada por 75 países.

• A terceira consiste mnos 22 países da OCDE com mais de um milhão de habitantes.

library(tidyverse)
library(foreign)
library(car)
library(lmtest)
library(stargazer)

mrw <- read.dta("mrw.dta")
head(mrw)

##   number      country n i o rgdpw60 rgdpw85 gdpgrowth popgrowth  i_y school
## 1      1      Algeria 1 1 0    2485    4371       4.8       2.6 24.1    4.5
## 2      2       Angola 1 0 0    1588    1171       0.8       2.1  5.8    1.8
## 3      3        Benin 1 0 0    1116    1071       2.2       2.4 10.8    1.8
## 4      4     Botswana 1 1 0     959    3671       8.6       3.2 28.3    2.9
## 5      5 Burkina Faso 1 0 0     529     857       2.9       0.9 12.7    0.4
## 6      6      Burundi 1 0 0     755     663       1.2       1.7  5.1    0.4

• Variáveis:
• number: número que identifica cada um dos 121 países
• country: nome do país
• n: dummy com valor um se petróleo não é a indústria dominante no país.
• i: dummy com valor um se o país foi incluído na amostra intermediária.
• o: dummy com valor um se o país faz parte da OCDE
• rgdpw60: PIB real por população em idade ativa no ano de 1960
• rgdpw85: PIB real por população em idade ativa no ano de 1985
• popgrowth: taxa média de crescimento da população entre 1960 e 1985
• i_y85: taxa de investimento, 1960 a 1985
• school: uma proxy para capital humano, a construção é explicada na página 419 de Mankiw, Romer e Weil (1992).

mrw <- mrw %>%
  mutate(rgpdw60 = as.numeric(rgdpw60),
         rgdpw85 = as.numeric(rgdpw85)) %>%
  mutate(lrgdpw60 = log(rgdpw60),
         lrgdpw85 = log(rgdpw85),
         ly_i = log(i_y/100),
         lngd = log(0.05+popgrowth/100),
         rdif = ly_i - lngd,
         lh = log(school/100),
         growth = log(rgdpw85/rgdpw60))

reg1 <- lm(lrgdpw85~ly_i + lngd, data=mrw, subset = n==1)
summary(reg1)

## 
## Call:
## lm(formula = lrgdpw85 ~ ly_i + lngd, data = mrw, subset = n == 
##     1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.79144 -0.39367  0.04124  0.43368  1.58046 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.4299     1.5839   3.428 0.000900 ***
## ly_i          1.4240     0.1431   9.951  < 2e-16 ***
## lngd         -1.9898     0.5634  -3.532 0.000639 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6891 on 95 degrees of freedom
## Multiple R-squared:  0.6009, Adjusted R-squared:  0.5925 
## F-statistic: 71.51 on 2 and 95 DF,  p-value: < 2.2e-16

reg2 <- lm(lrgdpw85~ly_i + lngd, data=mrw, subset = i==1)
reg3 <- lm(lrgdpw85~ly_i + lngd, data=mrw, subset = o==1)

stargazer(reg1,reg2,reg3, type="html", single.row = FALSE)

• O modelo indica que os coeficentes de ly_i e lngd possuem o mesmo módulo e sinais diferentes. Vamos testar essa hipótese.
• Esse tipo de teste envolve calcular uma estaística com distribuição F usando o \(R^2\) dos modelos com e sem restrição e os graus de liberdade das regressões.

reg1r <- lm(lrgdpw85~ rdif, data=mrw, subset = n==1)
reg2r <- lm(lrgdpw85~rdif, data=mrw, subset = i==1)
reg3r <- lm(lrgdpw85~rdif, data=mrw, subset = o==1)

linearHypothesis(reg1, "ly_i + lngd = 0")

## Linear hypothesis test
## 
## Hypothesis:
## ly_i  + lngd = 0
## 
## Model 1: restricted model
## Model 2: lrgdpw85 ~ ly_i + lngd
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     96 45.504                           
## 2     95 45.108  1   0.39612 0.8343 0.3634

linearHypothesis(reg2, "ly_i = -lngd")

## Linear hypothesis test
## 
## Hypothesis:
## ly_i  + lngd = 0
## 
## Model 1: restricted model
## Model 2: lrgdpw85 ~ ly_i + lngd
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     73 27.330                           
## 2     72 26.848  1   0.48226 1.2933 0.2592

linearHypothesis(reg3, "ly_i + lngd = 0")

## Linear hypothesis test
## 
## Hypothesis:
## ly_i  + lngd = 0
## 
## Model 1: restricted model
## Model 2: lrgdpw85 ~ ly_i + lngd
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     20 2.7147                           
## 2     19 2.7061  1 0.0085943 0.0603 0.8086

• Resultados favoráveis ao Modelo de Solow:
• Os coeficientes das taxas de poupança e crescimento da população possuem os sinais esperados e, em duas das três amostras, são significativos.
• A hipótese de que os coeficientes são iguais em valores absolutos e possuem sinais opostos não é rejeitada em nenhuma das regressões.
• O modelo explica boa parte da variação entre as rendas dos países (R2 em torno de 60% nas duas regressões com coeficientes significativos).

• Resultados desfavoráveis ao Modelo de Solow:
• Os valores de \(\alpha\) ímplicito nos coeficientes estimados (e significativos) é muito alto, próximo de 60% quando os valores normalmente estimados são próximos a 33%.

• A solução apresentada pelos autores foi incluir capital humano no modelo de Solow. Equação estimada passa a ser:

\[ \ln\left( \frac{Y}{L}\right) = a + \frac{\alpha}{1-\alpha - \beta} \ln(s_k) - \frac{\alpha + \beta}{1-\alpha-\beta} \ln(g + n + \delta) \\ +\frac{\beta}{1-\alpha-\beta} \ln(h^\star) + \varepsilon \]

mrw <- mrw %>% mutate(hdif = lh - lngd)

reg1.h <- lm(lrgdpw85~ly_i + lngd + lh, data=mrw, subset = n==1)
reg2.h <- lm(lrgdpw85~ly_i + lngd + lh, data=mrw, subset = i==1)
reg3.h <- lm(lrgdpw85~ly_i + lngd + lh, data=mrw, subset = o==1)

reg1r.h <- lm(lrgdpw85~ rdif + hdif, data=mrw, subset = n==1)
reg2r.h <- lm(lrgdpw85~rdif + hdif, data=mrw, subset = i==1)
reg3r.h <- lm(lrgdpw85~rdif + hdif, data=mrw, subset = o==1)

linearHypothesis(reg1.h, "ly_i + lngd  + lh= 0")

## Linear hypothesis test
## 
## Hypothesis:
## ly_i  + lngd  + lh = 0
## 
## Model 1: restricted model
## Model 2: lrgdpw85 ~ ly_i + lngd + lh
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     95 24.418                           
## 2     94 24.226  1   0.19185 0.7444 0.3904

linearHypothesis(reg2.h, "ly_i + lh = -lngd")

## Linear hypothesis test
## 
## Hypothesis:
## ly_i  + lngd  + lh = 0
## 
## Model 1: restricted model
## Model 2: lrgdpw85 ~ ly_i + lngd + lh
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     72 14.684                           
## 2     71 14.680  1 0.0045263 0.0219 0.8828

linearHypothesis(reg3.h, "ly_i + lngd + lh = 0")

## Linear hypothesis test
## 
## Hypothesis:
## ly_i  + lngd  + lh = 0
## 
## Model 1: restricted model
## Model 2: lrgdpw85 ~ ly_i + lngd + lh
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     19 1.9604                           
## 2     18 1.9603  1 0.0001471 0.0014 0.9711

• Conclusões:
• A medida de capital humano é significativa nas três amostras.
• A introdução de capital humano reduz os coeficientes do capital físico.
• Nas duas maiores amostras o modelo explica cerca de 80% da variação do PIB por trabalhador.
• O modelo ampliado para incluir capital humano prevê que a soma dos coeficientes das três variáveis é igual a zero, esse resultado não foi rejeitado em nenhuma das três amostras.
• Nas duas maiores amostras o valor de \(\alpha\) foi significativo e compativel com os valores estimados para a participação da renda do capital na renda total.
• A partir desses resultados os autores concluem que a introdução de capital humano melhora a performance do modelo.

• Na sequência os autores discutem a hipótese de convergência, qual seja, que maiores taxas de crescimento estão relacionadas a menores valores do PIB inicial do PIB por trabalhador.
• Regressão apenas com PIB por trabalhador inicial
• Regressão com PIB por trabalhador inicial, taxa de poupança e lngd
• Regressão com PIB por trabalhador inicial, taxa de poupança, lngd e capital humano

reg1.III.conv <- lm(growth~lrgdpw60, data=mrw, subset = n==1)
reg2.III.conv <- lm(growth~lrgdpw60, data=mrw, subset = i==1)
reg3.III.conv <- lm(growth~lrgdpw60, data=mrw, subset = o==1)

reg1.IV.conv <- lm(growth~lrgdpw60 +ly_i + lngd, data=mrw, subset = n==1)
reg2.IV.conv <- lm(growth~lrgdpw60 +ly_i + lngd, data=mrw, subset = i==1)
reg3.IV.conv <- lm(growth~lrgdpw60 +ly_i + lngd, data=mrw, subset = o==1)

reg1.V.conv <- lm(growth~lrgdpw60 +ly_i + lngd + lh, data=mrw, subset = n==1)
reg2.V.conv <- lm(growth~lrgdpw60 +ly_i + lngd + lh, data=mrw, subset = i==1)
reg3.V.conv <- lm(growth~lrgdpw60 +ly_i + lngd + lh, data=mrw, subset = o==1)

• Após a introdução de outras variáveis explicativas além do PIB por trabalhador inicial faz com que o coeficiente do PIB por trabalhador inicial seja negativo e significativo como previsto pela hipótese de convergência.

• Vamos aproveitar esse exemplo para falar de heterocedasticidade.

mrw.h <- mrw %>% filter(n==1)
reg <- lm(growth~ lrgdpw60 + ly_i + lngd, data=mrw.h)

## 
## Call:
## lm(formula = growth ~ lrgdpw60 + ly_i + lngd, data = mrw.h)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07648 -0.15215  0.01185  0.19595  0.96056 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.91938    0.83367   2.302  0.02352 *  
## lrgdpw60    -0.14090    0.05202  -2.709  0.00803 ** 
## ly_i         0.64724    0.08670   7.465 4.16e-11 ***
## lngd        -0.30235    0.30438  -0.993  0.32311    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3507 on 94 degrees of freedom
## Multiple R-squared:  0.4019, Adjusted R-squared:  0.3828 
## F-statistic: 21.05 on 3 and 94 DF,  p-value: 1.622e-10

• Um gráfico entre os resíduos dessa regressão e a variável lngd sugere que quanto maior lngd, maior será a variação dos resíduos.

ggplot(mrw.h, aes(lngd, resid(reg))) + geom_point()

• O teste de Breusch-Pagan testa a hipótese nula de homocedasticidade.
• p.valor baixo -> rejeita a hipótese nula
• A função bptest, do pacote, lmtest implementa o teste de Breusch-Pagan

bptest(reg)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg
## BP = 9.1463, df = 3, p-value = 0.02741

• A função coeftest, estima o erro padrão robusto para heterocedasticidade.
• O argumento vcov especifica o tipo de correção. Se o argumento não for definido a função retorna os valores sem correção, se for definido como vcov=hccm(reg, type=“hc0”) faz a correção como proposta em White (1980).

coeftest(reg)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  1.919380   0.833672  2.3023  0.023524 *  
## lrgdpw60    -0.140901   0.052018 -2.7087  0.008027 ** 
## ly_i         0.647238   0.086700  7.4653 4.159e-11 ***
## lngd        -0.302346   0.304383 -0.9933  0.323110    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(reg, vcov=hccm(reg, type="hc0"))

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  1.919380   0.804162  2.3868  0.018998 *  
## lrgdpw60    -0.140901   0.048231 -2.9214  0.004363 ** 
## ly_i         0.647238   0.101005  6.4080 5.819e-09 ***
## lngd        -0.302346   0.241106 -1.2540  0.212953    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Seja o modelo: \[Y_i = \beta_0 + \beta_1 X_i + u_i\]

• Suponha que o termo de erro, \(u_i\), seja correlacionado com \(X_i\), ou seja, a hipótese de exogeneidade não é observada.

• O estimador de mínimos quadrados ordinário não é consistente.

• Em casos assim, o estimador de variáveis instrumentais (IV), permite obter estimativas consistentes de \(\beta_1\).
• Para ser um instrumento válido a variável \(Z\) deve atender duas condições:
• X e Z devem ser correlacionados
• Z não pode ser correlacionado com o termo de erro

• Estimação em dois estágios:
• Primeiro estagio: a variação do regressor, X, é decomposta no termo explicado por Z (sem correlação com u) e no termo não explicado por Z (tem correlação com u). \[X_i = \pi_0 + \pi_1 Z_i + v_i\]
• Salve os valores estimados de X na regressão acima: \[\hat{X}_i = \hat{\pi}_0 + \hat{\pi}_1 Z_i\]
• Se \(Z\) é um instrumento válido, então \(\hat{X}\) é exógeno em relação a \(Y\).
• O segundo estágio consiste em uma regressão entre \(Y\) e \(\hat{X}\).

• Para o caso de apenas uma variável explicativa é possível mostrar que o valor estimado para \(\beta_1\) no segundo estágio é dado por: \[ \hat{\beta}_1 = \frac{Cov(Z,Y)}{Cov(Z,X)} \]

• Exemplo: demanda por cigarros
• O objetivo é estimar \(\beta_1\) em uma regressão entre a quantidade vendida e o preço do pacote de cigarros. \[\log(pack) = \beta_0 + \beta_1 \log(price) + u_i\]

library(AER)
data("CigarettesSW")

head(CigarettesSW)

##   state year   cpi population    packs    income  tax     price     taxs
## 1    AL 1985 1.076    3973000 116.4863  46014968 32.5 102.18167 33.34834
## 2    AR 1985 1.076    2327000 128.5346  26210736 37.0 101.47500 37.00000
## 3    AZ 1985 1.076    3184000 104.5226  43956936 31.0 108.57875 36.17042
## 4    CA 1985 1.076   26444000 100.3630 447102816 26.0 107.83734 32.10400
## 5    CO 1985 1.076    3209000 112.9635  49466672 31.0  94.26666 31.00000
## 6    CT 1985 1.076    3201000 109.2784  60063368 42.0 128.02499 51.48333

• instrumento: sales tax
• hipótese: essa taxa só influencia a quantidade por meio do preço

dados <- CigarettesSW %>%
  mutate(rprice = price/cpi,    #prego real
         salestax = (taxs - tax)/cpi  #sales tax
  ) %>%
  filter(year == 1995)

cor.test(dados$salestax, dados$rprice)

## 
##  Pearson's product-moment correlation
## 
## data:  dados$salestax and dados$rprice
## t = 6.3877, df = 46, p-value = 7.578e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4986121 0.8116363
## sample estimates:
##       cor 
## 0.6856138

• Primeiro estágio

reg1.cig <- lm(log(rprice) ~ salestax, data = dados)

• Primeiro estágio

coeftest(reg1.cig, vcov = vcovHC)

## 
## t test of coefficients:
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 4.616546   0.030847 149.6586 < 2.2e-16 ***
## salestax    0.030729   0.005142   5.9761 3.145e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Quanto da variação do log de rprice é explicado por saletax?

summary(reg1.cig)$r.squared

## [1] 0.4709961

• Guarde os valores estimados (fitted) no primeiro estágio

dados$lrprice.pred <- reg1.cig$fitted.values

• Segundo estágio

reg2.cig <- lm(log(packs) ~ lrprice.pred, data = dados)
coeftest(reg2.cig, vcov = vcovHC)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   9.71988    1.70304  5.7074 7.932e-07 ***
## lrprice.pred -1.08359    0.35563 -3.0469  0.003822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Usando a fórmula:

cov(dados$salestax,log(dados$packs))/cov(dados$salestax,log(dados$rprice))

## [1] -1.083587

• A função ivreg do pacote AER faz a estimação por variáveis instrumentais

reg_iv1.cig <- ivreg(log(packs) ~ log(rprice) | salestax, data = dados)
coeftest(reg_iv1.cig, vcov = vcovHC)

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  9.71988    1.60583  6.0529 2.413e-07 ***
## log(rprice) -1.08359    0.33508 -3.2338  0.002263 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Regressões múltiplas com variáveis instrumentais: \[Y = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k + \beta_{k+1} W_1 + \cdots + \beta_{k+r} W_r + u\]

• X: regressores endógenos (1,…,k)

• W: regressores exógenos (1,…,r)

• Z: variáveis instrumentais (1,…,m)

• Para usar o estimador VI precisamos que \(m \geq k\)

• Estimação por dois estágios:
• Faça regressão por OLS de cada variavel explicativa endógena (X_1,…,X_k) em todas as variaveis instrumentais (Z_1,…,Z_m), todas as variaveis exógenas (W_1,…,W_r) e um intercepto. Guarde os valores estimados para cada X.
• Faça a regressão da variável dependente nos valores estimados para X, nas variaveis exógenas e no intercepto usando OLS.

• Como exemplo considere a inclusão da renda per capita de cada estado como variável explicativa exógena.
• Vamos usar a função ivreg.

dados <- dados %>%
  mutate(rincome = income/population/cpi)   #renda real per capita

log(packs): variável dependente log(rprice): variável explicativa endógena log(rincome): variável explicativa exógena salestax: instrumento

reg_iv2.cig <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + 
                      salestax, data = dados)
coeftest(reg_iv2.cig, vcov = vcovHC)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   9.43066    1.34793  6.9964 1.032e-08 ***
## log(rprice)  -1.14338    0.40109 -2.8507  0.006562 ** 
## log(rincome)  0.21452    0.33067  0.6487  0.519807    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Para mais um exemplo vamos adicionar tax como segundo instrumento

dados <- dados %>%
  mutate(cigtax = tax/cpi)

reg_iv3.cig <- ivreg(log(packs) ~ log(rprice) + log(rincome) | 
                     log(rincome) + salestax + cigtax, data = dados)
coeftest(reg_iv3.cig, vcov = vcovHC)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   9.89496    1.02539  9.6500 1.569e-12 ***
## log(rprice)  -1.27742    0.26704 -4.7837 1.883e-05 ***
## log(rincome)  0.28040    0.26408  1.0618     0.294    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• O próximo exemplo diz respeito a estimação dos retornos da educação usando a distância para faculdade como variável instrumetal.
• Referência: Using Geographic Variation in College Proximity to Estimate the Return to Schooling, David Card (1993), NBER.
• Objetivo: estimar o impacto da educação no salário.

data(CollegeDistance)
head(CollegeDistance)

##   gender ethnicity score fcollege mcollege home urban unemp wage distance
## 1   male     other 39.15      yes       no  yes   yes   6.2 8.09      0.2
## 2 female     other 48.87       no       no  yes   yes   6.2 8.09      0.2
## 3   male     other 48.74       no       no  yes   yes   6.2 8.09      0.2
## 4   male      afam 40.40       no       no  yes   yes   6.2 8.09      0.2
## 5 female     other 40.48       no       no   no   yes   5.6 8.09      0.4
## 6   male     other 54.71       no       no  yes   yes   5.6 8.09      0.4
##   tuition education income region
## 1 0.88915        12   high  other
## 2 0.88915        12    low  other
## 3 0.88915        12    low  other
## 4 0.88915        12    low  other
## 5 0.88915        13    low  other
## 6 0.88915        12    low  other

reg1.wage <- lm(log(wage) ~ education, data = CollegeDistance)
reg2.wage <- lm(log(wage) ~ unemp + ethnicity + gender + urban + education, 
                data = CollegeDistance)

coeftest(reg1.wage, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 2.2132204  0.0163823 135.0982  < 2e-16 ***
## education   0.0020244  0.0011731   1.7257  0.08447 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(reg2.wage, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##                      Estimate  Std. Error  t value Pr(>|t|)    
## (Intercept)        2.15199991  0.01704415 126.2603  < 2e-16 ***
## unemp              0.01359382  0.00073822  18.4143  < 2e-16 ***
## ethnicityafam     -0.06191387  0.00611007 -10.1331  < 2e-16 ***
## ethnicityhispanic -0.05352044  0.00466547 -11.4716  < 2e-16 ***
## genderfemale      -0.00911503  0.00397305  -2.2942  0.02182 *  
## urbanyes           0.00893926  0.00466207   1.9174  0.05524 .  
## education          0.00067227  0.00111739   0.6016  0.54744    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Problema: educação não é distribuída de forma aleatória entre os indivíduos, trata-se de uma escolha. Assim é possível argumentar que educação não é exógena para explicar a renda, é possível, por exemplo, que indivíduos mais talentosos que a média escolham maiores níveis de educação e, por conta do talento, tenham renda acima da média.
• Hipótese: distância da faculdade afeta a decisão de obter mais educação e afeta a renda apenas por meio da educação.

cor.test(CollegeDistance$distance, CollegeDistance$education)

## 
##  Pearson's product-moment correlation
## 
## data:  CollegeDistance$distance and CollegeDistance$education
## t = -6.4414, df = 4737, p-value = 1.301e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.12133362 -0.06488278
## sample estimates:
##         cor 
## -0.09318309

reg.wage.iv2 <- ivreg(log(wage) ~ unemp + ethnicity + gender + urban + 
                        education | . - education + distance, data = CollegeDistance)
coeftest(reg.wage.iv2, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##                     Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        1.2171787  0.2008632  6.0597 1.469e-09 ***
## unemp              0.0142234  0.0009643 14.7500 < 2.2e-16 ***
## ethnicityafam     -0.0277621  0.0101397 -2.7380  0.006205 ** 
## ethnicityhispanic -0.0335043  0.0080922 -4.1403 3.529e-05 ***
## genderfemale      -0.0076101  0.0052689 -1.4444  0.148706    
## urbanyes           0.0064494  0.0063425  1.0168  0.309277    
## education          0.0673242  0.0143308  4.6979 2.703e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Quais as causas fundamentais para explicar as grandes diferenças entre as renda per capitas dos diversos países?
• Referência: The Colonial Origins of Comparative Development: An Empirical Investigation; Daron Acemoglu, Simon Johnson e James Robinson, The American Economic Review, Vol. 91, No. 5 (Dec., 2001), pp. 1369-1401.

• Instituições e direitos de propriedade são explicações promissoras.
• “Countries with better institutions, more secure property rights, and less distortionary policies will invest more in physical and human capital, and will use these factors more efficiently to achieve a greater level of income.”

• Problema: possível endogeneidade
• “It is quite likely that rich economies choose or can afford better institutions. Perhaps more important, economies that are different for a variety of reasons will differ both in their institutions and in their income per capita.”
• “To estimate the impact of institutions on economic performance, we need a source of exogenous variation in institutions. In this paper, we propose a theory of institutional differences among countries colonized by Europeans, and exploit this theory to derive a possible source of exogenous variation.”

• Premissas:
• Diferentes tipos de colonização criaram diferentes conjuntos de instituições. Em um extremo estão os estados extrativistas: o objetivo era extrair da colônia o máximo de recursos possível; não havia incentivos para criar instituições para proteção da propriedade privada e imposição de limites (checks and balances) para o governo. No outro extremo estão os casos onde os colonizadores tentavam reproduzir as instiuições da Europa, criação de Neo-Europas. Havia incentivos para criar instituições para proteger a propriedade privada e limitar a ação do governo.
• A estratégia de colonização depende do ambiente, em locais hostis não exisitiam incentivos para criar as Neo-Europas.
• As instituições coloniais persistem mesmo após a independência.

• Com essas premissas os autores usam a taxa de mortalidade esperada nas colônias como instrumento para as atuais instituições das ex-colônias.
• taxa de mortalidade potencial => colonização => instituições iniciais => instituições atuais => renda per capita atual

data2 <- read.dta("maketable2.dta")
data2.r <- data2 %>% filter(baseco == 1)

• Estimativas por mínimos quadrados

reg2.2 <- lm(logpgp95 ~ avexpr, data = data2.r)
reg2.5 <- lm(logpgp95 ~ avexpr + lat_abst, data = data2)

• Relação positiva e significativa entre instituições e PIB per capita

• “Overall, the results in Table 2 show a strong correlation between institutions and economic performance. Nevertheless, there are a number of important reasons for not interpreting this relationship as causal. First, rich economies may be able to afford, or perhaps prefer, better institutions. Arguably more important than this reverse causality problem, there are many omitted determinants of income differences that will naturally be correlated with institutions. Finally, the measures of institutions are constructed ex-post, and the analysts may have had a natural bias in seeing better institutions in richer places.”

• Após uma série de argumentos (e regressões) para defender a proposta de identificação (taxa de mortalidade potencial => colonização => instituições iniciais => instituições atuais => renda per capita atual), os autores refazem as estimações usando mortalidade nas colônias como variável instrumental.

data4 <- read.dta("maketable4.dta")
data4.r <- data4 %>% filter(baseco == 1)

• Estimativas por variáveis instrumentais

reg4.A1 <- ivreg(logpgp95 ~ avexpr| logem4, data = data4.r)
reg4.A2 <- ivreg(logpgp95 ~ avexpr + lat_abst| lat_abst + logem4, data = data4.r)

• O coeficiente da variável medindo instituições continua positivo e significativo, de fato, o estimador IV é maior do que o OLS.
• O coeficiente de latitude não é significativo.

• “This result suggests that many previous studies may have found latitude to be a significant determinant of economic performance because it is correlated with institutions (or with the exogenous component of institutions caused by early Colonial experience).”
• “Overall, the results in Table 4 show a large effect of institutions on economic performance. In the rest of the paper, we investigate the robustness of these results.”

• Dados em painel é o termo usado para base de dados onde diversas variáveis referentes a unidades específicas são observadas ao longo do tempo, estas unidade podem ser países, empresas, pessoas, estados, cidades ou qualquer unidade de interesse para a análise em questão.

• Modelo \[y_{it} = \beta_0 + \beta_1 X_{1it} + \beta_2 X_{2it} + \cdots + \beta_k X_{kit} + \varepsilon_{it}\]
• Onde \(\varepsilon_{it} = a_i + v_{it}\)

• Estimador de dados empilhados (pooling): Equivale a estimar por mínimos quadrados ignorando a estrutura de dados em painel. Pode ser usado quando é razoável supor que o intercepto é o mesmo para todas as unidades.
• Estimador de efeitos aleatórios: Deve ser usado quando o termo \(a_i\) é aleatório, mas não é correlacionado com os erros.

• Estimador de efeitos fixos: Deve ser usado quando não é razoável supor que o termo \(a_i\) não é correlacionado com o erro, especificamente quando é válido supor que cada unidade tem seu próprio intercepto.
• Estimador de primeiras diferenças: Também pode ser aplicado quando é válida a hipótese de interceptos diferentes para cada unidade.

• Sob hipóteses razoáveis e com uma amostra grande o estimador de efeitos aleatórios e preferível ao pooling e o estimador de efeitos fixos é preferível ao de primeiras diferenças.
• O teste de Hausman nos ajuda a escolher entre os estimadores de efeito fixos e de efeitos aleatórios.

• Taxar bebidas alcoolicas e mortes em acidentes de trânsito nas estradas.
• Stock e Watson (2007)

data("Fatalities", package = "AER")

Fatalities <- Fatalities %>%
  mutate(frate = fatal /pop * 10000)

• Para começar vamos analisar o primeiro e o último ano da amostra.

mod82 <- lm(frate ~ beertax, data = Fatalities, subset = year == 1982)

## 
## Call:
## lm(formula = frate ~ beertax, data = Fatalities, subset = year == 
##     1982)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9356 -0.4480 -0.1068  0.2295  2.1716 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.0104     0.1391  14.455   <2e-16 ***
## beertax       0.1485     0.1884   0.788    0.435    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6705 on 46 degrees of freedom
## Multiple R-squared:  0.01332,    Adjusted R-squared:  -0.008126 
## F-statistic: 0.6212 on 1 and 46 DF,  p-value: 0.4347

mod88 <- lm(frate ~ beertax, data = Fatalities, subset = year == 1988)

## 
## Call:
## lm(formula = frate ~ beertax, data = Fatalities, subset = year == 
##     1988)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72931 -0.36028 -0.07132  0.39938  1.35783 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.8591     0.1060  17.540   <2e-16 ***
## beertax       0.4388     0.1645   2.668   0.0105 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4903 on 46 degrees of freedom
## Multiple R-squared:  0.134,  Adjusted R-squared:  0.1152 
## F-statistic: 7.118 on 1 and 46 DF,  p-value: 0.0105

• Vamos fazer regressões para todos os anos!

library(purrr)
library(broom)

reg_purrr <- Fatalities %>%
  select(year, frate, beertax) %>%
  nest(data = -year) %>%
  mutate(fit = map(data, ~ lm(frate ~ beertax, data = .x)),
                            tid = map(fit, tidy)) %>%
  unnest(tid) %>%
  select(year, coef = term, est = estimate, p.valor = p.value) %>%
  filter(coef == "beertax")

## # A tibble: 7 × 4
##   year  coef      est p.valor
##   <fct> <chr>   <dbl>   <dbl>
## 1 1982  beertax 0.148 0.435  
## 2 1983  beertax 0.299 0.0822 
## 3 1984  beertax 0.400 0.0114 
## 4 1985  beertax 0.392 0.0148 
## 5 1986  beertax 0.480 0.00464
## 6 1987  beertax 0.483 0.00658
## 7 1988  beertax 0.439 0.0105

• Em todos os anos os valores estimados para o coeficiente de beertax foram positivos, em alguns anos foram significativos.
• E se fizessemos uma regressão com os dados de todos os anos (pooling)?

reg.pool.lm <- lm(frate ~ beertax, data = Fatalities)

## 
## Call:
## lm(formula = frate ~ beertax, data = Fatalities)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09060 -0.37768 -0.09436  0.28548  2.27643 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.85331    0.04357  42.539  < 2e-16 ***
## beertax      0.36461    0.06217   5.865 1.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5437 on 334 degrees of freedom
## Multiple R-squared:  0.09336,    Adjusted R-squared:  0.09065 
## F-statistic: 34.39 on 1 and 334 DF,  p-value: 1.082e-08

• A taxa sobre a cerveja é significativa para explicar mortes no trânsito, mas o sinal é positivo!

• O problema pode ser a presença de heterogeneidade não observada, especificamente em vez de estimar o modelo: \[\mbox{frate}_{it} = \alpha + \beta \ \mbox{beertax}_{it} + \varepsilon_{it}\]

• vamos estimar o modelo: \[\mbox{frate}_{it} = \alpha_{i} + \beta \ \mbox{beertax}_{it} + \varepsilon_{it}\] ## Dados em painel

• Para estimar o novo modelo usaremos o pacote plm. Esse pacote facilita o trabalho com dados em painel e permite o uso de diversos estimadores modificando o argumento da função plm.

• Por exemplo, para estimar o modelo pooling como fizemos anteriormente basta usar a função plm como o argumento definido como “pooling”.

library(plm)

## 
## Attaching package: 'plm'

## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead

reg.pool <- plm(frate ~ beertax, Fatalities, model = "pooling")

## Pooling Model
## 
## Call:
## plm(formula = frate ~ beertax, data = Fatalities, model = "pooling")
## 
## Balanced Panel: n = 48, T = 7, N = 336
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.09060 -0.37768 -0.09436  0.28548  2.27643 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 1.853308   0.043567 42.5391 < 2.2e-16 ***
## beertax     0.364605   0.062170  5.8647 1.082e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    108.92
## Residual Sum of Squares: 98.747
## R-Squared:      0.093363
## Adj. R-Squared: 0.090648
## F-statistic: 34.3943 on 1 and 334 DF, p-value: 1.0822e-08

• Uma maneira relativamente simples de retirar o parâmetro \(\alpha_i\) do modelo é trabalhar com diferenças, no caso: \[\mbox{frate}_{it} - \mbox{frate}_{it-\tau} = \alpha_{i} + \beta \ \mbox{beertax}_{it} + \varepsilon_{it} - \alpha_{i} - \beta \ \mbox{beertax}_{it-\tau} - \varepsilon_{it-\tau}\]
• \[\Rightarrow \Delta_\tau \mbox{frate}_i = \beta \Delta_\tau \mbox{beertax}_i + \Delta_\tau \varepsilon_i\]

reg.dif <- plm(diff(frate,1) ~ diff(beertax,1), Fatalities, model = "pooling")

## Pooling Model
## 
## Call:
## plm(formula = diff(frate, 1) ~ diff(beertax, 1), data = Fatalities, 
##     model = "pooling")
## 
## Balanced Panel: n = 48, T = 6, N = 288
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.697391 -0.106403  0.010073  0.109465  0.606420 
## 
## Coefficients:
##                    Estimate Std. Error t-value Pr(>|t|)
## (Intercept)      -0.0031368  0.0119115 -0.2633   0.7925
## diff(beertax, 1)  0.0136878  0.2852511  0.0480   0.9618
## 
## Total Sum of Squares:    11.213
## Residual Sum of Squares: 11.213
## R-Squared:      8.0509e-06
## Adj. R-Squared: -0.0034884
## F-statistic: 0.00230256 on 1 and 286 DF, p-value: 0.96176

reg.dif <- plm(frate ~ beertax, Fatalities, model = "fd")

## Oneway (individual) effect First-Difference Model
## 
## Call:
## plm(formula = frate ~ beertax, data = Fatalities, model = "fd")
## 
## Balanced Panel: n = 48, T = 7, N = 336
## Observations used in estimation: 288
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.697391 -0.106403  0.010073  0.109465  0.606420 
## 
## Coefficients:
##               Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -0.0031368  0.0119115 -0.2633   0.7925
## beertax      0.0136878  0.2852511  0.0480   0.9618
## 
## Total Sum of Squares:    11.213
## Residual Sum of Squares: 11.213
## R-Squared:      8.0509e-06
## Adj. R-Squared: -0.0034884
## F-statistic: 0.00230256 on 1 and 286 DF, p-value: 0.96176

• Outra abordagem para o problema é estimar os interceptos de cada estado.

reg.lsdv <- lm(frate ~ beertax + state -1, Fatalities)

## 
## Call:
## lm(formula = frate ~ beertax + state - 1, data = Fatalities)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58696 -0.08284 -0.00127  0.07955  0.89780 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## beertax -0.65587    0.18785  -3.491 0.000556 ***
## stateal  3.47763    0.31336  11.098  < 2e-16 ***
## stateaz  2.90990    0.09254  31.445  < 2e-16 ***
## statear  2.82268    0.13213  21.364  < 2e-16 ***
## stateca  1.96816    0.07401  26.594  < 2e-16 ***
## stateco  1.99335    0.08037  24.802  < 2e-16 ***
## statect  1.61537    0.08391  19.251  < 2e-16 ***
## statede  2.17003    0.07746  28.016  < 2e-16 ***
## statefl  3.20950    0.22151  14.489  < 2e-16 ***
## statega  4.00223    0.46403   8.625 4.43e-16 ***
## stateid  2.80861    0.09877  28.437  < 2e-16 ***
## stateil  1.51601    0.07848  19.318  < 2e-16 ***
## statein  2.01609    0.08867  22.736  < 2e-16 ***
## stateia  1.93370    0.10222  18.918  < 2e-16 ***
## stateks  2.25441    0.10863  20.753  < 2e-16 ***
## stateky  2.26011    0.08046  28.089  < 2e-16 ***
## statela  2.63051    0.16266  16.171  < 2e-16 ***
## stateme  2.36968    0.16006  14.805  < 2e-16 ***
## statemd  1.77119    0.08246  21.480  < 2e-16 ***
## statema  1.36788    0.08648  15.818  < 2e-16 ***
## statemi  1.99310    0.11663  17.089  < 2e-16 ***
## statemn  1.58042    0.09363  16.880  < 2e-16 ***
## statems  3.44855    0.20936  16.472  < 2e-16 ***
## statemo  2.18137    0.09252  23.576  < 2e-16 ***
## statemt  3.11724    0.09441  33.017  < 2e-16 ***
## statene  1.95545    0.10551  18.534  < 2e-16 ***
## statenv  2.87686    0.08106  35.492  < 2e-16 ***
## statenh  2.22318    0.14114  15.751  < 2e-16 ***
## statenj  1.37188    0.07333  18.709  < 2e-16 ***
## statenm  3.90401    0.10154  38.449  < 2e-16 ***
## stateny  1.29096    0.07563  17.070  < 2e-16 ***
## statenc  3.18717    0.25173  12.661  < 2e-16 ***
## statend  1.85419    0.10193  18.191  < 2e-16 ***
## stateoh  1.80321    0.10193  17.691  < 2e-16 ***
## stateok  2.93257    0.18428  15.913  < 2e-16 ***
## stateor  2.30963    0.08117  28.453  < 2e-16 ***
## statepa  1.71016    0.08648  19.776  < 2e-16 ***
## stateri  1.21258    0.07753  15.640  < 2e-16 ***
## statesc  4.03480    0.35479  11.372  < 2e-16 ***
## statesd  2.47391    0.14121  17.519  < 2e-16 ***
## statetn  2.60197    0.09162  28.398  < 2e-16 ***
## statetx  2.56016    0.10853  23.589  < 2e-16 ***
## stateut  2.31368    0.15453  14.972  < 2e-16 ***
## statevt  2.51159    0.13973  17.975  < 2e-16 ***
## stateva  2.18745    0.14664  14.917  < 2e-16 ***
## statewa  1.81811    0.08233  22.084  < 2e-16 ***
## statewv  2.58088    0.10767  23.971  < 2e-16 ***
## statewi  1.71836    0.07746  22.185  < 2e-16 ***
## statewy  3.24913    0.07233  44.922  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1899 on 287 degrees of freedom
## Multiple R-squared:  0.9931, Adjusted R-squared:  0.992 
## F-statistic: 847.8 on 49 and 287 DF,  p-value: < 2.2e-16

• O estimador de efeitos fixos (within) fornece resultados equivalentes a colocar uma constrante para cada estado, porém de forma mais compacta.
• Para estimar por efeitos fixos a função plm pode ser usada com o argumento model definido como “within” ou sem definir o argumento model, pois “within” é o valor padrão do argumento model.

reg.fe <- plm(frate ~ beertax, Fatalities, model = "within")

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = frate ~ beertax, data = Fatalities, model = "within")
## 
## Balanced Panel: n = 48, T = 7, N = 336
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.5869619 -0.0828376 -0.0012701  0.0795454  0.8977960 
## 
## Coefficients:
##         Estimate Std. Error t-value Pr(>|t|)    
## beertax -0.65587    0.18785 -3.4915 0.000556 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    10.785
## Residual Sum of Squares: 10.345
## R-Squared:      0.040745
## Adj. R-Squared: -0.11969
## F-statistic: 12.1904 on 1 and 287 DF, p-value: 0.00055597

reg.fe <- plm(frate ~ beertax, Fatalities)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = frate ~ beertax, data = Fatalities)
## 
## Balanced Panel: n = 48, T = 7, N = 336
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.5869619 -0.0828376 -0.0012701  0.0795454  0.8977960 
## 
## Coefficients:
##         Estimate Std. Error t-value Pr(>|t|)    
## beertax -0.65587    0.18785 -3.4915 0.000556 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    10.785
## Residual Sum of Squares: 10.345
## R-Squared:      0.040745
## Adj. R-Squared: -0.11969
## F-statistic: 12.1904 on 1 and 287 DF, p-value: 0.00055597

• A função fixef retorna os valores estimados para os efeitos fixos.
• Se o argumento type for definido como “dfirst” a função retorna a diferença para a primeiro unidade (estado)
• Se o argumento type for definido como “dmean” a função retorna o desvio da média

fixef(reg.fe)

##     al     az     ar     ca     co     ct     de     fl     ga     id     il 
## 3.4776 2.9099 2.8227 1.9682 1.9933 1.6154 2.1700 3.2095 4.0022 2.8086 1.5160 
##     in     ia     ks     ky     la     me     md     ma     mi     mn     ms 
## 2.0161 1.9337 2.2544 2.2601 2.6305 2.3697 1.7712 1.3679 1.9931 1.5804 3.4486 
##     mo     mt     ne     nv     nh     nj     nm     ny     nc     nd     oh 
## 2.1814 3.1172 1.9555 2.8769 2.2232 1.3719 3.9040 1.2910 3.1872 1.8542 1.8032 
##     ok     or     pa     ri     sc     sd     tn     tx     ut     vt     va 
## 2.9326 2.3096 1.7102 1.2126 4.0348 2.4739 2.6020 2.5602 2.3137 2.5116 2.1874 
##     wa     wv     wi     wy 
## 1.8181 2.5809 1.7184 3.2491

fixef(reg.fe, type = "dfirst")

##       az       ar       ca       co       ct       de       fl       ga 
## -0.56773 -0.65495 -1.50947 -1.48428 -1.86226 -1.30760 -0.26813  0.52460 
##       id       il       in       ia       ks       ky       la       me 
## -0.66902 -1.96162 -1.46154 -1.54393 -1.22322 -1.21752 -0.84712 -1.10795 
##       md       ma       mi       mn       ms       mo       mt       ne 
## -1.70644 -2.10975 -1.48453 -1.89721 -0.02908 -1.29626 -0.36039 -1.52218 
##       nv       nh       nj       nm       ny       nc       nd       oh 
## -0.60077 -1.25445 -2.10575  0.42637 -2.18667 -0.29047 -1.62344 -1.67442 
##       ok       or       pa       ri       sc       sd       tn       tx 
## -0.54506 -1.16800 -1.76747 -2.26505  0.55717 -1.00372 -0.87566 -0.91747 
##       ut       vt       va       wa       wv       wi       wy 
## -1.16395 -0.96604 -1.29018 -1.65952 -0.89675 -1.75927 -0.22850

fixef(reg.fe, type = "dmean")

##         al         az         ar         ca         co         ct         de 
##  1.1005552  0.5328284  0.4456036 -0.4089136 -0.3837252 -0.7617021 -0.2070468 
##         fl         ga         id         il         in         ia         ks 
##  0.8324250  1.6251582  0.4315327 -0.8610674 -0.3609864 -0.4433770 -0.1226608 
##         ky         la         me         md         ma         mi         mn 
## -0.1169618  0.2534390 -0.0073922 -0.6058845 -1.0091910 -0.3839717 -0.7966578 
##         ms         mo         mt         ne         nv         nh         nj 
##  1.0714754 -0.1957065  0.7401637 -0.4216227  0.4997802 -0.1538992 -1.0051943 
##         nm         ny         nc         nd         oh         ok         or 
##  1.5269302 -1.0861148  0.8100902 -0.5228841 -0.5738641  0.5554942 -0.0674447 
##         pa         ri         sc         sd         tn         tx         ut 
## -0.6669110 -1.1644991  1.6577289  0.0968344  0.2248966  0.1830818 -0.0633953 
##         vt         va         wa         wv         wi         wy 
##  0.1345114 -0.1896280 -0.5589686  0.2038012 -0.6587112  0.8720515

• Via de regra, é recomendável usar a função pdata.frame para definir o data.frame como um painel.
• O argumento index permite definir as dimensões dos indivíduos (estados, países, famílias, firmas, etc) e a dimensão de tempo. Se o argumento index não estiver definido o R entende que a primeira coluna define os indivíduos e a segunda o tempo.

pFatalities <- pdata.frame(Fatalities, index = c("state","year"))

pdim(pFatalities)

## Balanced Panel: n = 48, T = 7, N = 336

• Para usar o estimador de efeitos aleatórios o argumento model deve ser definido como “random”.

reg.re <- plm(frate ~ beertax, Fatalities, model = "random")

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = frate ~ beertax, data = Fatalities, model = "random")
## 
## Balanced Panel: n = 48, T = 7, N = 336
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.03605 0.18986 0.119
## individual    0.26604 0.51579 0.881
## theta: 0.8622
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.471090 -0.120045 -0.021530  0.091011  0.964350 
## 
## Coefficients:
##              Estimate Std. Error z-value Pr(>|z|)    
## (Intercept)  2.067141   0.099971 20.6773   <2e-16 ***
## beertax     -0.052016   0.124176 -0.4189   0.6753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    12.648
## Residual Sum of Squares: 12.642
## R-Squared:      0.00052508
## Adj. R-Squared: -0.0024674
## Chisq: 0.175467 on 1 DF, p-value: 0.6753

• Podemos usar a função map_df para estimar o modelo com vários estimadores e comparar os resultados.

models <- c("within", "random", "pooling", "fd")

map_df(models, ~ tidy(plm(frate ~ beertax, Fatalities, model=.x)), .id="modelo") %>%
  filter(term == "beertax") %>%
  mutate(modelo = models)

## # A tibble: 4 × 6
##   modelo  term    estimate std.error statistic      p.value
##   <chr>   <chr>      <dbl>     <dbl>     <dbl>        <dbl>
## 1 within  beertax  -0.656     0.188    -3.49   0.000556    
## 2 random  beertax  -0.0520    0.124    -0.419  0.675       
## 3 pooling beertax   0.365     0.0622    5.86   0.0000000108
## 4 fd      beertax   0.0137    0.285     0.0480 0.962

• O teste de Haussman é usado para escolher entre modelos de efeitos fixos e efeitos aleatórios
• \(H_0: E[a_i|X] = 0\)
• Não rejeitar \(H_0\): Estimador de efeitos aleatórios é consistente e assintóticamente eficiente, estimador de efeitos fixos é apenas consistente. Usar efeitos aleatórios.
• Rejeitar \(H_0\): Estimador de efeitos aleatórios é inconsistente, estimador de efeitos fixos é consistente. Usar efeitos fixos.

phtest(reg.fe, reg.re)

## 
##  Hausman Test
## 
## data:  frate ~ beertax
## chisq = 18.353, df = 1, p-value = 1.835e-05
## alternative hypothesis: one model is inconsistent

• No exemplo é recomendado usar efeitos fixos (p.valor < 0,05).

• Como último exemplo considere os dados de produção de ladrilhos de empresas em duas regiões do Egito.
• O objetivo é usar os dados de capital e trabalho para estimar a função de produção.

data(Tileries, package = "pder")

map_df(models, ~ tidy(plm(log(output) ~ log(labor) + machine, Tileries, model=.x, 
                          subset = area == "fayoum")), 
        .id="modelo") %>% 
  filter(term %in% c("log(labor)", "machine")) %>%
  mutate(modelo = rep(models,each=2))

## # A tibble: 8 × 6
##   modelo  term       estimate std.error statistic   p.value
##   <chr>   <chr>         <dbl>     <dbl>     <dbl>     <dbl>
## 1 within  log(labor) 0.917       0.0466  19.7     2.93e- 45
## 2 within  machine    0.000107    0.0124   0.00864 9.93e-  1
## 3 random  log(labor) 0.950       0.0405  23.4     1.69e-121
## 4 random  machine    0.00184     0.0107   0.172   8.63e-  1
## 5 pooling log(labor) 0.965       0.0382  25.3     3.99e- 60
## 6 pooling machine    0.00224     0.0100   0.224   8.23e-  1
## 7 fd      log(labor) 0.811       0.0402  20.2     2.07e- 46
## 8 fd      machine    0.00571     0.0124   0.462   6.45e-  1

• Repare que os valores estimados foram semelhantes para os diversos modelos usados, isso sugere que heterogenidade não observada não é um fator relevante.

reg.tl.fe <- plm(log(output) ~ log(labor) + machine, Tileries, model="within", 
                 subset = area == "fayoum")
reg.tl.re <- plm(log(output) ~ log(labor) + machine, Tileries, model="random", 
                 subset = area == "fayoum")

phtest(reg.tl.fe, reg.tl.re)

## 
##  Hausman Test
## 
## data:  log(output) ~ log(labor) + machine
## chisq = 2.4218, df = 2, p-value = 0.2979
## alternative hypothesis: one model is inconsistent

• O teste de Hausman não rejeita a hipótese nula, portanto é recomendável usar efeitos aleatórios.