Also available in Stata and Python versions
Load libraries
library(wooldridge)
library(AER)
library(stargazer)
library(plm)
library(sandwich)Demand for Air Travel
airfarep <- pdata.frame(airfare, index = c("id", "year"))
FE <- plm(log(passen) ~ log(fare) + log(dist) + log(dist)^2 + year, data=airfarep, model = "within")
RE <- plm(lpassen ~ lfare + ldist + ldistsq + year, data=airfarep, model = "random")
REIV <- plm(lpassen ~ lfare + ldist + ldistsq + year | concen + ldist + ldistsq + year, data=airfarep, model = "random")
FEIV <- plm(lpassen ~ lfare + ldist + ldistsq + year | concen + ldist + ldistsq + year, data=airfarep, model = "within")
stargazer(RE, FE, REIV, FEIV, column.labels=c("RE", "FE", "REIV", "FEIV"), title = "Table 11.1 Passenger demand model, USDR 1997-2000", no.space=TRUE, type="text")##
## Table 11.1 Passenger demand model, USDR 1997-2000
## =========================================================================
## Dependent variable:
## ------------------------------------------------------------
## lpassen log(passen) lpassen
## RE FE REIV FEIV
## (1) (2) (3) (4)
## -------------------------------------------------------------------------
## lfare -1.102*** -0.508** -0.302
## (0.022) (0.230) (0.277)
## ldist -1.971*** -1.505**
## (0.647) (0.693)
## ldistsq 0.171*** 0.118**
## (0.049) (0.055)
## log(fare) -1.155***
## (0.023)
## year1998 0.045*** 0.046*** 0.031*** 0.026***
## (0.006) (0.006) (0.009) (0.010)
## year1999 0.101*** 0.102*** 0.080*** 0.072***
## (0.006) (0.006) (0.010) (0.012)
## year2000 0.190*** 0.195*** 0.133*** 0.113***
## (0.006) (0.006) (0.023) (0.028)
## Constant 17.006*** 13.296***
## (2.133) (2.627)
## -------------------------------------------------------------------------
## Observations 4,596 4,596 4,596 4,596
## R2 0.378 0.451 0.339 0.319
## Adjusted R2 0.377 0.267 0.338 0.091
## F Statistic 2,787.704*** 706.353*** (df = 4; 3443) 231.099*** 179.428***
## =========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Robust SEs
coeftest(RE, vcovHC(RE, type="HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.0064030 2.3814920 7.1411 1.072e-12 ***
## lfare -1.1024574 0.1017023 -10.8400 < 2.2e-16 ***
## ldist -1.9707354 0.7037288 -2.8004 0.005125 **
## ldistsq 0.1709783 0.0535714 3.1916 0.001424 **
## year1998 0.0452090 0.0048053 9.4081 < 2.2e-16 ***
## year1999 0.1005163 0.0062177 16.1663 < 2.2e-16 ***
## year2000 0.1896112 0.0091983 20.6138 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(REIV, vcovHC(REIV, type="HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.296433 3.778260 3.5192 0.0004371 ***
## lfare -0.507876 0.498288 -1.0192 0.3081419
## ldist -1.504805 0.783967 -1.9195 0.0549863 .
## ldistsq 0.117601 0.067080 1.7532 0.0796428 .
## year1998 0.030736 0.013592 2.2614 0.0237800 *
## year1999 0.079655 0.020869 3.8169 0.0001369 ***
## year2000 0.132580 0.050712 2.6143 0.0089690 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(FE, vcovHC(FE, method="arellano", type="HC3"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## log(fare) -1.1550389 0.1098900 -10.5109 < 2.2e-16 ***
## year1998 0.0464889 0.0049217 9.4457 < 2.2e-16 ***
## year1999 0.1023612 0.0063251 16.1834 < 2.2e-16 ***
## year2000 0.1946548 0.0097579 19.9484 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(FEIV, vcovHC(FEIV, type="HC2"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## lfare -0.301576 0.613255 -0.4918 0.622919
## year1998 0.025715 0.016431 1.5650 0.117676
## year1999 0.072417 0.025129 2.8818 0.003979 **
## year2000 0.112791 0.062096 1.8164 0.069396 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1POLS
prisonp <- pdata.frame(prison, index = c("state", "year"))
OLS1 <- lm(gpris ~ 1 + final1 + final2 + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 + cag25_34 + cunem + cblack + cmetro + year, data=prisonp)
OLS1_robust <- vcovHC(OLS1, type="HC", cluster=c("group"))
robust.se <- sqrt(diag(OLS1_robust))
stargazer(OLS1, OLS1, se=list(NULL, robust.se),
column.labels=c("default","robust"), type= "text", single.row = T, no.space=T, align=T)##
## ==================================================================
## Dependent variable:
## -----------------------------------
## gpris
## default robust
## (1) (2)
## ------------------------------------------------------------------
## final1 -0.077*** (0.026) -0.077*** (0.016)
## final2 -0.053*** (0.018) -0.053*** (0.019)
## gpolpc -0.029 (0.044) -0.029 (0.033)
## gincpc 0.210 (0.131) 0.210 (0.191)
## cag0_14 2.617* (1.583) 2.617 (1.636)
## cag15_17 -1.609 (3.756) -1.609 (3.674)
## cag18_24 0.953 (1.731) 0.953 (1.602)
## cag25_34 -1.032 (1.763) -1.032 (1.782)
## cunem 0.162 (0.311) 0.162 (0.305)
## cblack -0.004 (0.026) -0.004 (0.025)
## cmetro -1.418* (0.786) -1.418* (0.847)
## year81 0.012 (0.014) 0.012 (0.015)
## year82 0.077*** (0.016) 0.077*** (0.017)
## year83 0.077*** (0.015) 0.077*** (0.017)
## year84 0.029 (0.018) 0.029 (0.019)
## year85 0.028* (0.016) 0.028 (0.017)
## year86 0.054*** (0.018) 0.054*** (0.021)
## year87 0.031* (0.017) 0.031* (0.018)
## year88 0.019 (0.017) 0.019 (0.017)
## year89 0.018 (0.017) 0.018 (0.017)
## year90 0.064*** (0.017) 0.064*** (0.017)
## year91 0.026 (0.017) 0.026 (0.017)
## year92 0.019 (0.018) 0.019 (0.017)
## year93 0.013 (0.019) 0.013 (0.020)
## Constant 0.027 (0.017) 0.027 (0.022)
## ------------------------------------------------------------------
## Observations 714 714
## R2 0.152 0.152
## Adjusted R2 0.123 0.123
## Residual Std. Error (df = 689) 0.062 0.062
## F Statistic (df = 24; 689) 5.153*** 5.153***
## ==================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01linearHypothesis(OLS1, vcov=vcovHC(OLS1, type="HC", cluster=c("group")), c("final1 = 0","final2 = 0"))## Linear hypothesis test
##
## Hypothesis:
## final1 = 0
## final2 = 0
##
## Model 1: restricted model
## Model 2: gpris ~ 1 + final1 + final2 + gpolpc + gincpc + cag0_14 + cag15_17 +
## cag18_24 + cag25_34 + cunem + cblack + cmetro + year
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df F Pr(>F)
## 1 691
## 2 689 2 15.432 2.778e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12SLS
summary(IV1 <- ivreg(gcriv ~ gpris + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 + cag25_34 + cunem + cblack + cmetro + year | final1 + final2 + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 + cag25_34 + cunem + cblack + cmetro + year, data=prisonp))##
## Call:
## ivreg(formula = gcriv ~ gpris + gpolpc + gincpc + cag0_14 + cag15_17 +
## cag18_24 + cag25_34 + cunem + cblack + cmetro + year | final1 +
## final2 + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 +
## cag25_34 + cunem + cblack + cmetro + year, data = prisonp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3863488 -0.0567073 -0.0002044 0.0509309 0.4133336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.014838 0.027520 0.539 0.58995
## gpris -1.031956 0.369963 -2.789 0.00543 **
## gpolpc 0.035315 0.067499 0.523 0.60101
## gincpc 0.910199 0.214327 4.247 2.47e-05 ***
## cag0_14 3.379384 2.634893 1.283 0.20008
## cag15_17 3.549945 5.766302 0.616 0.53834
## cag18_24 3.358348 2.680839 1.253 0.21073
## cag25_34 2.319993 2.706345 0.857 0.39161
## cunem 0.523696 0.478563 1.094 0.27420
## cblack -0.015848 0.040104 -0.395 0.69285
## cmetro -0.591517 1.298252 -0.456 0.64880
## year81 -0.056073 0.021735 -2.580 0.01009 *
## year82 0.028462 0.038477 0.740 0.45973
## year83 0.024703 0.037397 0.661 0.50911
## year84 0.012870 0.029334 0.439 0.66098
## year85 0.035403 0.027502 1.287 0.19844
## year86 0.092186 0.034388 2.681 0.00752 **
## year87 0.004771 0.029015 0.164 0.86944
## year88 0.053271 0.027322 1.950 0.05161 .
## year89 0.043086 0.027520 1.566 0.11790
## year90 0.144265 0.035462 4.068 5.29e-05 ***
## year91 0.061848 0.027650 2.237 0.02562 *
## year92 0.026657 0.028533 0.934 0.35050
## year93 0.022274 0.029610 0.752 0.45216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09547 on 690 degrees of freedom
## Multiple R-Squared: -0.1246, Adjusted R-squared: -0.1621
## Wald test: 6.075 on 23 and 690 DF, p-value: < 2.2e-16POLS on FD
reg1 <- lm(gcriv ~ gpris + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 + cag25_34 + cunem + cblack + cmetro + year, data=prisonp)
u1 <- resid(reg2 <- lm(gpris ~ final1 + final2 + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 + cag25_34 + cunem + cblack + cmetro + year, data=prisonp))
reg2 <- lm(gcriv ~ u1 + gpris + gpolpc + gincpc + cag0_14 + cag15_17 + cag18_24 + cag25_34 + cunem + cblack + cmetro + year, data=prisonp)
stargazer(reg1, reg2, no.space=TRUE, single.row = T, type="text")##
## ===================================================================
## Dependent variable:
## -----------------------------------------------
## gcriv
## (1) (2)
## -------------------------------------------------------------------
## u1 0.872*** (0.308)
## gpris -0.181*** (0.048) -1.032*** (0.304)
## gpolpc 0.051 (0.056) 0.035 (0.056)
## gincpc 0.738*** (0.166) 0.910*** (0.176)
## cag0_14 0.989 (2.007) 3.379 (2.168)
## cag15_17 4.984 (4.740) 3.550 (4.744)
## cag18_24 2.413 (2.191) 3.358 (2.205)
## cag25_34 2.880 (2.229) 2.320 (2.226)
## cunem 0.411 (0.394) 0.524 (0.394)
## cblack -0.015 (0.033) -0.016 (0.033)
## cmetro 0.538 (0.996) -0.592 (1.068)
## year81 -0.069*** (0.017) -0.056*** (0.018)
## year82 -0.041** (0.020) 0.028 (0.032)
## year83 -0.042** (0.020) 0.025 (0.031)
## year84 -0.014 (0.022) 0.013 (0.024)
## year85 0.009 (0.021) 0.035 (0.023)
## year86 0.044* (0.023) 0.092*** (0.028)
## year87 -0.024 (0.022) 0.005 (0.024)
## year88 0.035 (0.022) 0.053** (0.022)
## year89 0.025 (0.022) 0.043* (0.023)
## year90 0.087*** (0.021) 0.144*** (0.029)
## year91 0.039* (0.021) 0.062*** (0.023)
## year92 0.008 (0.023) 0.027 (0.023)
## year93 0.009 (0.024) 0.022 (0.024)
## Constant -0.006 (0.022) 0.015 (0.023)
## -------------------------------------------------------------------
## Observations 714 714
## R2 0.231 0.240
## Adjusted R2 0.206 0.214
## Residual Std. Error 0.079 (df = 690) 0.079 (df = 689)
## F Statistic 9.019*** (df = 23; 690) 9.065*** (df = 24; 689)
## ===================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01#coeftest(reg2, vcov=vcovHC(reg2, type="HC1", cluster="state"))Estimating a Dynamic Airfare Equation
airfarep <- pdata.frame(airfare, index = c("id", "year"))
airfarep['dlfare'] <- diff(airfarep$lfare, by="id")
airfarep['dconcen'] <- diff(airfarep$concen, by="id")
airfarep['dy00'] <- diff(airfarep$y00, by="id")
summary(PooledOLS <- plm(dlfare ~ lag(dlfare) + dconcen + dy00, data= airfarep, model="pooling"))## Pooling Model
##
## Call:
## plm(formula = dlfare ~ lag(dlfare) + dconcen + dy00, data = airfarep,
## model = "pooling")
##
## Balanced Panel: n = 1149, T = 2, N = 2298
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.8146005 -0.0480220 0.0056923 0.0573983 1.0622722
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 0.0150899 0.0035152 4.2928 1.838e-05 ***
## lag(dlfare) -0.1264673 0.0185911 -6.8026 1.307e-11 ***
## dconcen 0.0762671 0.0339566 2.2460 0.0248 *
## dy00 0.0473536 0.0048988 9.6664 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.477
## Residual Sum of Squares: 31.298
## R-Squared: 0.065095
## Adj. R-Squared: 0.063872
## F-statistic: 53.2414 on 3 and 2294 DF, p-value: < 2.22e-16Pooled IV
#Estimating separate reduced forms
airfarep['ivu1'] <- resid(lm(lag(dlfare) ~ lag(lfare, 2) , data=airfarep, subset=year=="1999"))
airfarep['ivu2'] <- resid(lm(lag(dlfare) ~ lag(lfare, 2) + lag(lfare, 3) , data=airfarep, subset=year=="2000"))
summary(PooledIV <- plm(dlfare ~ lag(dlfare) + dconcen + y00 | dconcen + ivu1 + ivu2 + dy00, data= airfarep, model="pooling")) ## Pooling Model
## Instrumental variable estimation
## (Balestra-Varadharajan-Krishnakumar's transformation)
##
## Call:
## plm(formula = dlfare ~ lag(dlfare) + dconcen + y00 | dconcen +
## ivu1 + ivu2 + dy00, data = airfarep, model = "pooling")
##
## Balanced Panel: n = 1149, T = 2, N = 2298
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.8033688 -0.0558699 0.0023602 0.0623822 0.9860717
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 0.0075120 0.0083956 0.8948 0.37092
## lag(dlfare) 0.2189715 0.3425326 0.6393 0.52265
## dconcen 0.1262795 0.0614620 2.0546 0.03992 *
## y00 0.0513845 0.0065979 7.7880 6.807e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.477
## Residual Sum of Squares: 36.008
## R-Squared: 0.001795
## Adj. R-Squared: 0.00048957
## Chisq: 99.0174 on 3 DF, p-value: < 2.22e-16[Note that the standard errors for both models are not the same as in the textbook.]
Arellano-Bond
Random Growth Model for Analyzing Enterprise Zones
ezunemp <- pdata.frame(ezunem, index = c("city", "year"))
ezunemp['dluclms'] <- diff(ezunemp$luclms, by="city")
ezunemp['dez'] <- diff(ezunemp$ez, by="city")
summary(plm(dluclms ~ dez + year, data=ezunemp, type="within"))## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = dluclms ~ dez + year, data = ezunemp, type = "within")
##
## Balanced Panel: n = 22, T = 8, N = 176
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.4240391 -0.1566907 -0.0080053 0.1472293 0.6747333
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## dez -0.191940 0.084991 -2.2584 0.02541 *
## year1982 0.778760 0.067579 11.5237 < 2.2e-16 ***
## year1983 -0.033119 0.067579 -0.4901 0.62481
## year1984 -0.014394 0.071443 -0.2015 0.84061
## year1985 0.324911 0.069323 4.6869 6.301e-06 ***
## year1986 0.292154 0.067579 4.3232 2.832e-05 ***
## year1987 0.053948 0.067579 0.7983 0.42599
## year1988 -0.017053 0.067579 -0.2523 0.80114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 20.221
## Residual Sum of Squares: 7.3344
## R-Squared: 0.63728
## Adj. R-Squared: 0.56523
## F-statistic: 32.0644 on 8 and 146 DF, p-value: < 2.22e-16#Alternatively,
summary(plm(guclms ~ cez + year, data=ezunemp, type="within"))## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = guclms ~ cez + year, data = ezunemp, type = "within")
##
## Balanced Panel: n = 22, T = 8, N = 176
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.4240391 -0.1566907 -0.0080053 0.1472293 0.6747333
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## cez -0.191940 0.084991 -2.2584 0.02541 *
## year1982 0.778760 0.067579 11.5237 < 2.2e-16 ***
## year1983 -0.033119 0.067579 -0.4901 0.62481
## year1984 -0.014394 0.071443 -0.2015 0.84061
## year1985 0.324911 0.069323 4.6869 6.301e-06 ***
## year1986 0.292154 0.067579 4.3232 2.832e-05 ***
## year1987 0.053948 0.067579 0.7983 0.42599
## year1988 -0.017053 0.067579 -0.2523 0.80114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 20.221
## Residual Sum of Squares: 7.3344
## R-Squared: 0.63728
## Adj. R-Squared: 0.56523
## F-statistic: 32.0644 on 8 and 146 DF, p-value: < 2.22e-16Testing for Coorelated Random Slopes in a Passenger Demand Equation
airfarep <- pdata.frame(airfare, index = c("id", "year"))
airfarep$concen_bar <- ave(airfarep$concen, airfarep$id)
concen2=airfarep$concen*airfarep$concen_bar
lfare2= airfarep$lfare*airfarep$concen_bar
FEIV <- plm(lpassen ~ lfare + lfare2 + year| concen + concen2 + year, data=airfarep, model = "within")
coeftest(FEIV, vcovHC(FEIV, type="HC2"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## lfare 7.928954 14.512462 0.5464 0.58486
## lfare2 -11.046300 18.588352 -0.5943 0.55238
## year1998 0.013847 0.041078 0.3371 0.73608
## year1999 0.072527 0.036570 1.9832 0.04742 *
## year2000 0.046838 0.184982 0.2532 0.80013
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1