Verbeek 5ed. Modern Econometrics Using R
Solomon Negash
May 26, 2020
Chapter 4. Heteroskedasticity and Autocorrelation
Figure 4.1
Load libraries
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-16Auxiliary 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-16Table 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-16Table 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-10Table 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-16Hetroskedasticity-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 ' ' 1Table 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.002876Table 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-16Figure 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.3876par(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-07dwtest(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 0Figure 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))
lines(OLS9$fitted.values,col="red")
Table 4.10
EGLS (iterative Cochrane–Orcutt) results
summary(cochrane.orcutt(OLS9))## Call:
## lm(formula = cons ~ price + income + temp, data = df)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.15714805 0.28962923 0.543 0.59222
## price -0.89239640 0.81085007 -1.101 0.28157
## income 0.00320274 0.00154606 2.072 0.04878 *
## temp 0.00355839 0.00055468 6.415 1.022e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0319 on 25 degrees of freedom
## Multiple R-squared: 0.6489 , Adjusted R-squared: 0.6068
## F-statistic: 15.4 on 3 and 25 DF, p-value: < 7.015e-06
##
## Durbin-Watson statistic
## (original): 1.02117 , p-value: 3.024e-04
## (transformed): 1.54884 , p-value: 5.061e-02Table 4.11
OLS results extended specification
df <- df %>%
mutate(temp1 = lag(temp, n = 1))
summary(OLS10 <- lm(cons ~ price + income + temp + temp1, data=df))##
## Call:
## lm(formula = cons ~ price + income + temp + temp1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.049070 -0.015391 -0.006745 0.014766 0.080892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1894822 0.2323169 0.816 0.42274
## price -0.8383021 0.6880205 -1.218 0.23490
## income 0.0028673 0.0010533 2.722 0.01189 *
## temp 0.0053321 0.0006704 7.953 3.5e-08 ***
## temp1 -0.0022039 0.0007307 -3.016 0.00597 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02987 on 24 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8285, Adjusted R-squared: 0.7999
## F-statistic: 28.98 on 4 and 24 DF, p-value: 7.099e-09Figure 4.5
Actual and fitted values when true model is: \(y_t = 0.5logt+\epsilon_i.\)
set.seed(123)
t=rep(1:40)
eps = rnorm(t, mean=0, sd=0.08)
y= .5 * log(t) + eps
model <- lm(y ~ t)
plot(t,y, pch=16, cex = .4, ylim= c(0, 2), xlim = c(1,40),xlab = "Time",
sub = "Figure 4.5 Actual and fitted values when true model is y = 0.5*logt + eps", cex.sub=0.8)
abline(model, col="red")
Figure 4.6
US$/EUR and US$/GBP exchange rates, January 1979–December 2001.
df <- read_dta("Data/Forward5.dta")
df$date = as.Date(df$dateid01)
plot(df$date, df$exusbp, type = "l", col="red", xlab="", ylab="",
sub="Figure 4.6 US$/EUR and US$/GBP exchange rates, January 1979 - December 2001.", ylim=c(0.5, 2.5), cex.sub=0.8)
lines(df$date, df$exuseur, col="blue")
Figure 4.7
df$dif1 <- log(df$exusbp)-log(df$f1usbp)
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.000108df$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.52764.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-09df$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.991Hetroscedasticy 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 ' ' 1coeftest(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