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.

• O método de controle sintético, apresentado por Abadie e Gardeazabal (2003) e Abadie, Diamond e Hainmueller (2010 e 2011), tenta resolver o problema do contrafactual comparando a tendência na região atingida pelo choque ou pela política com a tendência em uma região sintética composta a partir de diversas regiões observadas.
• Na definição em Abadie e Gardeazabal (2003) e Abadie , Diamond e Hainmueller (2010), a unidade de controle sintético é uma média ponderada das unidades de controle disponíveis que melhor aproxima as características, inclusive de tendência, da variável tratada antes do tratamento.

• Efeitos do terrorismo no país Basco
• Abadie, Alberto, Alexis Diamond and Jens Hainmueller. “Synth: An R Package for Synthetic Control Methods in Comparative Case Studies.” Journal of Statistical Software, June 2011, Volume 42, Issue 13, p.1-17.

• Unidade tratada: País Basco
• Unidades de controle: outras regiões da Espanha

library(Synth)

## Warning: package 'Synth' was built under R version 4.3.3

## ##
## ## Synth Package: Implements Synthetic Control Methods.

## ## See https://web.stanford.edu/~jhain/synthpage.html for additional information.

data("basque")

• Variável de interesse: PIB per capita

• Preditores:
• 1964 a 1969 averages for gross total investment/GDP (invest).
• 1964 a 1969 averages for the share of the working-age population that was illiterate (school.illit), the share with up to primary school education (school.prim), the with some high school (school.med), the share with high school (school.high), and the share with more than high school (school.post.high).
• 1961 a 1969 averages for six industrial-sector shares as a percentage of total production (these variables are named with a sec. prefix).
• 1960 a 1969 averages for real GDP per-capita (gdpcap) measured in thousands of 1986 USD.
• 1969 population density measured in persons per square kilometer (popdens).

• A função dataprep prepara os dados.
• X1: preditores da unidade tratada
• X0: preditores das unidades de controle
• Z0: variável de interesse para unidade tratada
• Z1: variável de interesse para as unidades de controle

dataprep.out <- dataprep(
  foo = basque,
  predictors = c("school.illit", "school.prim", "school.med",
                 "school.high", "school.post.high", "invest"),
  predictors.op = "mean",
  time.predictors.prior = 1964:1969,
  special.predictors = list(
    list("gdpcap", 1960:1969 , "mean"),
    list("sec.agriculture", seq(1961, 1969, 2), "mean"),
    list("sec.energy", seq(1961, 1969, 2), "mean"),
    list("sec.industry", seq(1961, 1969, 2), "mean"),
    list("sec.construction", seq(1961, 1969, 2), "mean"),
    list("sec.services.venta", seq(1961, 1969, 2), "mean"),
    list("sec.services.nonventa", seq(1961, 1969, 2), "mean"),
    list("popdens", 1969, "mean")),
  dependent = "gdpcap",
  unit.variable = "regionno",
  unit.names.variable = "regionname",
  time.variable = "year",
  treatment.identifier = 17,
  controls.identifier = c(2:16, 18),
  time.optimize.ssr = 1960:1969,
  time.plot = 1955:1997
)

class(dataprep.out)

## [1] "list"

names(dataprep.out)

## [1] "X0"                "X1"                "Z0"               
## [4] "Z1"                "Y0plot"            "Y1plot"           
## [7] "names.and.numbers" "tag"

dataprep.out$X0

