Introductory Econometrics Using R

Chapter 12. Serial Correlation and Heteroskedasticity in Time Series Regressions

Also covered using Python and Stata

library(wooldridge)
library(lmtest)
library(stargazer)
library(car)
library(dynlm)
library(prais)
library(sandwich)
options(width=120)

Example 12.1. Testing for AR(1) Serial Correlation in the Phillips Curve

ts_phillips <- ts(phillips)
uhat1 <- resid(dynlm(inf ~ unem + 1, data=ts_phillips))
phill_inf1 <- dynlm(uhat1 ~ L(uhat1,1))

ts_phillips2 <-ts(subset(phillips, phillips$year<1997))
uhat2 <- resid(dynlm(inf ~ unem + 1, data=ts_phillips2))
phill_inf2 <- dynlm(uhat2 ~ L(uhat2, 1))

uhat3<- resid(dynlm(cinf ~ unem + 1, data=ts_phillips2))
phill_cinf <- dynlm(uhat3 ~ L(uhat3, 1))

stargazer(phill_inf1, phill_inf2, phill_cinf, keep.stat=c("n","rsq", "adj.rsq"), no.space=TRUE, type="text")

## 
## ==========================================
##                   Dependent variable:     
##              -----------------------------
##                uhat1      uhat2    uhat3  
##                 (1)        (2)      (3)   
## ------------------------------------------
## L(uhat1, 1)   0.572***                    
##               (0.108)                     
## L(uhat2, 1)             0.573***          
##                          (0.116)          
## L(uhat3, 1)                        -0.036 
##                                   (0.124) 
## Constant       -0.112    -0.113    0.194  
##               (0.318)    (0.359)  (0.300) 
## ------------------------------------------
## Observations     55        48        47   
## R2             0.345      0.346    0.002  
## Adjusted R2    0.333      0.332    -0.020 
## ==========================================
## Note:          *p<0.1; **p<0.05; ***p<0.01

Example 12.2. Testing for AR(1) Serial Correlation in the Minimum Wage Equation

tsprminwge <- ts(prminwge)
minwg_reg <- dynlm(lprepop ~ lmincov + lprgnp + lusgnp + t, data=tsprminwge)
uhat <- resid(minwg_reg)
AR1c <- dynlm(uhat ~  lmincov + lprgnp + lusgnp + t + L(uhat, 1) , data=tsprminwge)
AR1s <- dynlm(uhat ~ L(uhat, 1))
stargazer(AR1c, AR1s, column.labels=c("AR1c", "AR1s"), no.space=TRUE, type="text")

## 
## ===========================================================
##                               Dependent variable:          
##                     ---------------------------------------
##                                      uhat                  
##                            AR1c                AR1s        
##                            (1)                 (2)         
## -----------------------------------------------------------
## lmincov                   0.038                            
##                          (0.035)                           
## lprgnp                    -0.078                           
##                          (0.071)                           
## lusgnp                    0.204                            
##                          (0.195)                           
## t                         -0.003                           
##                          (0.004)                           
## L(uhat, 1)               0.481***            0.417**       
##                          (0.166)             (0.159)       
## Constant                  -0.851              -0.001       
##                          (1.093)             (0.004)       
## -----------------------------------------------------------
## Observations                37                  37         
## R2                        0.242               0.165        
## Adjusted R2               0.120               0.141        
## Residual Std. Error  0.028 (df = 31)     0.027 (df = 35)   
## F Statistic         1.983 (df = 5; 31) 6.895** (df = 1; 35)
## ===========================================================
## Note:                           *p<0.1; **p<0.05; ***p<0.01

Example 12.3. Testing for AR(3) Serial Correlation

tsbarium<-ts(barium)
barium_reg <- dynlm(lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6, data=tsbarium)
u <- resid(barium_reg)
AR3 <- dynlm(u ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6 + L(u, 1) + L(u, 2) + L(u, 3) + 1, data = tsbarium)
summary(AR3)

