II Econometric Analysis Using R

Also available in Stata and Python versions

Chapter 11. More Topics in Linear Unobserved Effects Models

Example 11.1

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

HOME | Back to top

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

HOME | Back to top

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
## Adj. R-Squared: 0.063872
## 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
##    (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

HOME | Back to top

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

HOME | Back to top

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

HOME | Back to top