##                                                   2          3          4
## school.illit                             863.389160  73.121226  31.488423
## school.prim                             3062.424886 728.578929 670.909393
## school.med                               155.565318  44.215389  46.398482
## school.high                               57.266496  16.091676  14.799358
## school.post.high                          27.278924   8.684416   6.424505
## invest                                    19.320031  21.577486  22.769643
## special.gdpcap.1960.1969                   2.560747   3.699907   3.733876
## special.sec.agriculture.1961.1969         24.194000  21.726000  12.362000
## special.sec.energy.1961.1969               2.774000   6.278000  18.648000
## special.sec.industry.1961.1969            18.276000  22.780000  24.126000
## special.sec.construction.1961.1969         8.130000   7.832000   9.006000
## special.sec.services.venta.1961.1969      38.186000  34.289999  30.152000
## special.sec.services.nonventa.1961.1969    8.444000   7.096000   5.708000
## special.popdens.1969                      68.510002  24.040001  98.739998
##                                                  5          6          7
## school.illit                             47.903906 128.308287   9.394911
## school.prim                             300.813619 522.094955 289.571732
## school.med                               20.045204  45.447489  24.106414
## school.high                               5.921604  12.086849   7.304420
## school.post.high                          3.680154   5.844122   2.885214
## invest                                   24.441712  25.954247  29.071211
## special.gdpcap.1960.1969                  5.215974   3.051014   3.871173
## special.sec.agriculture.1961.1969        13.130000  19.944000  15.922000
## special.sec.energy.1961.1969              2.076000   7.818000   2.894000
## special.sec.industry.1961.1969           18.258000   9.816000  36.530000
## special.sec.construction.1961.1969        8.294000   8.670000   5.976000
## special.sec.services.venta.1961.1969     51.752000  45.278001  33.484000
## special.sec.services.nonventa.1961.1969   6.494000   8.482000   5.198000
## special.popdens.1969                    104.169998 148.250000  87.389999
##                                                   8           9          10
## school.illit                             105.508144  254.449707  277.935237
## school.prim                             1766.928691 1035.269368 2883.361735
## school.med                               102.299929   33.369543  215.165476
## school.high                               37.787317   16.539727   63.608315
## school.post.high                          19.173427    6.308165   32.374611
## invest                                    19.652509   17.970690   22.520380
## special.gdpcap.1960.1969                   2.807062    2.240688    5.147008
## special.sec.agriculture.1961.1969         29.718000   36.086001    6.936000
## special.sec.energy.1961.1969               7.818000    5.180000    2.880000
## special.sec.industry.1961.1969            16.474000   17.178000   40.056000
## special.sec.construction.1961.1969         7.144000    5.810000    6.706000
## special.sec.services.venta.1961.1969      30.844000   29.238000   39.068000
## special.sec.services.nonventa.1961.1969    8.000000    6.506000    4.356000
## special.popdens.1969                      28.770000   22.379999  153.119995
##                                                  11         12          13
## school.illit                             266.967601 177.727966  234.752001
## school.prim                             1695.117126 692.888550 1666.930420
## school.med                               103.959094  24.048450   77.464478
## school.high                               34.275772  11.271676   29.341941
## school.post.high                          16.822106   4.570566   13.543060
## invest                                    23.702696  20.239484   20.606269
## special.gdpcap.1960.1969                   3.843245   1.926478    2.516288
## special.sec.agriculture.1961.1969         18.582000  37.948000   28.862000
## special.sec.energy.1961.1969               2.698000   2.646000    5.328000
## special.sec.industry.1961.1969            28.046000  10.886000   17.882000
## special.sec.construction.1961.1969         6.170000   8.770000    7.032000
## special.sec.services.venta.1961.1969      39.064001  31.720000   33.714001
## special.sec.services.nonventa.1961.1969    5.444000   8.038000    7.188000
## special.popdens.1969                     128.699997  28.969999   91.019997
##                                                  14         15         16
## school.illit                             133.158962 105.877481  13.438050
## school.prim                             1856.074280 431.746806 280.002396
## school.med                               269.961647  27.907434  20.450936
## school.high                               62.458658   9.082540   6.550353
## school.post.high                          57.704985   4.428528   4.190368
## invest                                    16.234543  20.706271  19.955159
## special.gdpcap.1960.1969                   5.976692   2.899557   3.996087
## special.sec.agriculture.1961.1969          1.864000  19.352000  24.704000
## special.sec.energy.1961.1969               2.074000  10.228000   2.740000
## special.sec.industry.1961.1969            23.834000  19.644000  28.398000
## special.sec.construction.1961.1969         8.358000   6.290000   6.982000
## special.sec.services.venta.1961.1969      52.714000  36.177999  30.468000
## special.sec.services.nonventa.1961.1969   11.162000   8.312000   6.710000
## special.popdens.1969                     442.450012  73.360001  44.169998
##                                                 18
## school.illit                              9.151832
## school.prim                             152.269465
## school.med                                9.760035
## school.high                               3.377170
## school.post.high                          1.732763
## invest                                   18.055932
## special.gdpcap.1960.1969                  3.809205
## special.sec.agriculture.1961.1969        30.322001
## special.sec.energy.1961.1969              2.888000
## special.sec.industry.1961.1969           26.612000
## special.sec.construction.1961.1969        5.238000
## special.sec.services.venta.1961.1969     28.300000
## special.sec.services.nonventa.1961.1969   6.644000
## special.popdens.1969                     46.580002

dataprep.out$X1

##                                                  17
## school.illit                              39.888465
## school.prim                             1031.742299
## school.med                                90.358668
## school.high                               25.727525
## school.post.high                          13.479720
## invest                                    24.647383
## special.gdpcap.1960.1969                   5.285468
## special.sec.agriculture.1961.1969          6.844000
## special.sec.energy.1961.1969               4.106000
## special.sec.industry.1961.1969            45.082000
## special.sec.construction.1961.1969         6.150000
## special.sec.services.venta.1961.1969      33.754000
## special.sec.services.nonventa.1961.1969    4.072000
## special.popdens.1969                     246.889999

dataprep.out$Z0

##             2        3        4        5        6        7        8        9
## 1960 2.010140 2.881462 2.967295 4.058841 2.357684 3.137032 2.138817 1.667524
## 1961 2.129177 3.099543 3.143887 4.360254 2.445730 3.327621 2.239503 1.752428
## 1962 2.280348 3.359183 3.373536 4.646173 2.648243 3.555341 2.454227 1.920451
## 1963 2.431020 3.614182 3.597258 4.911525 2.844759 3.771423 2.672237 2.091902
## 1964 2.508855 3.680091 3.672594 5.050700 2.951157 3.839403 2.777778 2.182591
## 1965 2.584690 3.745287 3.743359 5.184662 3.054199 3.906098 2.882176 2.274707
## 1966 2.694444 3.883319 3.909383 5.466795 3.231791 4.032133 2.988075 2.378392
## 1967 2.802342 4.016138 4.073122 5.737646 3.403385 4.155955 3.094544 2.482362
## 1968 2.987361 4.243645 4.308626 6.161454 3.660312 4.375893 3.302271 2.709083
## 1969 3.179092 4.476221 4.549700 6.581691 3.912882 4.610826 3.520994 2.947444
##            10       11       12       13       14       15       16       18
## 1960 4.241788 3.219294 1.535847 1.983290 5.161097 2.118609 3.163525 2.969866
## 1961 4.575335 3.362468 1.596258 2.005784 5.632605 2.305484 3.335904 3.153171
## 1962 4.838046 3.569980 1.705584 2.185661 5.840831 2.521422 3.623393 3.404384
## 1963 5.081334 3.765210 1.817695 2.366395 6.024493 2.739074 3.894816 3.669238
## 1964 5.158098 3.823693 1.882819 2.458797 6.099329 2.851257 3.985147 3.803985
## 1965 5.223651 3.874179 1.948872 2.549700 6.152028 2.965938 4.072979 3.921808
## 1966 5.332477 3.978149 2.032633 2.669666 6.110469 3.099186 4.210011 4.032705
## 1967 5.429449 4.073408 2.117609 2.787846 6.057341 3.227292 4.352399 4.160311
## 1968 5.674379 4.279777 2.245501 2.978363 6.253142 3.461154 4.556984 4.373036
## 1969 5.915524 4.486290 2.381962 3.177378 6.435590 3.706155 4.765710 4.603542