## 
## Time series regression with "ts" data:
## Start = 4, End = 131
## 
## Call:
## dynlm(formula = u ~ lchempi + lgas + lrtwex + befile6 + affile6 + 
##     afdec6 + L(u, 1) + L(u, 2) + L(u, 3) + 1, data = tsbarium)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.89072 -0.32250  0.05873  0.36376  1.19650 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -14.36915   20.65567  -0.696   0.4880  
## lchempi      -0.14316    0.47203  -0.303   0.7622  
## lgas          0.62331    0.88597   0.704   0.4831  
## lrtwex        0.17867    0.39103   0.457   0.6486  
## befile6      -0.08592    0.25101  -0.342   0.7327  
## affile6      -0.12212    0.25470  -0.479   0.6325  
## afdec6       -0.06683    0.27437  -0.244   0.8080  
## L(u, 1)       0.22149    0.09166   2.417   0.0172 *
## L(u, 2)       0.13404    0.09216   1.454   0.1485  
## L(u, 3)       0.12554    0.09112   1.378   0.1709  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5704 on 118 degrees of freedom
## Multiple R-squared:  0.1159, Adjusted R-squared:  0.04848 
## F-statistic: 1.719 on 9 and 118 DF,  p-value: 0.09202

linearHypothesis(AR3, c("L(u, 1)=0 ", "L(u, 2)=0" , "L(u, 3)=0"))

## Linear hypothesis test
## 
## Hypothesis:
## L(u,0
## L(u, 2) = 0
## L(u, 3) = 0
## 
## Model 1: restricted model
## Model 2: u ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6 + L(u, 
##     1) + L(u, 2) + L(u, 3) + 1
## 
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)   
## 1    121 43.394                               
## 2    118 38.394  3    5.0005 5.1229 0.00229 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example 12.4. Prais-Winsten Estimation in the Event Study

