# 11 More Topics in Linear Unobserved Effects Models

Also available in Stata and Python versions

## 11.1 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``````

## 11.2 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"))``

## 11.3 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")
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
## F-statistic: 53.2414 on 3 and 2294 DF, p-value: < 2.22e-16``````

Pooled 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
##
## 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
## 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

## 11.4 Example 11.4

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
## 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
## F-statistic: 32.0644 on 8 and 146 DF, p-value: < 2.22e-16``````

## 11.5 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``````