dataprep.out$Z1

##            17
## 1960 4.285918
## 1961 4.574336
## 1962 4.898957
## 1963 5.197015
## 1964 5.338903
## 1965 5.465153
## 1966 5.545916
## 1967 5.614896
## 1968 5.852185
## 1969 6.081405

dataprep.out$Y0plot

##             2        3        4         5        6        7        8        9
## 1955 1.688732 2.288775 2.502928  3.143959 1.914382 2.559412 1.729149 1.327764
## 1956 1.758498 2.445159 2.615538  3.347758 2.071837 2.693873 1.838332 1.415096
## 1957 1.827621 2.603399 2.725793  3.549629 2.226078 2.820337 1.947658 1.503570
## 1958 1.852756 2.639032 2.751857  3.642673 2.220866 2.879035 1.971365 1.531420
## 1959 1.878035 2.677092 2.777421  3.734862 2.213439 2.943730 1.995144 1.559340
## 1960 2.010140 2.881462 2.967295  4.058841 2.357684 3.137032 2.138817 1.667524
## 1961 2.129177 3.099543 3.143887  4.360254 2.445730 3.327621 2.239503 1.752428
## 1962 2.280348 3.359183 3.373536  4.646173 2.648243 3.555341 2.454227 1.920451
## 1963 2.431020 3.614182 3.597258  4.911525 2.844759 3.771423 2.672237 2.091902
## 1964 2.508855 3.680091 3.672594  5.050700 2.951157 3.839403 2.777778 2.182591
## 1965 2.584690 3.745287 3.743359  5.184662 3.054199 3.906098 2.882176 2.274707
## 1966 2.694444 3.883319 3.909383  5.466795 3.231791 4.032133 2.988075 2.378392
## 1967 2.802342 4.016138 4.073122  5.737646 3.403385 4.155955 3.094544 2.482362
## 1968 2.987361 4.243645 4.308626  6.161454 3.660312 4.375893 3.302271 2.709083
## 1969 3.179092 4.476221 4.549700  6.581691 3.912882 4.610826 3.520994 2.947444
## 1970 3.354327 4.596258 4.631605  6.887032 4.224936 4.791417 3.670523 3.136890
## 1971 3.522922 4.723722 4.699657  7.168666 4.493573 4.969652 3.821837 3.319623
## 1972 3.756213 5.001285 5.022565  7.570694 4.781919 5.155384 4.081762 3.629177
## 1973 3.987718 5.283491 5.345615  7.956298 5.054913 5.338760 4.345402 3.946087
## 1974 4.051842 5.442374 5.502499  7.972722 4.967795 5.508355 4.464724 4.028135
## 1975 4.112182 5.600971 5.651957  7.974150 4.881534 5.675236 4.581120 4.111968
## 1976 4.221794 5.736575 5.592259  8.034847 4.948443 5.796058 4.732719 4.260711
## 1977 4.331691 5.866396 5.538703  8.080548 5.007926 5.902956 4.884676 4.412168
## 1978 4.305199 5.925021 5.614610  8.041988 5.040917 5.916595 4.946801 4.446301
## 1979 4.226935 5.952799 5.704299  8.004856 5.233505 5.906455 4.920380 4.408741
## 1980 4.211511 5.974722 5.796987  8.259783 5.437089 5.899100 4.893745 4.328763
## 1981 4.207012 6.011140 5.935661  8.531134 5.658241 5.914882 4.881034 4.261282
## 1982 4.249072 6.133605 5.849757  8.722508 5.687304 5.916881 4.923093 4.343187
## 1983 4.291631 6.260854 5.769066  8.925307 5.719866 5.941660 4.970722 4.424664
## 1984 4.358683 6.372894 5.887318  9.275921 5.801342 6.028849 5.113468 4.550057
## 1985 4.426593 6.495501 6.011283  9.652242 5.885604 6.138318 5.261997 4.677664
## 1986 4.663239 6.926521 6.234790 10.257783 6.256784 6.420451 5.647315 4.980648
## 1987 4.900671 7.358612 6.465296 10.823336 6.612682 6.713225 6.044630 5.295559
## 1988 5.159597 7.802770 6.688160 11.120395 6.977007 7.023422 6.354399 5.677878
## 1989 5.417738 8.242645 6.913525 11.408169 7.337903 7.333619 6.674022 6.065339
## 1990 5.585261 8.458298 6.983148 11.512425 7.345044 7.450729 6.870323 6.279420
## 1991 5.749214 8.668238 7.040631 11.679520 7.347187 7.596401 7.063196 6.474507
## 1992 5.641245 8.466866 6.922808 11.319623 7.220080 7.462154 7.045487 6.330691
## 1993 5.534918 8.256927 6.798343 10.969723 7.092188 7.327906 7.027635 6.188589
## 1994 5.638817 8.573979 6.954156 11.419594 7.410740 7.550700 7.074264 6.230934
## 1995 5.720723 8.846758 7.116467 11.773779 7.616395 7.777064 7.282919 6.328763
## 1996 5.995930 9.096687 7.217224 11.926592 7.817052 7.907741 7.611397 6.614396
## 1997 6.300986 9.518709 7.475721 12.350043 8.060554 8.226935 7.888460 6.865396
##             10       11       12       13        14       15        16
## 1955  3.546630 2.575978 1.243430 1.634676  4.594473 1.679520  2.555127
## 1956  3.690446 2.738503 1.332548 1.725578  4.786632 1.764282  2.698158
## 1957  3.826835 2.899886 1.422451 1.816481  4.963439 1.850328  2.839831
## 1958  3.875678 2.963510 1.440231 1.840903  4.906170 1.887389  2.881891
## 1959  3.921737 3.026207 1.458083 1.865396  4.846401 1.924093  2.930877
## 1960  4.241788 3.219294 1.535847 1.983290  5.161097 2.118609  3.163525
## 1961  4.575335 3.362468 1.596258 2.005784  5.632605 2.305484  3.335904
## 1962  4.838046 3.569980 1.705584 2.185661  5.840831 2.521422  3.623393
## 1963  5.081334 3.765210 1.817695 2.366395  6.024493 2.739074  3.894816
## 1964  5.158098 3.823693 1.882819 2.458797  6.099329 2.851257  3.985147
## 1965  5.223651 3.874179 1.948872 2.549700  6.152028 2.965938  4.072979
## 1966  5.332477 3.978149 2.032633 2.669666  6.110469 3.099186  4.210011
## 1967  5.429449 4.073408 2.117609 2.787846  6.057341 3.227292  4.352399
## 1968  5.674379 4.279777 2.245501 2.978363  6.253142 3.461154  4.556984
## 1969  5.915524 4.486290 2.381962 3.177378  6.435590 3.706155  4.765710
## 1970  6.066838 4.654741 2.518495 3.361968  6.543345 3.906384  4.979434
## 1971  6.227649 4.817124 2.654028 3.534776  6.674736 4.086475  5.199657
## 1972  6.539060 5.138889 2.846401 3.783205  7.086904 4.379463  5.466795
## 1973  6.837975 5.449372 3.045630 4.029991  7.475007 4.666952  5.730506
## 1974  6.987361 5.557984 3.103042 4.132391  7.655670 4.697801  5.984790
## 1975  7.124893 5.655955 3.156741 4.231005  7.816338 4.722722  6.222936
## 1976  7.135390 5.761354 3.262568 4.359040  7.707084 4.784347  6.297772
## 1977  7.142959 5.861040 3.365396 4.485647  7.599972 4.844544  6.384462
## 1978  7.019352 5.810911 3.487075 4.521708  7.390032 4.864610  6.324193
## 1979  7.010997 5.775493 3.515853 4.557483  7.321480 4.836404  6.239646
## 1980  7.078835 5.790988 3.539560 4.581691  7.417166 4.763210  6.208440
## 1981  7.182234 5.899672 3.572765 4.618752  7.532134 4.708083  6.177878
## 1982  7.287204 6.018138 3.610968 4.710797  7.542845 4.811125  6.362039
## 1983  7.397886 6.139817 3.648957 4.806341  7.558555 4.906884  6.544702
## 1984  7.484290 6.236861 3.737146 4.899314  7.699228 5.031991  6.797201
## 1985  7.569980 6.336118 3.828478 4.997786  7.839189 5.154599  7.047772
## 1986  8.077692 6.739360 4.161739 5.277921  8.347615 5.471508  7.449300
## 1987  8.583976 7.144387 4.498572 5.565767  8.849615 5.788560  7.879178
## 1988  9.057412 7.560697 4.769494 5.909526  9.254499 6.124893  8.349758
## 1989  9.525850 7.969152 5.051414 6.254499  9.657955 6.450443  8.803913
## 1990  9.785062 8.138389 5.234076 6.453228  9.806484 6.578192  9.197372
## 1991 10.050700 8.306198 5.398315 6.641603  9.963582 6.710368  9.591545
## 1992  9.837903 8.080548 5.365467 6.544487  9.840046 6.662882  9.345187
## 1993  9.625107 7.857041 5.332619 6.447515  9.718652 6.616324  9.117395
## 1994 10.006427 8.068409 5.439874 6.556484  9.882177 6.784847  9.365895
## 1995 10.339903 8.289061 5.501357 6.688660 10.098543 6.885818  9.758640
## 1996 10.576264 8.429734 5.905813 6.862468 10.322765 7.045416 10.060697
## 1997 11.045416 8.725364 6.224579 7.138532 10.732648 7.295058 10.522708
##             18
## 1955  2.390460
## 1956  2.535204
## 1957  2.680020
## 1958  2.726435
## 1959  2.772851
## 1960  2.969866
## 1961  3.153171
## 1962  3.404384
## 1963  3.669238
## 1964  3.803985
## 1965  3.921808
## 1966  4.032705
## 1967  4.160311
## 1968  4.373036
## 1969  4.603542
## 1970  4.793416
## 1971  4.983362
## 1972  5.230077
## 1973  5.474650
## 1974  5.575836
## 1975  5.676093
## 1976  5.816552
## 1977  5.955798
## 1978  6.066552
## 1979  6.101043
## 1980  6.108683
## 1981  6.139960
## 1982  6.309769
## 1983  6.502142
## 1984  6.626893
## 1985  6.775564
## 1986  7.165096
## 1987  7.580691
## 1988  8.002713
## 1989  8.453299
## 1990  8.858897
## 1991  9.229506
## 1992  9.180948
## 1993  9.132391
## 1994  9.498000
## 1995  9.752213
## 1996 10.056413
## 1997 10.476292

