Verbeek5. Modern Econometrics Using R - Chap 4

### Chapter 4. Heteroskedasticity and Autocorrelation

#### Figure 4.1

library(haven)
library(AER)
library(orcutt)
library(dplyr)

An Engel curve with heteroskedasticity, hypothetical data

set.seed(123)
x=rep(50:500, 1)
sigma2 = x^1.35
eps = rnorm(x,mean=0,sd=sqrt(sigma2))
y= 30 + 0.75*x + eps
model <- lm(y  ~  x)
plot(x,y, pch=16, cex = .4, ylim = c(0,500), xlim = c(0,600),
main = "Figure 4.1 An Engel curve with heteroskedasticity",  cex.main=0.8)
abline(model, col="red")

#### Table 4.1

OLS results linear model

df <- read_dta("Data/labour2.dta")

summary(OLS1 <- lm(labor ~ wage + output + capital, data=df))
##
## Call:
## lm(formula = labor ~ wage + output + capital, data = df)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -1267.04   -54.11   -14.02    37.20  1560.48
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 287.7186    19.6418   14.65   <2e-16 ***
## wage         -6.7419     0.5014  -13.45   <2e-16 ***
## output       15.4005     0.3556   43.30   <2e-16 ***
## capital      -4.5905     0.2690  -17.07   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 156.3 on 565 degrees of freedom
## Multiple R-squared:  0.9352, Adjusted R-squared:  0.9348
## F-statistic:  2716 on 3 and 565 DF,  p-value: < 2.2e-16

Auxiliary regression Breuschâ€“Pagan test

u2 <- resid(OLS1)^2
summary(OLS2 <- lm(u2 ~ wage + output + capital, data=df))
##
## Call:
## lm(formula = u2 ~ wage + output + capital, data = df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -500023  -12448    2722   13354 1193685
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22719.5    11838.9  -1.919   0.0555 .
## wage           228.9      302.2   0.757   0.4492
## output        5362.2      214.4  25.016   <2e-16 ***
## capital      -3543.5      162.1 -21.858   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 94180 on 565 degrees of freedom
## Multiple R-squared:  0.5818, Adjusted R-squared:  0.5796
## F-statistic:   262 on 3 and 565 DF,  p-value: < 2.2e-16

#### Table 4.3

OLS results loglinear model

summary(OLS3 <- lm(log(labor) ~ log(wage) + log(output) + log(capital), data=df))
##
## Call:
## lm(formula = log(labor) ~ log(wage) + log(output) + log(capital),
##     data = df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -3.9744 -0.1641  0.1079  0.2595  1.9466
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)   6.177290   0.246211  25.089   <2e-16 ***
## log(wage)    -0.927764   0.071405 -12.993   <2e-16 ***
## log(output)   0.990047   0.026410  37.487   <2e-16 ***
## log(capital) -0.003697   0.018770  -0.197    0.844
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4653 on 565 degrees of freedom
## Multiple R-squared:  0.843,  Adjusted R-squared:  0.8421
## F-statistic:  1011 on 3 and 565 DF,  p-value: < 2.2e-16

#### Table 4.4

Auxiliary regression White test

