Para resolução das questões da lista é necessário instalar e carregar os seguintes pacotes:

library(tidyverse)
library(wooldridge)
library(AER)
library(car)
library(estimatr)
library(sampleSelection)
library(plm)
library(modelsummary)

Questão 1

Utilize a base de dados crime1 do pacote wooldridge para esta questão.

  1. Defina uma variável binária, chamada \(arr86\), igual a 1 se um homem foi preso pelo menos uma vez durante 1986, e zero caso contrário. Estime um Modelo de Probabilidade Linear (MPL) relacionando \(arr86\) a \(pcnv\), \(avgsen\), \(tottime\), \(ptime86\), \(inc86\), \(black\), \(hispan\) e \(born60\). Reporte os erros-padrão usuais e robustos à heterocedasticidade. Qual é o efeito estimado na probabilidade de prisão se \(pcnv\) for de 0.25 para 0.75?

  2. Teste a significância conjunta de \(avgsen\) e \(tottime\), usando um teste não robusto e um teste robusto.

  3. Agora estime o modelo por probit. Nos valores médios de \(avgsen\), \(tottime\), \(inc86\) e \(ptime86\) na amostra, e com \(black = 1\), \(hispan = 0\) e \(born60 = 1\), qual é o efeito estimado na probabilidade de prisão se \(pcnv\) for de 0.25 para 0.75? Compare este resultado com a resposta da parte a.

  4. Para o modelo probit estimado na parte c, obtenha o percentual corretamente previsto. Qual é o percentual corretamente previsto quando \(narr86 = 0\)? \(Quando narr86 = 1\)? Comente os resultados obtidos.

  5. No modelo probit, adicione os termos \(pcnv^2\), \(ptime86^2\) e \(inc86^2\) ao modelo. Estes são individualmente ou conjuntamente significativos? Descreva a relação estimada entre a probabilidade de prisão e \(pcnv\). Em particular, em que ponto a probabilidade de condenação tem um efeito negativo sobre a probabilidade de prisão?

Solução

data(crime1)
  1. Para estimar os modelos de probabalidade linear com erros-padrão não robustos e robustos à heterocedasticidade, utiliza-se o comando lm_robust do pacote estimatr, pois o mesmo já informa os resultados desejados no output da estimação principal, sem a necessidade de realizar algum procedimento de correção após a estimação. Assim, o modelo com erros não robustos pode ser estimado da seguinte forma:
crime1 <- crime1 |>  
  mutate(arr86=ifelse(narr86>0,1,0))