dataprep.out$Y1plot

##             17
## 1955  3.853185
## 1956  3.945658
## 1957  4.033562
## 1958  4.023422
## 1959  4.013782
## 1960  4.285918
## 1961  4.574336
## 1962  4.898957
## 1963  5.197015
## 1964  5.338903
## 1965  5.465153
## 1966  5.545916
## 1967  5.614896
## 1968  5.852185
## 1969  6.081405
## 1970  6.170094
## 1971  6.283633
## 1972  6.555555
## 1973  6.810769
## 1974  7.105184
## 1975  7.377892
## 1976  7.232934
## 1977  7.089831
## 1978  6.786704
## 1979  6.639817
## 1980  6.562839
## 1981  6.500785
## 1982  6.545059
## 1983  6.595330
## 1984  6.761497
## 1985  6.937161
## 1986  7.332191
## 1987  7.742788
## 1988  8.120537
## 1989  8.509711
## 1990  8.776778
## 1991  9.025279
## 1992  8.873893
## 1993  8.718224
## 1994  9.018138
## 1995  9.440874
## 1996  9.686518
## 1997 10.170666

dataprep.out$names.and.numbers

##                      unit.names unit.numbers
## 1   Basque Country (Pais Vasco)           17
## 2                     Andalucia            2
## 3                        Aragon            3
## 4        Principado De Asturias            4
## 5              Baleares (Islas)            5
## 6                      Canarias            6
## 7                     Cantabria            7
## 8               Castilla Y Leon            8
## 9            Castilla-La Mancha            9
## 10                     Cataluna           10
## 11         Comunidad Valenciana           11
## 12                  Extremadura           12
## 13                      Galicia           13
## 14        Madrid (Comunidad De)           14
## 15           Murcia (Region de)           15
## 16 Navarra (Comunidad Foral De)           16
## 17                   Rioja (La)           18