u2 <- resid(OLS3)**2
summary(OLS4 <- lm(u2 ~ log(wage) + log(output) + log(capital) +
I(log(wage)^2) + I(log(output)^2) + I(log(capital)^2) +
I(log(wage)*log(output)) + I(log(wage)*log(capital)) + I(log(output)*log(capital)), data=df))
##
## Call:
## lm(formula = u2 ~ log(wage) + log(output) + log(capital) + I(log(wage)^2) +
##     I(log(output)^2) + I(log(capital)^2) + I(log(wage) * log(output)) +
##     I(log(wage) * log(capital)) + I(log(output) * log(capital)),
##     data = df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -2.2664 -0.1650 -0.0724  0.0212 15.2247
##
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)
## (Intercept)                    2.54460    3.00278   0.847 0.397126
## log(wage)                     -1.29900    1.75274  -0.741 0.458929
## log(output)                   -0.90372    0.55985  -1.614 0.107045
## log(capital)                   1.14205    0.37582   3.039 0.002486 **
## I(log(wage)^2)                 0.19274    0.25895   0.744 0.457003
## I(log(output)^2)               0.13820    0.03565   3.877 0.000118 ***
## I(log(capital)^2)              0.08954    0.01399   6.401 3.27e-10 ***
## I(log(wage) * log(output))     0.13804    0.16256   0.849 0.396168
## I(log(wage) * log(capital))   -0.25178    0.10497  -2.399 0.016782 *
## I(log(output) * log(capital)) -0.19160    0.03687  -5.197 2.84e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8507 on 559 degrees of freedom
## Multiple R-squared:  0.1029, Adjusted R-squared:  0.08845
## F-statistic: 7.124 on 9 and 559 DF,  p-value: 8.334e-10

#### Table 4.5

OLS results loglinear model with White standard errors

summary(OLS5 <- lm(log(labor) ~ log(wage) + log(output) + log(capital), data=df))
##
## Call:
## lm(formula = log(labor) ~ log(wage) + log(output) + log(capital),
##     data = df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -3.9744 -0.1641  0.1079  0.2595  1.9466
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)   6.177290   0.246211  25.089   <2e-16 ***
## log(wage)    -0.927764   0.071405 -12.993   <2e-16 ***
## log(output)   0.990047   0.026410  37.487   <2e-16 ***
## log(capital) -0.003697   0.018770  -0.197    0.844
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4653 on 565 degrees of freedom
## Multiple R-squared:  0.843,  Adjusted R-squared:  0.8421
## F-statistic:  1011 on 3 and 565 DF,  p-value: < 2.2e-16

Hetroskedasticity-robust standard errors

coeftest(OLS5, vcovHC(OLS5, type="HC1"))
##
## t test of coefficients:
##
##                Estimate Std. Error  t value Pr(>|t|)
## (Intercept)   6.1772896  0.2938869  21.0193   <2e-16 ***
## log(wage)    -0.9277642  0.0866604 -10.7057   <2e-16 ***
## log(output)   0.9900474  0.0467902  21.1593   <2e-16 ***
## log(capital) -0.0036975  0.0378770  -0.0976   0.9223
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#### Table 4.6

Auxiliary regression multiplicative heteroskedasticity

summary(OLS6 <- lm(log(u2) ~ log(wage) + log(output) + log(capital), data=df))
##
## Call:
## lm(formula = log(u2) ~ log(wage) + log(output) + log(capital),
##     data = df)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -11.7447  -0.7645   0.3281   1.1430   6.7871
##
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -3.25382    1.18545  -2.745 0.006247 **
## log(wage)    -0.06105    0.34380  -0.178 0.859112
## log(output)   0.26695    0.12716   2.099 0.036232 *
## log(capital) -0.33069    0.09037  -3.659 0.000277 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.24 on 565 degrees of freedom
## Multiple R-squared:  0.02449,    Adjusted R-squared:  0.01931
## F-statistic: 4.728 on 3 and 565 DF,  p-value: 0.002876

#### Table 4.7

EGLS results loglinear model

wt <- 1/exp(predict(OLS6))
summary(OLS7 <- lm(log(labor) ~ log(wage) + log(output) + log(capital), weights=wt, data=df))
##
## Call:
## lm(formula = log(labor) ~ log(wage) + log(output) + log(capital),
##     data = df, weights = wt)
##
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max
## -29.1086  -0.7875   0.4852   1.2394  10.4219
##
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)   5.89536    0.24764  23.806  < 2e-16 ***
## log(wage)    -0.85558    0.07188 -11.903  < 2e-16 ***
## log(output)   1.03461    0.02731  37.890  < 2e-16 ***
## log(capital) -0.05686    0.02158  -2.636  0.00863 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.509 on 565 degrees of freedom
## Multiple R-squared:  0.8509, Adjusted R-squared:  0.8501
## F-statistic:  1074 on 3 and 565 DF,  p-value: < 2.2e-16