Q1_formula <- as.formula("arr86 ~ pcnv + avgsen + tottime + ptime86 + 
                            inc86 + black + hispan + born60")
Q1_LPM <- lm_robust(Q1_formula, data = crime1, se_type = "classical")

Para estimar o modelo com erros robustos à heterocedasticidade é necessário alterar o argumento se_type de classical para HC1 ou stata. Pode-se verificar as diferentes formulas para o cálculo dos erros-padrão aqui

Q1_LPM.2 <- update(Q1_LPM, se_type="stata")

Os resultados dos modelos estimados podem ser verificados na tabela abaixo:

LPM 1 LPM 2
(Intercept) 0.361 0.361
(0.016) (0.017)
pcnv -0.154 -0.154
(0.021) (0.019)
avgsen 0.004 0.004
(0.006) (0.006)
tottime -0.002 -0.002
(0.005) (0.004)
ptime86 -0.022 -0.022
(0.004) (0.003)
inc86 -0.001 -0.001
(0.000) (0.000)
black 0.162 0.162
(0.024) (0.026)
hispan 0.089 0.089
(0.021) (0.021)
born60 0.003 0.003
(0.017) (0.017)
Num.Obs. 2725 2725
R2 0.082 0.082
R2 Adj. 0.080 0.080
AIC 3137.3 3137.3
BIC 3196.4 3196.4
RMSE 0.43 0.43

os erros-padrão robustos nas variáveis não demográficas são frequentemente menores do que os erros-padrão usuais. A significância estatística dos coeficientes das regressões estimadas é a mesma usando qualquer um dos conjuntos de erros-padrão.

Quando a proporção de condenações anteriores (\(pcnv\)) vai de 0.25 para 0.75 a probabilidade estimada de ser preso cai cerca de 7,7 pontos percentuais, como pode ser visto abaixo:

.5*coef(Q1_LPM)["pcnv"]
##       pcnv 
## -0.0771901
  1. Para testar a significância conjunta das variáveis \(avgsen\) e \(tottime\) será utilizada função linearHypothesis do pacote car.
linearHypothesis(Q1_LPM, c("avgsen=0","tottime=0"), test = "F")
## Linear hypothesis test
## 
## Hypothesis:
## avgsen = 0
## tottime = 0
## 
## Model 1: restricted model
## Model 2: arr86 ~ pcnv + avgsen + tottime + ptime86 + inc86 + black + hispan + 
##     born60
## 
##   Res.Df Df      F Pr(>F)
## 1   2718                 
## 2   2716  2 0.1791  0.836
linearHypothesis(Q1_LPM.2, c("avgsen=0","tottime=0"), test = "F")
## Linear hypothesis test
## 
## Hypothesis:
## avgsen = 0
## tottime = 0
## 
## Model 1: restricted model
## Model 2: arr86 ~ pcnv + avgsen + tottime + ptime86 + inc86 + black + hispan + 
##     born60
## 
##   Res.Df Df      F Pr(>F)
## 1   2718                 
## 2   2716  2 0.1839  0.832

Os resultados são virtualmente idênticos independentemente da forma de cálculo do erro-padrão, indicando que não se pode rejeitar a hipotese nula não possuirem significância estatistica conjuntamente.

  1. O modelo probit é estimado da seguinte forma:
Q1_probit <- glm(Q1_formula, family = binomial(link = "probit"), data = crime1)

Os resultados são comparados com aqueles obtidos através do LPM:

LPM 1 LPM 2 PROBIT
(Intercept) 0.361 0.361 -0.314
(0.016) (0.017) (0.051)
pcnv -0.154 -0.154 -0.553
(0.021) (0.019) (0.072)
avgsen 0.004 0.004 0.013
(0.006) (0.006) (0.021)
tottime -0.002 -0.002 -0.008
(0.005) (0.004) (0.017)
ptime86 -0.022 -0.022 -0.081
(0.004) (0.003) (0.017)
inc86 -0.001 -0.001 -0.005
(0.000) (0.000) (0.000)
black 0.162 0.162 0.467
(0.024) (0.026) (0.072)
hispan 0.089 0.089 0.291
(0.021) (0.021) (0.065)
born60 0.003 0.003 0.011
(0.017) (0.017) (0.056)
Num.Obs. 2725 2725 2725
R2 0.082 0.082
R2 Adj. 0.080 0.080
AIC 3137.3 3137.3 2985.3
BIC 3196.4 3196.4 3038.5
Log.Lik. -1483.641
F 27.446
RMSE 0.43 0.43 0.43

No modelo probit, o efeito na probabilidade de prisão quando (\(pcnv\)) vai de 0.25 para 0.75 é obtido da seguinte forma:

efct <- predict(Q1_probit, newdata = data.frame(
  "black"=1,
  "hispan"=0,
  "born60"=1,
  "avgsen"=mean(crime1$avgsen),
  "tottime"=mean(crime1$tottime),
  "inc86"=mean(crime1$inc86),
  "ptime86"=mean(crime1$ptime86),
  "pcnv"=c(0.25,0.75)),
  type = "response")

diff(efct)
##          2 
## -0.1016607

O resultado acima mostra que a probabilidade cai em cerca de 0,102 o que é um pouco maior do que o efeito obtido no LPM.

  1. Os percentuais corretamente previstos (total e quando \(narr86 = 0\)? e \(narr86 = 1\)? ) são calculados com base nos valores preditos de Q1_probit. Abaixo é apresentado como calcular tais percentuais manualmente e também utilizando o comando ConfusionMatrix do pacote caret
probit_predicted <- predict(Q1_probit, type = "response")
probit_class <- ifelse(probit_predicted>=0.5, 1,0)

cm <- table(Previsto=probit_class,Real=crime1$arr86)

total <- colSums(cm)
cm <- rbind(cm, total)

total <- rowSums(cm)
cm <- cbind(cm, total)

cm
##          0   1 total
## 0     1903 677  2580
## 1       67  78   145
## total 1970 755  2725
cm.2 <- caret::confusionMatrix(factor(probit_class),factor(crime1$arr86))
cm.2
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1903  677
##          1   67   78
##                                           
##                Accuracy : 0.727           
##                  95% CI : (0.7098, 0.7436)
##     No Information Rate : 0.7229          
##     P-Value [Acc > NIR] : 0.3275          
##                                           
##                   Kappa : 0.0923          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.9660          
##             Specificity : 0.1033          
##          Pos Pred Value : 0.7376          
##          Neg Pred Value : 0.5379          
##              Prevalence : 0.7229          
##          Detection Rate : 0.6983          
##    Detection Prevalence : 0.9468          
##       Balanced Accuracy : 0.5347          
##                                           
##        'Positive' Class : 0               
##

Para os homens que não foram presos, o modelo probit acerta cerca de 96,6% das vezes. Para os homens que foram presos, o modelo probit acerta apenas cerca de 10,3% das vezes. A porcentagem total de previsões corretas é consideravalmente alta 72,7%.

  1. O modelo com termos quadráticos das variáveis \(pcnv\), \(ptime86\) e \(inc86\) é estimado como se segue e seus resultados são comparados com as demais estimações na tabela abaixo:
Q1_probit.2 <-update(Q1_probit, . ~ . + pcnvsq + pt86sq + inc86sq)

tab.3 <- modelsummary(list("LPM 1"=Q1_LPM,
                           "LPM 2"=Q1_LPM.2,
                           "PROBIT"=Q1_probit,
                           "PROBIT 2"=Q1_probit.2), output = "gt")
tab.3
LPM 1 LPM 2 PROBIT PROBIT 2
(Intercept) 0.361 0.361 -0.314 -0.337
(0.016) (0.017) (0.051) (0.056)
pcnv -0.154 -0.154 -0.553 0.217
(0.021) (0.019) (0.072) (0.259)
avgsen 0.004 0.004 0.013 0.014
(0.006) (0.006) (0.021) (0.026)
tottime -0.002 -0.002 -0.008 -0.018
(0.005) (0.004) (0.017) (0.021)
ptime86 -0.022 -0.022 -0.081 0.745
(0.004) (0.003) (0.017) (0.144)
inc86 -0.001 -0.001 -0.005 -0.006
(0.000) (0.000) (0.000) (0.001)
black 0.162 0.162 0.467 0.437
(0.024) (0.026) (0.072) (0.073)
hispan 0.089 0.089 0.291 0.266
(0.021) (0.021) (0.065) (0.067)
born60 0.003 0.003 0.011 -0.015
(0.017) (0.017) (0.056) (0.057)
pcnvsq -0.857
(0.270)
pt86sq -0.104
(0.022)
inc86sq 0.000
(0.000)
Num.Obs. 2725 2725 2725 2725
R2 0.082 0.082
R2 Adj. 0.080 0.080
AIC 3137.3 3137.3 2985.3 2903.6
BIC 3196.4 3196.4 3038.5 2974.5
Log.Lik. -1483.641 -1439.800
F 27.446 23.228
RMSE 0.43 0.43 0.43 0.42

Os termos quadráticos são significativos individualmente e conjuntamente, como pode ser visto no teste de significância conjunta abaixo:

linearHypothesis(Q1_probit.2, c("pcnvsq=0","pt86sq=0","inc86sq=0"), test = "F")
## Linear hypothesis test
## 
## Hypothesis:
## pcnvsq = 0
## pt86sq = 0
## inc86sq = 0
## 
## Model 1: restricted model
## Model 2: arr86 ~ pcnv + avgsen + tottime + ptime86 + inc86 + black + hispan + 
##     born60 + pcnvsq + pt86sq + inc86sq
## 
##   Res.Df Df      F    Pr(>F)    
## 1   2716                        
## 2   2713  3 12.985 2.023e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ao inserir o termo quadrático na regressão temos que em baixos níveis de de \(pcnv\) existe na verdade uma relação positiva entre a probabilidade de prisão e \(pcnv\), o que não faz muito sentido. O ponto de inflexão dessa relação ocorre em \(pcnv = 0.127\), conforme calculado abaixo:

abs(Q1_probit.2$coefficients["pcnv"]/(2*Q1_probit.2$coefficients["pcnvsq"]))
##      pcnv 
## 0.1264604

Questão 2

A base de dados ShipAccidents disponível no pacote AER contém o número de acidentes por de meses de serviço (incidents) para uma determinada amostra de navios. Adicionalmente, são informados a quantidade agregada de meses de serviço de um navio (service), seu tipo (type), o período em que foi construído (construction) e o período de operação (operation). Faça o que se pede:

  1. Estime um modelo de regressão Poisson nos seguintes termos: \[\ln E [incidents] = \mathbf{x^{'}}\beta.\] Considere as seguintes situações: (i) modelo com o conjunto completo de variáveis explicativas; (ii) modelo sem considerar os efeitos da variável \(type\) e (iii) modelo sem a inclusão de $construction$.

  2. Realize os testes de hipotese adequados para avaliar se as restrições de exclusão aplicadas no item anterior são válidas.

  3. Teste se ocorre sobredispersão no modelo estimado.

Solução

data("ShipAccidents")

ShipAccidents <- ShipAccidents |> 
  filter(service>0)
  1. Os parâmetros estimados dos modelos de poisson descritos no enunciado são reportados na tabela abaixo:
Q2_formula<- as.formula("incidents ~ type + construction + operation")

Q2_poisson <- glm(Q2_formula, family = poisson, data = ShipAccidents,
                  offset = log(service))

Q2_poisson.2 <- update(Q2_poisson, . ~ . - type)

Q2_poisson.3 <- update(Q2_poisson, . ~ . - construction)
Poisson 1 Poisson 2 Poisson 3
(Intercept) -6.403 -6.947 -5.800
(0.218) (0.127) (0.178)
typeB -0.545 -0.744
(0.178) (0.169)
typeC -0.689 -0.755
(0.329) (0.328)
typeD -0.074 -0.184
(0.291) (0.288)
typeE 0.321 0.384
(0.236) (0.235)
construction1965-69 0.696 0.754
(0.150) (0.149)
construction1970-74 0.817 1.050
(0.170) (0.158)
construction1975-79 0.445 0.700
(0.233) (0.220)
operation1975-79 0.384 0.387 0.500
(0.118) (0.118) (0.112)
Num.Obs. 34 34 34
AIC 154.8 170.4 180.2
BIC 168.6 178.0 189.4
Log.Lik. -68.415 -80.201 -84.115
F 13.430 18.784 17.351
RMSE 2.55 3.70 6.46
  1. Para testar se as restrições de exclusão dos modelos 2 e 3 são válidas, será utilizado o teste da Razão de Verossimilhança (RV).
# Opção 1 - Calculando o Log-likelihood para cada modelo estimado

logL_M1 <- logLik(Q2_poisson)
logL_M2 <- logLik(Q2_poisson.2)
logL_M3 <- logLik(Q2_poisson.3)

rv_1 <- 2*(logL_M1-logL_M2)
rv_1 
## 'log Lik.' 23.57334 (df=9)
qchisq(df=3, p=.05, lower.tail = F)
## [1] 7.814728
rv_2 <- 2*(logL_M1-logL_M3)
rv_2 
## 'log Lik.' 31.40118 (df=9)
qchisq(df=4, p=.05, lower.tail = F)
## [1] 9.487729
# Opção 2 - Utilizar a função lrtest do pacote lmtest

lmtest::lrtest(Q2_poisson,Q2_poisson.2)
## Likelihood ratio test
## 
## Model 1: incidents ~ type + construction + operation
## Model 2: incidents ~ construction + operation
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   9 -68.415                         
## 2   5 -80.201 -4 23.573  9.725e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lmtest::lrtest(Q2_poisson, Q2_poisson.3)
## Likelihood ratio test
## 
## Model 1: incidents ~ type + construction + operation
## Model 2: incidents ~ type + operation
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   9 -68.415                         
## 2   6 -84.115 -3 31.401  6.998e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As hipóteses que o efeito do tipo do navio e do período de construção são irrelevantes para a estimação são rejeitadas.

  1. Para testar a ocorrência de sobredispersão serão aplicados os testes de Cameron e Trivedi (1990) e o teste do Multiplicador de Lagrange.
# Cameron e Trivedi (1990) Test

mu <- predict(Q2_poisson, type = "response")
z <- ((ShipAccidents$incidents - mu)^2 - ShipAccidents$incidents)/(mu*sqrt(2))
summary(zscore <- lm(z~ 1))
## 
## Call:
## lm(formula = z ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2698 -0.7757 -0.3301  0.3837  6.7109 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.2248     0.2406   0.934    0.357
## 
## Residual standard error: 1.403 on 33 degrees of freedom
summary(zscore <- lm(z~ 0 + mu^2/mu))
## 
## Call:
## lm(formula = z ~ 0 + mu^2/mu)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0408 -0.4059  0.0804  0.6193  6.9474 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)
## mu -0.00791    0.01291  -0.613    0.544
## 
## Residual standard error: 1.413 on 33 degrees of freedom
## Multiple R-squared:  0.01126,    Adjusted R-squared:  -0.01871 
## F-statistic: 0.3756 on 1 and 33 DF,  p-value: 0.5441
dispersiontest(Q2_poisson)
## 
##  Overdispersion test
## 
## data:  Q2_poisson
## z = 0.93429, p-value = 0.1751
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   1.317918
dispersiontest(Q2_poisson, trafo = 2)
## 
##  Overdispersion test
## 
## data:  Q2_poisson
## z = -0.6129, p-value = 0.73
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
##      alpha 
## -0.0111868
#LM Test
(sum((ShipAccidents$incidents - mu)^2-ShipAccidents$incidents)/sqrt(2*sum(mu^2)))^2
## [1] 0.7504421
qchisq(df=1, p=.05, lower.tail = F)
## [1] 3.841459

As estatísticas t para as duas regressões nos testes baseados em regressão de Cameron e Trivedi são 0,934 e -0,613, respectivamente, então, com base nesses testes, não rejeitamos H0: ausência de sobredispersão. A estatística LM para a mesma hipótese é 0,75044 com um grau de liberdade. O valor crítico da tabela é 3,84, então, novamente, a hipótese do modelo de Poisson não é rejeitada. No entanto, o teste de momento condicional é contraditório;então a hipótese agora é rejeitada. Este teste é muito mais geral, pois a forma da sobredispersão não é especificada, o que pode explicar a diferença.

Questão 3

Utilize o dataframe smoke para responder os itens abaixo:

  1. Utilize um modelo de regressão linear para explicar \(cigs\), o número de cigarros fumados por dia. Use como variáveis explicativas \(log(cigpric)\), \(log(income)\), \(restaurn\), \(white\), \(educ\), \(age\) e \(age2\). As variáveis de preço e renda são significativas? O uso de erros-padrão robustos à heterocedasticidade altera suas conclusões?

  2. Agora estime um modelo de regressão Poisson para \(cigs\), com uma média condicional exponencial e as mesmas variáveis explicativas do item a. Usando os erros-padrão usuais do método de Máxima Verossimilhança (MLE), as variáveis de preço e renda são significativas ao nível de 5%? Interprete seus coeficientes.

  3. Há evidência de sobredispersão?

Solução

data("smoke")

Q3.formula <- as.formula("cigs ~ lcigpric+lincome + restaurn+
                         white + educ + age + agesq")
  1. A estimação por OLS com erros-padrão robustos e não robustos apresenta os seguintes resultados:
Q3.lm <- lm_robust(Q3.formula,data = smoke, se_type = "classical")

Q3.lm.robust <- update(Q3.lm, se_type="stata")
OLS 1 OLS 2
(Intercept) -2.682 -2.682
(24.221) (25.902)
lcigpric -0.851 -0.851
(5.782) (6.054)
lincome 0.869 0.869
(0.729) (0.598)
restaurn -2.866 -2.866
(1.117) (1.017)
white -0.559 -0.559
(1.459) (1.378)
educ -0.502 -0.502
(0.167) (0.162)
age 0.775 0.775
(0.161) (0.138)
agesq -0.009 -0.009
(0.002) (0.001)
Num.Obs. 807 807
R2 0.053 0.053
R2 Adj. 0.045 0.045
AIC 6490.3 6490.3
BIC 6532.5 6532.5
RMSE 13.35 13.35

Nem a variável preço nem a variável renda são significativas em nenhum nível razoável de significância, embora as estimativas dos coeficientes sejam do sinal esperado. Não importa se usamos os erros-padrão usuais ou robustos.

  1. Como pode ser visto abaixo, ao estimar o modelo de regressão poisson a variável de preço ainda é insignificante (valor-p = 0,46), a variável de renda, com base nos erros-padrão Poisson usuais, é muito significativa: t = 5,11. Ambas as estimativas são elasticidades: a elasticidade estimada do preço é -0,106 e a elasticidade estimada da renda é 0,104.
Q3.poisson <- glm(Q3.formula, family = poisson, data = smoke)
summary(Q3.poisson)
## 
## Call:
## glm(formula = Q3.formula, family = poisson, data = smoke)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -6.329  -4.224  -3.275   2.245  13.976  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.964e-01  6.139e-01   0.646    0.518    
## lcigpric    -1.060e-01  1.434e-01  -0.739    0.460    
## lincome      1.037e-01  2.028e-02   5.115 3.14e-07 ***
## restaurn    -3.636e-01  3.122e-02 -11.646  < 2e-16 ***
## white       -5.520e-02  3.742e-02  -1.475    0.140    
## educ        -5.942e-02  4.256e-03 -13.961  < 2e-16 ***
## age          1.143e-01  4.969e-03  22.994  < 2e-16 ***
## agesq       -1.371e-03  5.695e-05 -24.070  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15821  on 806  degrees of freedom
## Residual deviance: 14752  on 799  degrees of freedom
## AIC: 16239
## 
## Number of Fisher Scoring iterations: 6
dispersiontest(Q3.poisson)
## 
##  Overdispersion test
## 
## data:  Q3.poisson
## z = 13.097, p-value < 2.2e-16
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   20.12705
dispersiontest(Q3.poisson, trafo = 2)
## 
##  Overdispersion test
## 
## data:  Q3.poisson
## z = 12.429, p-value < 2.2e-16
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
##    alpha 
## 1.970601

Questão 4

Replique a Tabela 22.7 do livro texto ECONOMETRIC ANALYSIS, p. 786

Solução

data(PSID1976)

PSID1976 <- PSID1976 |> 
  mutate(kids=factor(ifelse(youngkids | oldkids >0,"yes","no")))
LFP_formula <- as.formula(participation ~ age + I(age^2) + fincome + education + kids)

Wage_formula <- as.formula(wage ~ experience + I(experience^2) + education + city)

Q4_OLS <- lm_robust(Wage_formula, data = PSID1976 |> filter(participation=="yes"),
                    se_type = "classical")

Q4_ML <- selection(LFP_formula, Wage_formula, data = PSID1976, method = "ml")

Q4_heck <- update(Q4_ML, method="2step")

As estimativas de ML estão de acordo com a correção do texto realizada.

Questão 5

Considere o seguinte modelo para estimar os efeitos de várias variáveis, incluindo o consumo de cigarros, no peso dos recém-nascidos: \[ \log(bwgth) = \beta_{0} + \beta_{1}male + \beta_{2}parity + \beta_{3}\log{faminc} + \beta_{packs} + u.\]Onde \(male\) é um indicador binário igual a um se a criança for do sexo masculino, \(parity\) é a ordem de nascimento desta criança, \(faminc\) é a renda familiar e \(packs\) é o número médio de maços de cigarros fumados por dia durante a gravidez.

  1. Por que você poderia esperar que \(packs\) esteja correlacionado com \(u\)?

  2. Suponha que você tenha dados sobre o preço médio dos cigarros no estado de residência de cada mulher. Discuta se essa informação provavelmente satisfaz as propriedades de uma boa variável instrumental para \(packs\).

  3. Use os dados em bwght para estimar a equação do enunciado. Primeiro, use Mínimos Quadrados Ordinários (OLS). Em seguida, use 2SLS, onde \(cigprice\) é um instrumento para \(packs\). Discuta quaisquer diferenças importantes nas estimativas OLS e 2SLS.

  4. Estime a forma reduzida para \(packs\). O que você conclui sobre a identificação do modelo estimado usando \(cigprice\) como instrumento para \(packs\)? Que impacto essa conclusão tem em sua resposta da parte c?

Solução

data(bwght)

Q5_formula <- as.formula("lbwght~ male + parity + lfaminc + packs")
IV_formula <- as.formula("lbwght ~ male + parity + lfaminc + packs | cigprice + male + parity + lfaminc")
red_formula <- as.formula("packs ~ male + parity + lfaminc + cigprice")

  1. Pode haver fatores de saúde não observados correlacionados com o comportamento de fumar que afetam o peso ao nascer do bebê. Por exemplo, mulheres que fumam durante a gravidez podem, em média, consumir mais café ou álcool, ou ter uma alimentação menos nutritiva.

  2. A economia básica diz que os maços de cigarro devem estar negativamente correlacionados com o preço do cigarro, embora a correlação possa ser pequena (especialmente porque o preço é agregado ao nível estadual). À primeira vista, parece que o preço do cigarro deve ser exógeno na equação do enunciado, mas devemos ter um pouco de cuidado. Um componente do preço do cigarro é o imposto estadual sobre cigarros. Estados que têm impostos mais baixos sobre cigarros também podem ter uma qualidade mais baixa dos cuidados de saúde, em média. A qualidade dos cuidados de saúde está em \(u\), e assim talvez o preço do cigarro não atenda ao requisito de exogeneidade para uma IV.

Q5_OLS <- lm_robust(Q5_formula, data = bwght, se_type = "classical")

Q5_2SLS <-ivreg(formula = IV_formula, data = bwght)
OLS 2SLS
(Intercept) 4.676 4.468
(0.022) (0.259)
male 0.026 0.030
(0.010) (0.018)
parity 0.015 -0.001
(0.006) (0.022)
lfaminc 0.018 0.064
(0.006) (0.057)
packs -0.084 0.797
(0.017) (1.086)
Num.Obs. 1388 1388
R2 0.035 -1.812
R2 Adj. 0.032 -1.820
AIC -700.1 784.4
BIC -668.7 815.8
RMSE 0.19 0.32

A diferença entre OLS e IV no efeito estimado de pacotes no peso ao nascer é enorme. Com a estimativa de OLS, um pacote adicional de cigarros é estimado para reduzir o peso ao nascer em cerca de 8,4%, e é estatisticamente significativo. A estimativa de IV tem o sinal oposto, é de grande magnitude e não é estatisticamente significativa. O sinal e o tamanho do efeito do fumo não são realistas.

Q5_red.f <- lm_robust(red_formula, data = bwght, se_type = "classical" )
summary(Q5_red.f)
## 
## Call:
## lm_robust(formula = red_formula, data = bwght, se_type = "classical")
## 
## Standard error type:  classical 
## 
## Coefficients:
##              Estimate Std. Error t value  Pr(>|t|)   CI Lower CI Upper   DF
## (Intercept)  0.137408  0.1040005  1.3212 1.866e-01 -0.0666084  0.34142 1383
## male        -0.004726  0.0158539 -0.2981 7.657e-01 -0.0358264  0.02637 1383
## parity       0.018149  0.0088802  2.0438 4.116e-02  0.0007291  0.03557 1383
## lfaminc     -0.052637  0.0086991 -6.0509 1.852e-09 -0.0697023 -0.03557 1383
## cigprice     0.000777  0.0007763  1.0009 3.171e-01 -0.0007459  0.00230 1383
## 
## Multiple R-squared:  0.03045 ,   Adjusted R-squared:  0.02765 
## F-statistic: 10.86 on 4 and 1383 DF,  p-value: 1.137e-08

As estimativas da forma reduzida mostram que o preço do cigarro não afeta significativamente \(packs\). Na verdade, o coeficiente sobre o preço do cigarro não tem o sinal que esperamos. Portanto, o preço do cigarro falha como um IV para os pacotes porque o preço do cigarro não está parcialmente correlacionado com os pacotes com um sinal adequado para a correlação.

Questão 6

Para os itens abaixo, utilize os dados base card.

  1. Estime uma equação de \(log(wage)\), por OLS, com \(educ\), \(exper\), \(exper2\), \(black\), \(south\), \(smsa\), \(reg661\) até \(reg668\) e \(smsa66\) como variáveis explicativas. Compare seus resultados com a Tabela 2, Coluna (2) em Card (1995).

  2. Estime uma equação de forma reduzida para \(educ\) contendo todas as variáveis explicativas da parte a e a variável dummy \(nearc4\). \(Educ\) e \(nearc4\) têm uma correlação parcial estatisticamente significativa?

  3. Estime a equação \(log(wage)\) por IV, usando \(nearc4\) como instrumento para \(educ\). Compare o intervalo de confiança de 95% para o retorno à educação com aquele obtido da parte a.

  4. Agora use \(nearc2\) juntamente com \(nearc4\) como instrumentos para \(educ\). Primeiro, estime a forma reduzida para \(educ\) e comente se \(nearc2\) ou \(nearc4\) está mais fortemente relacionado com \(educ\). Como as estimativas 2SLS se comparam com as estimativas anteriores?

Questão 7

Use os dados em crime4 para esta questão.

  1. Considere um modelo que relaciona a taxa de criminalidade ($ crmrte$) a \(prbarr\), \(prbconv\), \(prbpris\), \(avgsen\) e \(polpc\), em que todas as variáveis estão em sua forma logarítmica. Alem disso, inclua variáveis dummy para cada ano da amostra. Estime uma versão com Efeitos Aleatórios (RE) e uma com Efeitos Fixos (FE) do modelo. Utilize o Teste de Hausman para compararar as estimações com RE e FE.

  2. Adicione as variáveis salariais (em forma logarítmica) e teste a significância conjunta após a estimação pelo método de Efeitos Fixos (FE).

  3. Estime a equação pelo método de Primeira Diferença (FD) e comente sobre quaisquer mudanças notáveis. Os erros-padrão mudam muito entre FE e FD?

Solução

d_vars <- paste(paste0("d8", 2:7), collapse = " + ")
formula_string.q7 <- paste("lcrmrte ~ lprbarr + lprbconv + lprbpris + lavgsen + lpolpc +",d_vars)
Q7_formula <- as.formula(formula_string.q7)
  1. Embora haja diferenças nas estimativas de efeitos aleatórios (RE) e de efeitos fixos (FE), os sinais são os mesmos e as magnitudes são semelhantes. O teste de Hausman oferece suporte para rejeição da suposição de RE. Portanto, evemos preferir as estimativas de efeitos fixos (FE).
Q7_RE <- plm(Q7_formula, model="random", data = crime4)

G <- length(unique(crime4$county))
c <- G/(G - 1)

coeftest(Q7_RE,c*vcovHC(Q7_RE,type="HC1",cluster="group"))
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -1.6726322  0.5683384 -2.9430 0.0033722 ** 
## lprbarr     -0.4252097  0.0629647 -6.7531 3.342e-11 ***
## lprbconv    -0.3271464  0.0499984 -6.5431 1.265e-10 ***
## lprbpris    -0.1793507  0.0457911 -3.9167 9.977e-05 ***
## lavgsen     -0.0083696  0.0322633 -0.2594 0.7954022    
## lpolpc       0.4294148  0.0879358  4.8833 1.330e-06 ***
## d82          0.0137442  0.0164988  0.8330 0.4051412    
## d83         -0.0753880  0.0194987 -3.8663 0.0001222 ***
## d84         -0.1130975  0.0217197 -5.2071 2.615e-07 ***
## d85         -0.1057261  0.0254789 -4.1495 3.798e-05 ***
## d86         -0.0795307  0.0239331 -3.3230 0.0009430 ***
## d87         -0.0424581  0.0246604 -1.7217 0.0856219 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Q7_FE <- update(Q7_RE, model="within")

coeftest(Q7_FE,c*vcovHC(Q7_FE,type="HC1",cluster="group"))
## 
## t test of coefficients:
## 
##            Estimate Std. Error t value  Pr(>|t|)    
## lprbarr  -0.3597944  0.0594670 -6.0503 2.738e-09 ***
## lprbconv -0.2858733  0.0515213 -5.5486 4.557e-08 ***
## lprbpris -0.1827812  0.0452805 -4.0366 6.223e-05 ***
## lavgsen  -0.0044879  0.0333495 -0.1346  0.893001    
## lpolpc    0.4241142  0.0849040  4.9952 7.996e-07 ***
## d82       0.0125802  0.0160064  0.7859  0.432249    
## d83      -0.0792813  0.0195637 -4.0525 5.829e-05 ***
## d84      -0.1177281  0.0217115 -5.4224 8.957e-08 ***
## d85      -0.1119561  0.0256579 -4.3634 1.540e-05 ***
## d86      -0.0818268  0.0236273 -3.4632  0.000577 ***
## d87      -0.0404704  0.0241762 -1.6740  0.094726 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
phtest(Q7_RE, Q7_FE)
## 
##  Hausman Test
## 
## data:  Q7_formula
## chisq = 36.924, df = 11, p-value = 0.0001187
## alternative hypothesis: one model is inconsistent
Q7_FE.2 <- update(Q7_FE, .~. + lwcon+lwtuc+lwtrd+lwfir+lwser+
                    lwmfg+lwfed+lwsta+lwloc)
summary(Q7_FE.2)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = lcrmrte ~ lprbarr + lprbconv + lprbpris + lavgsen + 
##     lpolpc + d82 + d83 + d84 + d85 + d86 + d87 + lwcon + lwtuc + 
##     lwtrd + lwfir + lwser + lwmfg + lwfed + lwsta + lwloc, data = crime4, 
##     model = "within")
## 
## Balanced Panel: n = 90, T = 7, N = 630
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.5678088 -0.0640856 -0.0026452  0.0712846  0.5109494 
## 
## Coefficients:
##             Estimate  Std. Error  t-value  Pr(>|t|)    
## lprbarr  -0.35635153  0.03215909 -11.0809 < 2.2e-16 ***
## lprbconv -0.28595393  0.02105128 -13.5837 < 2.2e-16 ***
## lprbpris -0.17513548  0.03234035  -5.4154  9.36e-08 ***
## lavgsen  -0.00287394  0.02621081  -0.1096  0.912732    
## lpolpc    0.42290002  0.02639421  16.0225 < 2.2e-16 ***
## d82       0.01889147  0.02512435   0.7519  0.452440    
## d83      -0.05528597  0.03302871  -1.6739  0.094756 .  
## d84      -0.06151625  0.04108050  -1.4975  0.134881    
## d85      -0.03971153  0.05616351  -0.7071  0.479840    
## d86      -0.00011328  0.06801243  -0.0017  0.998672    
## d87       0.05370417  0.07989535   0.6722  0.501767    
## lwcon    -0.03454480  0.03916160  -0.8821  0.378125    
## lwtuc     0.04597472  0.01903400   2.4154  0.016062 *  
## lwtrd    -0.02017663  0.04060728  -0.4969  0.619489    
## lwfir    -0.00354453  0.02833305  -0.1251  0.900491    
## lwser     0.01012645  0.01919154   0.5277  0.597966    
## lwmfg    -0.30056909  0.10940682  -2.7473  0.006218 ** 
## lwfed    -0.33312264  0.17644804  -1.8879  0.059591 .  
## lwsta     0.02152087  0.11306485   0.1903  0.849116    
## lwloc     0.18102146  0.11806435   1.5332  0.125824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    17.991
## Residual Sum of Squares: 9.7606
## R-Squared:      0.45747
## Adj. R-Squared: 0.34374
## F-statistic: 21.9232 on 20 and 520 DF, p-value: < 2.22e-16
coeftest(Q7_FE.2,c*vcovHC(Q7_FE.2,type="HC1",cluster="group"))
## 
## t test of coefficients:
## 
##             Estimate  Std. Error t value  Pr(>|t|)    
## lprbarr  -0.35635153  0.06151863 -5.7926 1.201e-08 ***
## lprbconv -0.28595393  0.05076341 -5.6331 2.903e-08 ***
## lprbpris -0.17513548  0.04576161 -3.8271 0.0001454 ***
## lavgsen  -0.00287394  0.03339857 -0.0860 0.9314600    
## lpolpc    0.42290002  0.08221507  5.1438 3.820e-07 ***
## d82       0.01889147  0.02211552  0.8542 0.3933778    
## d83      -0.05528597  0.03060077 -1.8067 0.0713893 .  
## d84      -0.06151625  0.04060185 -1.5151 0.1303522    
## d85      -0.03971153  0.06033892 -0.6581 0.5107386    
## d86      -0.00011328  0.07202127 -0.0016 0.9987456    
## d87       0.05370417  0.08477270  0.6335 0.5266806    
## lwcon    -0.03454480  0.02533385 -1.3636 0.1732890    
## lwtuc     0.04597472  0.01611031  2.8537 0.0044930 ** 
## lwtrd    -0.02017663  0.03131229 -0.6444 0.5196211    
## lwfir    -0.00354453  0.01300897 -0.2725 0.7853701    
## lwser     0.01012645  0.02022431  0.5007 0.6167894    
## lwmfg    -0.30056909  0.10637182 -2.8256 0.0049000 ** 
## lwfed    -0.33312264  0.22457267 -1.4834 0.1385840    
## lwsta     0.02152087  0.10517275  0.2046 0.8379460    
## lwloc     0.18102146  0.16298603  1.1107 0.2672296    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
linearHypothesis(Q7_FE.2, c("lwcon=0","lwtuc=0","lwtrd=0","lwfir=0","lwser=0",
                            "lwmfg=0","lwfed=0","lwsta=0","lwloc=0"), 
                 test = "F", 
                 white.adjust="hc1")
## Linear hypothesis test
## 
## Hypothesis:
## lwcon = 0
## lwtuc = 0
## lwtrd = 0
## lwfir = 0
## lwser = 0
## lwmfg = 0
## lwfed = 0
## lwsta = 0
## lwloc = 0
## 
## Model 1: restricted model
## Model 2: lcrmrte ~ lprbarr + lprbconv + lprbpris + lavgsen + lpolpc + 
##     d82 + d83 + d84 + d85 + d86 + d87 + lwcon + lwtuc + lwtrd + 
##     lwfir + lwser + lwmfg + lwfed + lwsta + lwloc
## 
##   Res.Df Df      F   Pr(>F)   
## 1    529                      
## 2    520  9 2.4739 0.009047 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Q7_FD <- update(Q7_FE.2, model="fd")
summary(Q7_FD)
## Oneway (individual) effect First-Difference Model
## 
## Call:
## plm(formula = lcrmrte ~ lprbarr + lprbconv + lprbpris + lavgsen + 
##     lpolpc + d82 + d83 + d84 + d85 + d86 + d87 + lwcon + lwtuc + 
##     lwtrd + lwfir + lwser + lwmfg + lwfed + lwsta + lwloc, data = crime4, 
##     model = "fd")
## 
## Balanced Panel: n = 90, T = 7, N = 630
## Observations used in estimation: 540
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.6654347 -0.0774249  0.0056945  0.0680134  0.6815969 
## 
## Coefficients: (1 dropped because of singularities)
##               Estimate Std. Error  t-value  Pr(>|t|)    
## (Intercept)  0.0067651  0.0135424   0.4996 0.6176024    
## lprbarr     -0.3230993  0.0300195 -10.7630 < 2.2e-16 ***
## lprbconv    -0.2402885  0.0182474 -13.1683 < 2.2e-16 ***
## lprbpris    -0.1693859  0.0261700  -6.4725 2.231e-10 ***
## lavgsen     -0.0156167  0.0224126  -0.6968 0.4862508    
## lpolpc       0.3977221  0.0269870  14.7375 < 2.2e-16 ***
## d82          0.0130871  0.0160325   0.8163 0.4147123    
## d83         -0.0846911  0.0208459  -4.0627 5.599e-05 ***
## d84         -0.1090144  0.0214083  -5.0922 4.958e-07 ***
## d85         -0.0965130  0.0201721  -4.7845 2.237e-06 ***
## d86         -0.0519503  0.0153372  -3.3872 0.0007596 ***
## lwcon       -0.0442368  0.0304142  -1.4545 0.1464167    
## lwtuc        0.0253997  0.0142093   1.7875 0.0744321 .  
## lwtrd       -0.0290309  0.0307907  -0.9428 0.3461978    
## lwfir        0.0091220  0.0212318   0.4296 0.6676373    
## lwser        0.0219549  0.0144342   1.5210 0.1288588    
## lwmfg       -0.1402482  0.1019317  -1.3759 0.1694437    
## lwfed        0.0174221  0.1716065   0.1015 0.9191738    
## lwsta       -0.0517891  0.0957109  -0.5411 0.5886710    
## lwloc       -0.0305153  0.1021028  -0.2989 0.7651599    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    22.197
## Residual Sum of Squares: 12.329
## R-Squared:      0.44454
## Adj. R-Squared: 0.42425
## F-statistic: 21.9034 on 19 and 520 DF, p-value: < 2.22e-16

Questão 8

A base Grunfeld traz parte dos dados usados em um estudo clássico sobre a demanda de investimento. Os dados consistem em séries temporais de 20 observações anuais para cinco empresas e três variáveis:

  • \(I_{it}\) = Investimento bruto;
  • \(F_{it}\) = Valor de mercado da empresa no final do ano anterior;
  • \(C_{it}\) = Valor do estoque de instalações e equipamentos no final do ano anterior.

O modelo a ser estimado com esses dados é:

\[I_{it} = \beta_{1} + \beta{2}F_{it} + \beta{3}C_{it} + \epsilon_{it}.\]

Estime modelos considerando RE e FE e utilize o Teste de Hausman para determinar qual modelo é preferível para esses dados.