• Pode ser útil em certas ocasiões manipular e modificar X0 e X1 sem precisar voltar ao conjunto de dados original. Para demonstrar, trabalhamos com as cinco variáveis educacionais diferentes (school.illit, school.prim, school.med, school.high, school.post.high) que representam os números, em milhares de indivíduos, com diferentes níveis de escolaridade.
• Abadie e Gardeazabal (2003) consolidam as duas variáveis mais altas (school.high e school.post.high) para representar todos aqueles com educação além do ensino médio e utilizam a participação percentual de cada preditor em vez do número total de indivíduos.

• O código a seguir ilustra como consolidar essas variáveis tanto em X0 quanto em X1 e transformar os valores necessários em participações percentuais:

dataprep.out$X1["school.high",] <- dataprep.out$X1["school.high",] + 
  dataprep.out$X1["school.post.high",]
dataprep.out$X1 <- as.matrix(dataprep.out$X1[-which(rownames(dataprep.out$X1) == 
                                                      "school.post.high"),])

dataprep.out$X0["school.high",] <- dataprep.out$X0["school.high",] + 
  dataprep.out$X0["school.post.high",]
dataprep.out$X0 <- dataprep.out$X0[-which(rownames(dataprep.out$X0) == 
                                            "school.post.high"),]

lowest <- which(rownames(dataprep.out$X0) == "school.illit")
highest <- which(rownames(dataprep.out$X0) == "school.high")

dataprep.out$X1[lowest:highest,] <- 
  (100 * dataprep.out$X1[lowest:highest,]) /sum(dataprep.out$X1[lowest:highest,])
dataprep.out$X0[lowest:highest,] <- 
  100 * scale(dataprep.out$X0[lowest:highest,], 
              center = FALSE, scale = colSums(dataprep.out$X0[lowest:highest,]))

dataprep.out$X1

##                                               [,1]
## school.illit                              3.320727
## school.prim                              85.892870
## school.med                                7.522387
## school.high                               3.264015
## invest                                   24.647383
## special.gdpcap.1960.1969                  5.285468
## special.sec.agriculture.1961.1969         6.844000
## special.sec.energy.1961.1969              4.106000
## special.sec.industry.1961.1969           45.082000
## special.sec.construction.1961.1969        6.150000
## special.sec.services.venta.1961.1969     33.754000
## special.sec.services.nonventa.1961.1969   4.072000
## special.popdens.1969                    246.889999

• A função synth estima o modelo.
• A função synth.tab permite analisar os resultados

synth.out <- synth(data.prep.obj = dataprep.out, method = "BFGS")

## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.008864629 
## 
## solution.v:
##  0.01556808 0.001791073 0.04417159 0.03409436 8.45034e-05 0.2009837 0.09484593 0.007689228 0.1339499 0.008723843 0.009680725 0.1081258 0.3402913 
## 
## solution.w:
##  4.92e-08 5.17e-08 1.352e-07 2.85e-08 5.32e-08 5.177e-07 5.24e-08 7.29e-08 0.8507986 2.274e-07 4.03e-08 9.51e-08 0.1491998 5.61e-08 9.02e-08 1.061e-07

synth.tables <- synth.tab(dataprep.res = dataprep.out, synth.res = synth.out)

names(synth.tables)

## [1] "tab.pred" "tab.v"    "tab.w"    "tab.loss"

synth.tables$tab.pred

##                                         Treated Synthetic Sample Mean
## school.illit                              3.321     7.645      10.983
## school.prim                              85.893    82.285      80.911
## school.med                                7.522     6.965       5.427
## school.high                               3.264     3.105       2.679
## invest                                   24.647    21.583      21.424
## special.gdpcap.1960.1969                  5.285     5.271       3.581
## special.sec.agriculture.1961.1969         6.844     6.179      21.353
## special.sec.energy.1961.1969              4.106     2.760       5.310
## special.sec.industry.1961.1969           45.082    37.636      22.425
## special.sec.construction.1961.1969        6.150     6.952       7.276
## special.sec.services.venta.1961.1969     33.754    41.104      36.528
## special.sec.services.nonventa.1961.1969   4.072     5.371       7.111
## special.popdens.1969                    246.890   196.288      99.414

synth.tables$tab.w

##    w.weights                   unit.names unit.numbers
## 2      0.000                    Andalucia            2
## 3      0.000                       Aragon            3
## 4      0.000       Principado De Asturias            4
## 5      0.000             Baleares (Islas)            5
## 6      0.000                     Canarias            6
## 7      0.000                    Cantabria            7
## 8      0.000              Castilla Y Leon            8
## 9      0.000           Castilla-La Mancha            9
## 10     0.851                     Cataluna           10
## 11     0.000         Comunidad Valenciana           11
## 12     0.000                  Extremadura           12
## 13     0.000                      Galicia           13
## 14     0.149        Madrid (Comunidad De)           14
## 15     0.000           Murcia (Region de)           15
## 16     0.000 Navarra (Comunidad Foral De)           16
## 18     0.000                   Rioja (La)           18