#### Figure 4.2

df <- read_dta("Data/icecream.dta")
summary(OLS8 <- lm(cons ~ income + price, data=df))
##
## Call:
## lm(formula = cons ~ income + price, data = df)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -0.10951 -0.03865 -0.01045  0.03665  0.15635
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.9002396  0.4550343   1.978   0.0582 .
## income       0.0002135  0.0019687   0.108   0.9144
## price       -2.0300367  1.4738935  -1.377   0.1797
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06583 on 27 degrees of freedom
## Multiple R-squared:  0.0678, Adjusted R-squared:  -0.001257
## F-statistic: 0.9818 on 2 and 27 DF,  p-value: 0.3876
par(mar=c(4, 4.1, 1.1, 2))
plot(df$cons, type="p", col="blue", las=1, pch=19, xlab="Time", ylab="Consumption", main="Figure 4.2 Actual and fitted consumption of ice cream, March 1951â€“July 1953.", cex.main=0.8) lines(OLS8$fitted.values, col="red")

#### Figure 4.3

Ice cream consumption, price and temperature/100

plot(df$temp/100, type="o", las=1, xlab="Time", ylab="", pch=19, xaxs="i", main="Figure 4.3 Ice cream consumption, price and temperature/100.", cex.main=0.8) lines(df$cons, type="o", pch=17,  col="red")
lines(df$price,type="o", pch=18, col="blue") legend(7.5, .7, legend=c("Consumption", "Temp/100", "Price"), col=c("red","black" ,"blue"), lty=1:2, cex=0.75, box.lty=0) #### Table 4.9 OLS results summary(OLS9 <- lm(cons ~ price + income + temp, data=df)) ## ## Call: ## lm(formula = cons ~ price + income + temp, data = df) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.065302 -0.011873 0.002737 0.015953 0.078986 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.1973149 0.2702161 0.730 0.47179 ## price -1.0444131 0.8343570 -1.252 0.22180 ## income 0.0033078 0.0011714 2.824 0.00899 ** ## temp 0.0034584 0.0004455 7.762 3.1e-08 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.03683 on 26 degrees of freedom ## Multiple R-squared: 0.719, Adjusted R-squared: 0.6866 ## F-statistic: 22.17 on 3 and 26 DF, p-value: 2.451e-07 dwtest(cons ~ price + income + temp, data=df) ## ## Durbin-Watson test ## ## data: cons ~ price + income + temp ## DW = 1.0212, p-value = 0.0003024 ## alternative hypothesis: true autocorrelation is greater than 0 #### Figure 4.4 Actual and fitted values (connected) of ice cream consumption plot(df$cons, type="p",las=1, col="blue", ylab="Consumption",
main="Figure 4.4 Actual and fitted values (connected) of ice cream consumption",
cex.main=0.8, xlab="Time", pch=19, ylim=c(0.2, 0.6), xlim = c(0, 30))
df$dif2 <- log(df$exuseur)-log(df$f1useur) plot(df$date, df$dif1, sub="Figure 4.7 Forward discount, US$/EUR and US$/GBP, January 1979â€“December 2001", type="l", cex.sub=0.8, las=1, xlab="", ylab="", ylim=c(-0.02,0.01), col="red") lines(df$date,df$dif2,col="blue") #### Table 4.12 OLS results USD/PoundSterling df$exusbp2 <- as.vector(log(df$exusbp)) df$f1usbp2 <- as.vector(log(df$f1usbp)) st <- embed(df$exusbp2, 2)
ft<- embed(df$f1usbp2, 2) summary(lm(I(st[,1]-ft[,(2)]) ~ I(st[,(2)]-ft[,(2)]))) ## ## Call: ## lm(formula = I(st[, 1] - ft[, (2)]) ~ I(st[, (2)] - ft[, (2)])) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.14766 -0.01909 0.00073 0.02082 0.12527 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.005112 0.002365 -2.162 0.031514 * ## I(st[, (2)] - ft[, (2)]) 3.212170 0.817474 3.929 0.000108 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.03154 on 273 degrees of freedom ## Multiple R-squared: 0.05353, Adjusted R-squared: 0.05006 ## F-statistic: 15.44 on 1 and 273 DF, p-value: 0.000108 df$exuseur2 <- as.vector(log(df$exuseur)) df$f1useur2 <- as.vector(log(df$f1useur)) st <- embed(df$exuseur2, 2)
ft<- embed(df$f1useur2, 2) summary(lm(I(st[,1]-ft[,(2)]) ~ I(st[,(2)]-ft[,(2)]))) ## ## Call: ## lm(formula = I(st[, 1] - ft[, (2)]) ~ I(st[, (2)] - ft[, (2)])) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.103024 -0.021487 -0.000015 0.020975 0.088699 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.002280 0.003149 -0.724 0.470 ## I(st[, (2)] - ft[, (2)]) 0.484791 0.766435 0.633 0.528 ## ## Residual standard error: 0.03368 on 273 degrees of freedom ## Multiple R-squared: 0.001463, Adjusted R-squared: -0.002194 ## F-statistic: 0.4001 on 1 and 273 DF, p-value: 0.5276 #### 4.11.3 Tests for Risk Premia Using Overlapping Samples df$f3usbp2 <- as.vector(log(df$f3usbp)) st <- embed(df$exusbp2, 4)
ft<- embed(df$f3usbp2, 4) summary(USD_BP_3m <- lm(I(st[,1]-ft[,(4)]) ~ I(st[,(4)]-ft[,(4)]))) ## ## Call: ## lm(formula = I(st[, 1] - ft[, (4)]) ~ I(st[, (4)] - ft[, (4)])) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.285511 -0.025561 0.001782 0.029698 0.176615 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -0.013566 0.004216 -3.218 0.00145 ** ## I(st[, (4)] - ft[, (4)]) 3.135215 0.529277 5.924 9.53e-09 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.05647 on 271 degrees of freedom ## Multiple R-squared: 0.1146, Adjusted R-squared: 0.1114 ## F-statistic: 35.09 on 1 and 271 DF, p-value: 9.525e-09 df$f3useur2 <- as.vector(log(df$f3useur)) st <- embed(df$exuseur2, 4)
ft<- embed(df\$f3useur2, 4)
summary(USD_Euro_3m <- lm(I(st[,1] - ft[,(4)]) ~ I(st[,(4)] - ft[,(4)])))
##
## Call:
## lm(formula = I(st[, 1] - ft[, (4)]) ~ I(st[, (4)] - ft[, (4)]))
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -0.15097 -0.04598 -0.00358  0.04268  0.15541
##
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)              -0.010506   0.005983  -1.756   0.0802 .
## I(st[, (4)] - ft[, (4)])  0.006050   0.534784   0.011   0.9910
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06059 on 271 degrees of freedom
## Multiple R-squared:  4.722e-07,  Adjusted R-squared:  -0.00369
## F-statistic: 0.000128 on 1 and 271 DF,  p-value: 0.991

Hetroscedasticy robust errors

coeftest(USD_BP_3m, vcovHC(USD_BP_3m, type="HC1"))
##
## t test of coefficients:
##
##                            Estimate Std. Error t value  Pr(>|t|)
## (Intercept)              -0.0135664  0.0038136 -3.5574  0.000442 ***
## I(st[, (4)] - ft[, (4)])  3.1352149  0.7279863  4.3067 2.319e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(USD_Euro_3m, vcovHC(USD_Euro_3m, type="HC1"))
##
## t test of coefficients:
##
##                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)              -0.0105060  0.0058511 -1.7956  0.07368 .
## I(st[, (4)] - ft[, (4)])  0.0060495  0.5404237  0.0112  0.99108
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1