II Econometric Analysis Using R

Also available in Stata and Python versions

### Chapter 11. More Topics in Linear Unobserved Effects Models

#### Example 11.1

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.01

Robust 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 ' ' 1
coeftest(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 ' ' 1
coeftest(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 ' ' 1
coeftest(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 ' ' 1

#### Example 11.2

POLS

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
## 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.01
linearHypothesis(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 ' ' 1

2SLS

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-16

POLS 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
## 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"))

#### Example 11.3

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")
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-16 #### Example 11.5 Testing 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