• A função path.plot cria um gráfico comparando a variável de interesse na unidade tratada e na unidade sintética

path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
          Ylab = "real per-capita GDP (1986 USD, thousand)", Xlab = "year",
          Ylim = c(0, 12), Legend = c("Basque country",
                                      "synthetic Basque country"), Legend.position = "bottomright")

• A função gaps.plot mostra a diferença entra variável de interesse na unidade tratada e sintética.

gaps.plot(synth.res = synth.out, dataprep.res = dataprep.out,
          Ylab = "gap in real per-capita GDP (1986 USD, thousand)", Xlab = "year",
          Ylim = c(-1.5, 1.5), Main = NA)

• Como os resultados estão armazenados, podemos calcular o “gaps” ilustrado no gráfico anterior

gaps <- dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)

head(gaps)

##                17
## 1955  0.150218957
## 1956  0.091663836
## 1957  0.037147612
## 1958 -0.006003578
## 1959 -0.045912501
## 1960 -0.093028265

tail(gaps)

##              17
## 1992 -0.9643264
## 1993 -0.9208371
## 1994 -0.9697474
## 1995 -0.8630142
## 1996 -0.8519196
## 1997 -0.8280809

• Um benefício dos métodos de controle sintético é que eles se prestam a testes de placebo. Esses testes envolvem aplicar o método de controle sintético após reatribuir a intervenção nos dados para unidades e períodos onde a intervenção não ocorreu.
• Vamos repetir o exercício acima para Catalunha (região 10)

dataprep.out <- dataprep(
  foo = basque,
  predictors = c("school.illit", "school.prim", "school.med",
                 "school.high", "school.post.high", "invest"),
  predictors.op = "mean",
  time.predictors.prior = 1964:1969,
  special.predictors = list(
    list("gdpcap", 1960:1969 , "mean"),
    list("sec.agriculture", seq(1961, 1969, 2), "mean"),
    list("sec.energy", seq(1961, 1969, 2), "mean"),
    list("sec.industry", seq(1961, 1969, 2), "mean"),
    list("sec.construction", seq(1961, 1969, 2), "mean"),
    list("sec.services.venta", seq(1961, 1969, 2), "mean"),
    list("sec.services.nonventa", seq(1961, 1969, 2), "mean"),
    list("popdens", 1969, "mean")),
  dependent = "gdpcap",
  unit.variable = "regionno",
  unit.names.variable = "regionname",
  time.variable = "year",
  treatment.identifier = 10,
  controls.identifier = c(2:9, 11:16, 18),
  time.optimize.ssr = 1960:1969,
  time.plot = 1955:1997
)

dataprep.out$X1["school.high",] <- dataprep.out$X1["school.high",] + 
  dataprep.out$X1["school.post.high",]
dataprep.out$X1 <- as.matrix(dataprep.out$X1[-which(rownames(dataprep.out$X1) == 
                                                      "school.post.high"),])

dataprep.out$X0["school.high",] <- dataprep.out$X0["school.high",] + 
  dataprep.out$X0["school.post.high",]
dataprep.out$X0 <- dataprep.out$X0[-which(rownames(dataprep.out$X0) ==
                                            "school.post.high"),]

lowest <- which(rownames(dataprep.out$X0) == "school.illit")
highest <- which(rownames(dataprep.out$X0) == "school.high")

dataprep.out$X1[lowest:highest,] <- 
  (100 * dataprep.out$X1[lowest:highest,]) /sum(dataprep.out$X1[lowest:highest,])
dataprep.out$X0[lowest:highest,] <- 
  100 * scale(dataprep.out$X0[lowest:highest,], 
              center = FALSE, scale = colSums(dataprep.out$X0[lowest:highest,]))

synth.out <- synth(data.prep.obj = dataprep.out, method = "BFGS")

## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0003093978 
## 
## solution.v:
##  0.01405653 0.01342696 0.005037431 0.001792396 0.07939944 0.5166658 0.281674 0.08745683 9.4099e-06 1.64076e-05 1.97882e-05 0.0004392198 5.8367e-06 
## 
## solution.w:
##  1.09272e-05 7.4858e-06 0.03591013 0.2715622 1.49049e-05 0.2574583 3.0584e-06 2.667e-06 2.53151e-05 2.2189e-06 4.9246e-06 0.4349737 1.81706e-05 3.5542e-06 2.4455e-06

path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
          Ylab = "real per-capita GDP (1986 USD, thousand)", Xlab = "year",
          Ylim = c(0, 12), Legend = c("Catalonia",
                                      "synthetic Catalonia"), Legend.position = "bottomright")

• Para encerrar vamos usar um loop para aplicar o método de controle sintético para todas as regiões e fazer um gráfico com todos os “gaps”.
• Seguindo Abadie et al. (2010) as regiões com aproximações muito ruins serão retiradas (mspe, mean square prediction erros, muito alto).

ctr <- 2:18
Gaps <- matrix(0,43,18) # matriz de zeros
Gaps[,1] <- 1955:1997
mspe <- 1:17