ols <- dynlm(lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6, data=tsbarium)
stargazer(ols, column.labels=c("OLS"), no.space=TRUE, single.row = TRUE, type="text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               lchnimp          
##                                 OLS            
## -----------------------------------------------
## lchempi                  3.117*** (0.479)      
## lgas                       0.196 (0.907)       
## lrtwex                    0.983** (0.400)      
## befile6                    0.060 (0.261)       
## affile6                   -0.032 (0.264)       
## afdec6                    -0.565* (0.286)      
## Constant                 -17.803 (21.045)      
## -----------------------------------------------
## Observations                    131            
## R2                             0.305           
## Adjusted R2                    0.271           
## Residual Std. Error      0.597 (df = 124)      
## F Statistic           9.064*** (df = 6; 124)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

pw <- prais_winsten(lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6, data=tsbarium)

## Iteration 0: rho = 0
## Iteration 1: rho = 0.2708
## Iteration 2: rho = 0.291
## Iteration 3: rho = 0.293
## Iteration 4: rho = 0.2932
## Iteration 5: rho = 0.2932
## Iteration 6: rho = 0.2932
## Iteration 7: rho = 0.2932

summary(pw)

## 
## Call:
## prais_winsten(formula = lchnimp ~ lchempi + lgas + lrtwex + befile6 + 
##     affile6 + afdec6, data = tsbarium)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.99386 -0.32219  0.03747  0.40226  1.50281 
## 
## AR(1) coefficient rho after 7 Iterations: 0.2932
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -37.07770   22.77830  -1.628   0.1061    
## lchempi       2.94095    0.63284   4.647 8.46e-06 ***
## lgas          1.04638    0.97734   1.071   0.2864    
## lrtwex        1.13279    0.50666   2.236   0.0272 *  
## befile6      -0.01648    0.31938  -0.052   0.9589    
## affile6      -0.03316    0.32181  -0.103   0.9181    
## afdec6       -0.57681    0.34199  -1.687   0.0942 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5733 on 124 degrees of freedom
## Multiple R-squared:  0.2021, Adjusted R-squared:  0.1635 
## F-statistic: 5.235 on 6 and 124 DF,  p-value: 7.764e-05
## 
## Durbin-Watson statistic (original): 1.458 
## Durbin-Watson statistic (transformed): 2.087

Example 12.5. Static Phillips Curve

phillips2 <- ts(subset(phillips, phillips$year<1997))
ols <- lm(inf ~ unem, data=phillips2)
stargazer(ols, column.labels=c("OLS"), no.space=TRUE, single.row = TRUE, type="text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 inf            
##                                 OLS            
## -----------------------------------------------
## unem                       0.468 (0.289)       
## Constant                   1.424 (1.719)       
## -----------------------------------------------
## Observations                    49             
## R2                             0.053           
## Adjusted R2                    0.033           
## Residual Std. Error       3.131 (df = 47)      
## F Statistic             2.616 (df = 1; 47)     
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

pw <- prais_winsten(inf ~ unem, data=phillips2)

## Iteration 0: rho = 0
## Iteration 1: rho = 0.5727
## Iteration 2: rho = 0.7307
## Iteration 3: rho = 0.7719
## Iteration 4: rho = 0.7792
## Iteration 5: rho = 0.7803
## Iteration 6: rho = 0.7805
## Iteration 7: rho = 0.7805
## Iteration 8: rho = 0.7805
## Iteration 9: rho = 0.7805

summary(pw)

## 
## Call:
## prais_winsten(formula = inf ~ unem, data = phillips2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5470 -2.6764 -0.3579  1.4108 10.2853 
## 
## AR(1) coefficient rho after 9 Iterations: 0.7805
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.2959     2.2314   3.718 0.000535 ***
## unem         -0.7157     0.3135  -2.283 0.026988 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.267 on 47 degrees of freedom
## Multiple R-squared:  0.1358, Adjusted R-squared:  0.1174 
## F-statistic: 7.386 on 1 and 47 DF,  p-value: 0.009174
## 
## Durbin-Watson statistic (original): 0.8027 
## Durbin-Watson statistic (transformed):  1.91

Example 12.6. Differencing the Interest Rate Equation

tsintdef <- ts(intdef)
reg1 <- dynlm(i3 ~ inf + def, data=tsintdef)
u <- resid(reg1)
AR1u <- dynlm(u ~ L(u, 1), data=tsintdef)

reg2 <-dynlm(ci3 ~ cinf + cdef, data=tsintdef)
e <- resid(reg2)
AR1e <- dynlm(e ~ L(e, 1), data=tsintdef)

stargazer(reg1, reg2, AR1u, AR1e, column.labels=c("reg1", "reg2", "AR1u", "AR1e"), no.space=TRUE, keep.stat=c("n","rsq", "adj.rsq"), type="text")

## 
## ==============================================
##                     Dependent variable:       
##              ---------------------------------
##                 i3      ci3      u        e   
##                reg1    reg2     AR1u    AR1e  
##                (1)      (2)     (3)      (4)  
## ----------------------------------------------
## inf          0.606***                         
##              (0.082)                          
## def          0.513***                         
##              (0.118)                          
## cinf                   0.149                  
##                       (0.092)                 
## cdef                  -0.181                  
##                       (0.148)                 
## L(u, 1)                       0.623***        
##                               (0.110)         
## L(e, 1)                                 0.072 
##                                        (0.134)
## Constant     1.733***  0.042   0.015   -0.041 
##              (0.432)  (0.171) (0.190)  (0.166)
## ----------------------------------------------
## Observations    56      55       55      54   
## R2            0.602    0.176   0.377    0.005 
## Adjusted R2   0.587    0.145   0.366   -0.014 
## ==============================================
## Note:              *p<0.1; **p<0.05; ***p<0.01

Example 12.7. The Puerto Rican Minimum Wage

OLS <- lm(lprepop ~ lmincov + lprgnp + lusgnp + t, data=prminwge)
Newey <- coeftest(OLS, vcov.=NeweyWest(OLS, lag=2, prewhite=FALSE, adjust=TRUE, verbose=TRUE))

## 
## Lag truncation parameter chosen: 2

stargazer(OLS, Newey, column.labels=c("OLS", "Newey"), no.space=TRUE, keep.stat=c("n","rsq", "adj.rsq"), type="text")

## 
## =========================================
##                  Dependent variable:     
##              ----------------------------
##                 lprepop                  
##                   OLS        coefficient 
##                                 test     
##                   OLS           Newey    
##                   (1)            (2)     
## -----------------------------------------
## lmincov        -0.212***      -0.212***  
##                 (0.040)        (0.046)   
## lprgnp          0.285***      0.285***   
##                 (0.080)        (0.100)   
## lusgnp          0.486**        0.486*    
##                 (0.222)        (0.279)   
## t              -0.027***      -0.027***  
##                 (0.005)        (0.006)   
## Constant       -6.663***      -6.663***  
##                 (1.258)        (1.536)   
## -----------------------------------------
## Observations       38                    
## R2               0.889                   
## Adjusted R2      0.876                   
## =========================================
## Note:         *p<0.1; **p<0.05; ***p<0.01

summary(prais_winsten(lprepop ~  lmincov + lprgnp + lusgnp + t , data=tsprminwge))

## Iteration 0: rho = 0
## Iteration 1: rho = 0.4197
## Iteration 2: rho = 0.5325
## Iteration 3: rho = 0.5796
## Iteration 4: rho = 0.5999
## Iteration 5: rho = 0.6086
## Iteration 6: rho = 0.6123
## Iteration 7: rho = 0.6139
## Iteration 8: rho = 0.6146
## Iteration 9: rho = 0.6149
## Iteration 10: rho = 0.615
## Iteration 11: rho = 0.6151
## Iteration 12: rho = 0.6151
## Iteration 13: rho = 0.6151
## Iteration 14: rho = 0.6151
## Iteration 15: rho = 0.6151
## Iteration 16: rho = 0.6151

## 
## Call:
## prais_winsten(formula = lprepop ~ lmincov + lprgnp + lusgnp + 
##     t, data = tsprminwge)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.079574 -0.028121 -0.004816  0.005918  0.075978 
## 
## AR(1) coefficient rho after 16 Iterations: 0.6151
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -4.652853   1.376470  -3.380  0.00188 **
## lmincov     -0.147711   0.045842  -3.222  0.00286 **
## lprgnp       0.251383   0.116462   2.158  0.03826 * 
## lusgnp       0.255711   0.231750   1.103  0.27784   
## t           -0.020502   0.005856  -3.501  0.00135 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02852 on 33 degrees of freedom
## Multiple R-squared:  0.7509, Adjusted R-squared:  0.7207 
## F-statistic: 24.87 on 4 and 33 DF,  p-value: 1.47e-09
## 
## Durbin-Watson statistic (original): 1.014 
## Durbin-Watson statistic (transformed): 1.736

Example 12.8. Heteroskedasticity and the Efficient Markets Hypothesis

tsnyse <-ts(nyse)
u2 <- resid(dynlm(return ~ return_1, data=tsnyse))**2
summary(dynlm(u2 ~ return_1 + 1, data = tsnyse))

## 
## Time series regression with "ts" data:
## Start = 3, End = 691
## 
## Call:
## dynlm(formula = u2 ~ return_1 + 1, data = tsnyse)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -9.689  -3.929  -2.021   0.960 223.730 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.6565     0.4277  10.888  < 2e-16 ***
## return_1     -1.1041     0.2014  -5.482  5.9e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.18 on 687 degrees of freedom
## Multiple R-squared:  0.04191,    Adjusted R-squared:  0.04052 
## F-statistic: 30.05 on 1 and 687 DF,  p-value: 5.905e-08

Example 12.9. ARCH in Stock Returns

tsnyse <-ts(nyse)
u <- resid(dynlm(return ~ return_1, data=tsnyse))
u2 <- u*u
u2_reg <- dynlm(u2 ~ L(u2, 1) + 1, data = tsnyse)
u_reg <- dynlm(u ~ L(u, 1) + 1)
stargazer(u2_reg, u_reg,  no.space=TRUE, type="text")

## 
## ===========================================================
##                                    Dependent variable:     
##                                ----------------------------
##                                      u2             u      
##                                      (1)           (2)     
## -----------------------------------------------------------
## L(u2, 1)                          0.337***                 
##                                    (0.036)                 
## L(u, 1)                                           0.001    
##                                                  (0.038)   
## Constant                          2.947***        -0.001   
##                                    (0.440)       (0.081)   
## -----------------------------------------------------------
## Observations                         688           688     
## R2                                  0.114        0.00000   
## Adjusted R2                         0.112         -0.001   
## Residual Std. Error (df = 686)     10.759         2.112    
## F Statistic (df = 1; 686)         87.923***       0.001    
## ===========================================================
## Note:                           *p<0.1; **p<0.05; ***p<0.01