for (i in ctr) {
  dataprep.out <- dataprep(
    foo = basque,
    predictors = c("school.illit", "school.prim", "school.med",
                   "school.high", "school.post.high", "invest"),
    predictors.op = "mean",
    time.predictors.prior = 1964:1969,
    special.predictors = list(
      list("gdpcap", 1960:1969 , "mean"),
      list("sec.agriculture", seq(1961, 1969, 2), "mean"),
      list("sec.energy", seq(1961, 1969, 2), "mean"),
      list("sec.industry", seq(1961, 1969, 2), "mean"),
      list("sec.construction", seq(1961, 1969, 2), "mean"),
      list("sec.services.venta", seq(1961, 1969, 2), "mean"),
      list("sec.services.nonventa", seq(1961, 1969, 2), "mean"),
      list("popdens", 1969, "mean")),
    dependent = "gdpcap",
    unit.variable = "regionno",
    unit.names.variable = "regionname",
    time.variable = "year",
    treatment.identifier = i,
    controls.identifier = ctr[-(i-1)],
    time.optimize.ssr = 1960:1969,
    time.plot = 1955:1997)
  
  synth.out <- synth(data.prep.obj = dataprep.out, method = "BFGS")
  Gaps[,i] <- dataprep.out$Y1plot - (dataprep.out$Y0plot %*% synth.out$solution.w)
  mspe[i-1] <- synth.out$loss.v
}

## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0004029626 
## 
## solution.v:
##  0.07386376 0.00272891 0.000207713 1.58814e-05 3.91957e-05 0.01190587 0.3607668 0.01333243 0.07576577 0.1451656 0.01529151 0.0256263 0.1005784 0.1747118 
## 
## solution.w:
##  3.111e-07 5.17e-08 3.504e-07 1.628e-07 1.042e-07 2.241e-07 0.001126068 0.002187036 0.0464731 0.1801132 0.7437621 0.02633669 1.167e-07 2.019e-07 5.84e-08 1.481e-07 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0004724927 
## 
## solution.v:
##  0.002642844 0.002284832 0.01639297 0.004553756 0.005658241 5.327e-05 0.004016764 0.3677382 0.07505774 0.3265803 6.8676e-06 0.1867112 0.008126692 0.0001762674 
## 
## solution.w:
##  0.02156546 0.1748896 0.09309231 1.67749e-05 0.007708914 0.05514567 0.0004566911 0.006421838 0.004765852 0.1315823 0.006229713 0.04396864 0.0778147 0.3702593 0.006079837 2.4633e-06 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 9.150648e-05 
## 
## solution.v:
##  0.01262663 0.0003027476 0.0007161697 0.001071258 5.1717e-06 0.04701083 0.0746889 0.1150623 0.05464994 0.07921913 0.3159775 0.00808504 0.24245 0.04813437 
## 
## solution.w:
##  1.29e-08 0.1818503 5.7792e-06 0.548692 3.03e-08 2.31e-08 1.73e-08 0.2694516 2.91e-08 1.42e-08 2.41e-08 1.4e-09 2.14e-08 2.6e-08 1.546e-07 1.45e-08 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.1029748 
## 
## solution.v:
##  0.02362584 6.75209e-05 2.46698e-05 6.49735e-05 0.02796197 0.004184091 0.6837365 1.50247e-05 0.02660968 5.76913e-05 0.002206854 0.009646694 0.2217722 2.63611e-05 
## 
## solution.w:
##  6.7843e-06 1.62849e-05 4.1984e-06 0.1309941 3.1741e-06 6.8273e-06 8.9681e-06 3.3059e-06 1.65549e-05 1.27143e-05 8.4667e-06 0.2436678 7.3661e-06 0.006267506 0.618828 0.0001480355 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.00132584 
## 
## solution.v:
##  0.0007907415 0.007949306 0.008451933 0.02155552 0.1149236 0.09864217 0.134748 0.0001619193 0.1124956 0.1302721 4.25e-08 5.3043e-06 0.3690693 0.0009345612 
## 
## solution.w:
##  4.6043e-06 3.637e-07 2.278e-07 0.176953 6.608e-07 5.8504e-06 6.14e-08 1.861e-07 4.386e-07 0.2631039 7.22e-07 3.00636e-05 0.5598997 7.12e-08 1.185e-07 4.47e-08 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0004556373 
## 
## solution.v:
##  0.0001801145 7.33009e-05 0.0417314 0.004366259 0.008976707 0.04583749 0.4423431 0.0001902437 0.1700635 7.47777e-05 7.71648e-05 0.1026893 0.04411983 0.1392767 
## 
## solution.w:
##  1.03735e-05 0.1297808 0.0001567025 0.06413347 0.000375877 9.26583e-05 0.0001075941 9.5351e-06 0.2525474 0.2462958 6.67252e-05 6.826e-07 9.51604e-05 7.6406e-06 0.3061585 0.0001610162 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0005743882 
## 
## solution.v:
##  0.1203314 0.05174992 0.00476962 0.004021745 0.0001260651 0.1878667 0.07559746 0.0576073 0.1399424 0.02228812 0.0175486 0.1135976 0.007462705 0.1970904 
## 
## solution.w:
##  7.67e-08 0.1955284 0.1544311 1.082e-07 6.07e-08 6.27e-08 0.4636655 2.7794e-06 1.527e-07 1.316e-07 0.1863652 8.549e-07 2.745e-07 1.6549e-06 1.766e-07 3.4813e-06 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.004213696 
## 
## solution.v:
##  0.177189 0.02239107 0.06514578 0.05188416 0.04243621 0.005301616 0.3965772 0.06071386 0.134774 2.43307e-05 0.004677235 0.000450996 0.0006268229 0.03780769 
## 
## solution.w:
##  0.1149473 2.94674e-05 2e-10 9.483e-06 4.8114e-06 1.21247e-05 0.03661307 5.8732e-06 1.77048e-05 0.6256898 6.62212e-05 1.6736e-06 0.2208592 4.18094e-05 6.7385e-06 0.001694781 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.001439206 
## 
## solution.v:
##  9.63019e-05 0.0001303607 0.07803746 0.004760696 0.1110528 0.009505668 0.002258596 0.1019908 0.05699398 0.1301058 0.003409267 0.2446441 0.2005504 0.05646384 
## 
## solution.w:
##  3.65762e-05 3.3466e-06 4.3727e-06 1.86e-08 7.741e-07 3.1119e-06 6.1878e-06 3.1549e-06 0.1733004 1.156e-06 3.8351e-06 0.1598184 1.2021e-06 2.557e-06 0.6668128 2.1619e-06 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.001059888 
## 
## solution.v:
##  0.07334557 1.59812e-05 0.05794287 0.06598513 3.1532e-06 0.1520786 0.124526 1.44564e-05 0.1543537 0.04329191 0.01525443 0.005884243 6.949e-07 0.3073033 
## 
## solution.w:
##  0.225633 7.2265e-06 7.162e-07 0.0005551303 8.4643e-06 0.4553902 1.9966e-05 0.08330879 0.1129504 1.06115e-05 4.9358e-06 0.1218524 2.0811e-06 9.7441e-06 1e-10 0.0002463404 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.114639 
## 
## solution.v:
##  0.003627175 0.02519924 0.06378065 0.01481177 0.03098362 0.004408407 0.139872 0.4151944 0.1211422 0.04577488 0.0001632541 0.04698581 0.02180939 0.06624725 
## 
## solution.w:
##  5.86e-08 6.6e-09 6e-10 3.01e-08 4.2e-08 7.9e-09 3.4e-09 0.9999998 2.4e-09 1.07e-08 2.41e-08 7.4e-09 4.7e-09 1.5e-08 2.2e-09 1.95e-08 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0005603041 
## 
## solution.v:
##  2.78728e-05 0.1164736 0.04072458 0.03835988 0.04222253 5.11271e-05 0.4046217 0.1316852 0.05491416 0.09097371 0.0231572 0.05678729 1.033e-06 9.38e-08 
## 
## solution.w:
##  0.3524079 0.002318552 0.04825064 0.001381329 0.0317231 0.002445645 4.0802e-06 0.4452373 0.001372995 0.002610457 0.03595204 0.0005092839 0.07133017 0.00165655 0.001216286 0.001583646 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.5308468 
## 
## solution.v:
##  0.1915117 0.01258113 0.006014191 0.01493367 0.02118893 0.01335143 0.2489348 0.1472647 0.01235817 0.03482911 0.009433851 0.0002934643 0.01595775 0.2713471 
## 
## solution.w:
##  2.529e-07 5e-09 8.3e-09 1.22e-08 6.52e-08 2.4e-09 8.9e-09 5.9e-09 5.244e-07 3.17e-08 6.9e-09 1.93e-08 9.8e-09 3.8e-09 0.999999 3.6e-09 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.001468835 
## 
## solution.v:
##  3.73216e-05 5.08413e-05 0.06502392 0.0003014844 0.05564236 0.001601267 0.4329044 1.11608e-05 0.04089693 0.0001131968 0.01435637 0.0006470147 0.2737504 0.1146633 
## 
## solution.w:
##  9.29e-08 1.9609e-06 0.0008957282 1.69e-07 0.5441035 3.644e-07 0.2612477 0.1432925 1.385e-07 2.103e-07 7e-09 3.501e-07 2.28e-08 2.41e-08 2.016e-07 0.050457 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0002624747 
## 
## solution.v:
##  0.01433404 0.003357316 0.0002533356 0.0002779596 0.0002326023 0.01556371 0.003633353 0.3762792 0.007697547 0.3426853 0.002334471 0.2023423 0.02959998 0.001408978 
## 
## solution.w:
##  0.0003640501 0.1859778 0.008069571 1.3572e-06 6.15178e-05 0.03890214 0.0006712699 1.003e-07 6.13724e-05 5.34634e-05 0.002186609 0.0001337051 0.01722144 0.0002545825 0.1170119 0.6290291 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.008864606 
## 
## solution.v:
##  0.02773094 1.194e-07 1.60609e-05 0.0007163836 1.486e-07 0.002423908 0.0587055 0.2651997 0.02851006 0.291276 0.007994382 0.004053188 0.009398579 0.303975 
## 
## solution.w:
##  2.53e-08 4.63e-08 6.44e-08 2.81e-08 3.37e-08 4.844e-07 4.2e-08 4.69e-08 0.8508145 9.75e-08 3.2e-08 5.54e-08 0.1491843 4.86e-08 9.89e-08 1.162e-07 
## 
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0007161077 
## 
## solution.v:
##  0.1276237 0.114429 0.07064796 0.08790274 0.07068704 0.01296239 0.119247 0.01685659 0.07696669 0.08874159 0.02770201 0.05267171 0.06492279 0.06863878 
## 
## solution.w:
##  1.6038e-06 7.0426e-06 1.20374e-05 5.4303e-06 9.7175e-06 4.9246e-06 1.00274e-05 0.1048676 3.4142e-06 7.9418e-06 4.66e-08 8.6145e-06 2.9816e-06 0.0008517394 0.8941913 1.55603e-05

mspe

##  [1] 4.029626e-04 4.724927e-04 9.150648e-05 1.029748e-01 1.325840e-03
##  [6] 4.556373e-04 5.743882e-04 4.213696e-03 1.439206e-03 1.059888e-03
## [11] 1.146390e-01 5.603041e-04 5.308468e-01 1.468835e-03 2.624747e-04
## [16] 8.864606e-03 7.161077e-04

mspe1 <- mspe/mspe[16]

Gaps1 <- Gaps[1:43,]
ctr.p <- c(2:4, 6:11, 13,15,16,18)

Gaps1 <- Gaps[1:43,]
ctr.p <- c(2:4, 6:11, 13,15,16,18)

plot(Gaps1[,1],Gaps1[,17], type="l", lwd=5, ylim = c(-2,2),
     xlab = " ",
     ylab = "Diferença entre observado e unidade sintética",
     main="Teste de Placebo\nDiferença entre observado e unidade sintética, país Basco em destaque")
for(i in ctr.p) {
  lines(Gaps1[,1],Gaps1[,i], lwd=2, col=16)
}
abline(v=1970, lty=2)
abline(h=0, lty=2)