env <- new.env()
load("../data.rda", envir = env)
Since this data is entirely constructed I will define a function that will build the true regression output and the error terms given a set of distrubances (u). It will print out the error terms against time (Xt) while accepting additional plotting parameters.
Note that the interested reader could try other disturbances from the standard normal distribution by using rnorm, which defaults to mean = 0 and sd = 1. Depending on the size required, the function will build a dataframe of that size.
error <- function(u, estart = 3.0, plot = FALSE, ...)
{
e <- vector(mode = "numeric", length = length(u))
y <- vector(mode = "numeric", length = length(u))
Xt <- 0:10
for(i in seq(u))
if (i != 1) e[i] <- e[i-1] + u[i] else e[i] <- estart
y <- 2 + 0.5 * Xt + e # Define the response from systematic Xt and error term u
if (plot)
plot(e, pch = 19, ...)
return(data.frame(u = u, e = e, y = y, x = 0:10))
}
u <- c(0, 0.5, -0.7, 0.3, 0, -2.3, -1.9, 0.2, -0.3, 0.2, -0.1)
df <- error(u, xlab = "Time", plot = TRUE, ylim = c(-3, 4))
abline(0, 0)
print(df)
## u e y x
## 1 0.0 3.0 5.0 0
## 2 0.5 3.5 6.0 1
## 3 -0.7 2.8 5.8 2
## 4 0.3 3.1 6.6 3
## 5 0.0 3.1 7.1 4
## 6 -2.3 0.8 5.3 5
## 7 -1.9 -1.1 3.9 6
## 8 0.2 -0.9 4.6 7
## 9 -0.3 -1.2 4.8 8
## 10 0.2 -1.0 5.5 9
## 11 -0.1 -1.1 5.9 10
Since © uses “different disturbances” the results will differ each time.
plot(y ~ x, df, ylim = c(0, 10), main = "(a) True Regression Line e0 = 3", pch = 19)
abline(2, 0.5)
plot(y ~ x, df, ylim = c(0, 10), main = "(b) Fitted Regression Line e0 = 3", pch = 19)
abline(lm(y ~ x, df))
# Create new data set
df <- error(rnorm(11), estart = -0.2)
plot(y ~ x, df, ylim = c(0, 10), main = "(c) Fitted Regression Line e0 = -0.2", pch = 19)
abline(lm(y ~ x, df))
Since R has extensive facilities for dealing with time series data, I will convert the data to a time series (ts) object with the appropriate time interval. For more see the Vito Ricci's reference card
It will be noted now that throughout the rest of this chapter there are two references to many of the terms, either by t or t-1. Since R is inherently vectorized and each of these vectors include n-1 terms from the n-sized vector, we can unambiguously define t = 2:20 and t-1 = 1:19 by the vectorized arithmetic. This will make a lot of the notation clearer.
df <- get("CH12TA02", envir = env)
names(df) <- c("Y", "X")
fit <- lm(Y ~ X, df)
t <- seq(2, 20)
plot(resid(fit), xlab = "Time", ylim = c(-0.5, .5), pch = 19)
abline(0,0)
R contains a Durban-Watson test dwtest (lmtest). See below. It is possible to cbind the table and the data object, but it would be superfluous. I then calculate D which can be checked against Table B.7 for n = 20, alpha = 0.1, and p = 2. The test conclusion is the same as dwtest.
library(lmtest)
tab <- as.table(cbind(
'(1)' = df$Y,
'(2)' = df$X,
'(3)' = resid(fit),
'(4)' = c(NA, resid(fit)[t] - resid(fit)[t-1]),
'(5)' = c(NA, (resid(fit)[t] - resid(fit)[t-1])^2),
'(6)' = resid(fit)^2))
round(tab, 4)
## (1) (2) (3) (4) (5) (6)
## 1 20.9600 127.3000 -0.0261 0.0007
## 2 21.4000 130.0000 -0.0620 -0.0360 0.0013 0.0038
## 3 21.9600 132.7000 0.0220 0.0840 0.0071 0.0005
## 4 21.5200 129.4000 0.1638 0.1417 0.0201 0.0268
## 5 22.3900 135.0000 0.0466 -0.1172 0.0137 0.0022
## 6 22.7600 137.1000 0.0464 -0.0002 0.0000 0.0022
## 7 23.4800 141.2000 0.0436 -0.0028 0.0000 0.0019
## 8 23.6600 142.8000 -0.0584 -0.1021 0.0104 0.0034
## 9 24.1000 145.5000 -0.0944 -0.0360 0.0013 0.0089
## 10 24.0100 145.3000 -0.1491 -0.0547 0.0030 0.0222
## 11 24.5400 148.3000 -0.1480 0.0012 0.0000 0.0219
## 12 24.3000 146.4000 -0.0531 0.0949 0.0090 0.0028
## 13 25.0000 150.2000 -0.0229 0.0301 0.0009 0.0005
## 14 25.6400 153.1000 0.1059 0.1288 0.0166 0.0112
## 15 26.3600 157.3000 0.0855 -0.0204 0.0004 0.0073
## 16 26.9800 160.7000 0.1061 0.0206 0.0004 0.0113
## 17 27.5200 164.2000 0.0291 -0.0770 0.0059 0.0008
## 18 27.7800 165.6000 0.0423 0.0132 0.0002 0.0018
## 19 28.2400 168.7000 -0.0442 -0.0865 0.0075 0.0020
## 20 28.7800 171.7000 -0.0330 0.0112 0.0001 0.0011
tab <- (resid(fit)[t] - resid(fit)[t-1])^2 /
sum(resid(fit)^2)
sum(tab) # D = 0.735 (p 488)
## [1] 0.7347
dwtest(fit) # alternative hypotheses defaults to rho > 0
##
## Durbin-Watson test
##
## data: fit
## DW = 0.7347, p-value = 0.0001748
## alternative hypothesis: true autocorrelation is greater than 0
fit.ordinary <- fit # Save for end of chapter
# An NA value is prefixed to each vector to fill out the t+1 = n observations
as.table(cbind(
'(1)' = resid(fit),
'(2)' = c(NA, resid(fit)[t-1]),
'(3)' = c(NA, resid(fit)[t-1] * resid(fit)[t]),
'(4)' = c(NA, resid(fit)[t-1]^2)))
## (1) (2) (3) (4)
## 1 -0.0260519
## 2 -0.0620154 -0.0260519 0.0016156 0.0006787
## 3 0.0220210 -0.0620154 -0.0013656 0.0038459
## 4 0.1637542 0.0220210 0.0036060 0.0004849
## 5 0.0465705 0.1637542 0.0076261 0.0268155
## 6 0.0463766 0.0465705 0.0021598 0.0021688
## 7 0.0436171 0.0463766 0.0020228 0.0021508
## 8 -0.0584354 0.0436171 -0.0025488 0.0019024
## 9 -0.0943990 -0.0584354 0.0055162 0.0034147
## 10 -0.1491425 -0.0943990 0.0140789 0.0089112
## 11 -0.1479909 -0.1491425 0.0220717 0.0222435
## 12 -0.0530536 -0.1479909 0.0078514 0.0219013
## 13 -0.0229282 -0.0530536 0.0012164 0.0028147
## 14 0.1058516 -0.0229282 -0.0024270 0.0005257
## 15 0.0854638 0.1058516 0.0090465 0.0112046
## 16 0.1061022 0.0854638 0.0090679 0.0073041
## 17 0.0291124 0.1061022 0.0030889 0.0112577
## 18 0.0423165 0.0291124 0.0012319 0.0008475
## 19 -0.0441603 0.0423165 -0.0018687 0.0017907
## 20 -0.0330087 -0.0441603 0.0014577 0.0019501
c("rho" = sum(resid(fit)[t-1] * resid(fit)[t]) / sum(resid(fit)[t-1]^2))
## rho
## 0.6312
library(lmtest)
cochrane.orcutt <- function(model)
{
x <- model.matrix(model)[, -1]
y <- model.response(model.frame(model))
e <- resid(model)
n <- length(e)
t <- 2:n
r <- sum(e[t-1] * e[t]) / sum(e[t-1]^2) # (12.22) Step 1
y <- y[t] - r * y[t-1] # (12.18a)
x <- x[t] - r * x[t-1] # (12.18b)
model <- lm(y ~ x) # (12.19) Step 2
model$rho <- r
return(model)
}
fit <- lm(Y ~ X, df)
fit <- cochrane.orcutt(fit)
cbind(df,
rbind(NA, model.frame(fit)))
## Y X y x
## 1 20.96 127.3 NA NA
## 2 21.40 130.0 8.171 49.65
## 3 21.96 132.7 8.453 50.65
## 4 21.52 129.4 7.660 45.64
## 5 22.39 135.0 8.807 53.33
## 6 22.76 137.1 8.628 51.89
## 7 23.48 141.2 9.115 54.67
## 8 23.66 142.8 8.840 53.68
## 9 24.10 145.5 9.167 55.37
## 10 24.01 145.3 8.799 53.47
## 11 24.54 148.3 9.386 56.59
## 12 24.30 146.4 8.811 52.80
## 13 25.00 150.2 9.663 57.80
## 14 25.64 153.1 9.861 58.30
## 15 26.36 157.3 10.177 60.67
## 16 26.98 160.7 10.343 61.42
## 17 27.52 164.2 10.491 62.77
## 18 27.78 165.6 10.410 61.96
## 19 28.24 168.7 10.706 64.18
## 20 28.78 171.7 10.956 65.22
# Generated in function above and stored with model to be used below.
c('rho' = fit$rho)
## rho
## 0.6312
summary(fit)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0970 -0.0568 0.0099 0.0345 0.1250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.39411 0.16723 -2.36 0.031 *
## x 0.17376 0.00296 58.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0672 on 17 degrees of freedom
## Multiple R-squared: 0.995, Adjusted R-squared: 0.995
## F-statistic: 3.45e+03 on 1 and 17 DF, p-value: <2e-16
anova(fit)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 15.57 15.6 3454 <2e-16 ***
## Residuals 17 0.08 0.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Durbin-Watson test on Cochrane-Orcutt model after 1 iteration
dwtest(fit) # DW = 1.6502, fail to reject H0
##
## Durbin-Watson test
##
## data: fit
## DW = 1.65, p-value = 0.1517
## alternative hypothesis: true autocorrelation is greater than 0
# (12.23) Convert transformed coefficients back to original
cat("Y = ", coef(fit)[1] / (1 - fit$rho), " + ", coef(fit)[2], "X\n", sep = "")
## Y = -1.069 + 0.1738X
fit.cochrane <- fit # Save for end of chapter
rho <- c(0.1, 0.3, 0.5, 0.7, 0.9, 0.92, 0.94, 0.95, 0.96, 0.97, 0.98)
hildreth.lu <- function(rho, model)
{
x <- model.matrix(model)[, -1]
y <- model.response(model.frame(model))
n <- length(y)
t <- 2:n
y <- y[t] - rho * y[t-1]
x <- x[t] - rho * x[t-1]
return(lm(y ~ x))
}
fit <- lm(Y ~ X, df)
# Minimum SSE = 0.07167118 for rho = 0.96. See plot below.
tab <- data.frame('rho' = rho,
'SSE' = sapply(rho, function(r) {deviance(hildreth.lu(r, fit))}))
round(tab, 4)
## rho SSE
## 1 0.10 0.1170
## 2 0.30 0.0938
## 3 0.50 0.0805
## 4 0.70 0.0758
## 5 0.90 0.0728
## 6 0.92 0.0723
## 7 0.94 0.0718
## 8 0.95 0.0717
## 9 0.96 0.0717
## 10 0.97 0.0717
## 11 0.98 0.0720
plot(SSE ~ rho, tab, type = 'l')
abline(v = tab[tab$SSE == min(tab$SSE), 'rho'], lty = 3)
fit <- hildreth.lu(0.96, fit)
summary(fit)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1149 -0.0440 0.0111 0.0397 0.1395
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07117 0.05797 1.23 0.24
## x 0.16045 0.00684 23.46 2.2e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0649 on 17 degrees of freedom
## Multiple R-squared: 0.97, Adjusted R-squared: 0.968
## F-statistic: 550 on 1 and 17 DF, p-value: 2.18e-14
dwtest(fit)
##
## Durbin-Watson test
##
## data: fit
## DW = 1.725, p-value = 0.2816
## alternative hypothesis: true autocorrelation is greater than 0
# (12.28) Convert transformed coefficients back to original
cat("Y = ", coef(fit)[1] / (1 - 0.96), " + ", coef(fit)[2], "X", sep = "")
## Y = 1.779 + 0.1605X
fit.hildreth <- fit # Save for end of chapter
# One could also just do lm(diff(Y) ~ 0 + diff(X), df), but the coefficients would be oddly named diff(X).
fit <- lm(y ~ x - 1, data.frame(y = diff(df$Y), x = diff(df$X)))
cbind(
'(1)' = df$Y,
'(2)' = df$X,
'(3)' = c(NA, diff(df$Y)),
'(4)' = c(NA, diff(df$X)))
## (1) (2) (3) (4)
## [1,] 20.96 127.3 NA NA
## [2,] 21.40 130.0 0.44 2.7
## [3,] 21.96 132.7 0.56 2.7
## [4,] 21.52 129.4 -0.44 -3.3
## [5,] 22.39 135.0 0.87 5.6
## [6,] 22.76 137.1 0.37 2.1
## [7,] 23.48 141.2 0.72 4.1
## [8,] 23.66 142.8 0.18 1.6
## [9,] 24.10 145.5 0.44 2.7
## [10,] 24.01 145.3 -0.09 -0.2
## [11,] 24.54 148.3 0.53 3.0
## [12,] 24.30 146.4 -0.24 -1.9
## [13,] 25.00 150.2 0.70 3.8
## [14,] 25.64 153.1 0.64 2.9
## [15,] 26.36 157.3 0.72 4.2
## [16,] 26.98 160.7 0.62 3.4
## [17,] 27.52 164.2 0.54 3.5
## [18,] 27.78 165.6 0.26 1.4
## [19,] 28.24 168.7 0.46 3.1
## [20,] 28.78 171.7 0.54 3.0
summary(fit)
##
## Call:
## lm(formula = y ~ x - 1, data = data.frame(y = diff(df$Y), x = diff(df$X)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0896 -0.0323 0.0241 0.0534 0.1514
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x 0.1685 0.0051 33.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0694 on 18 degrees of freedom
## Multiple R-squared: 0.984, Adjusted R-squared: 0.983
## F-statistic: 1.09e+03 on 1 and 18 DF, p-value: <2e-16
anova(fit)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 5.26 5.26 1093 <2e-16 ***
## Residuals 18 0.09 0.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dwtest(fit)
##
## Durbin-Watson test
##
## data: fit
## DW = 1.739, p-value = 0.3279
## alternative hypothesis: true autocorrelation is greater than 0
# (12.33) Convert transformed coefficient back to original
cat("Y = ", mean(df$Y) - coef(fit)*mean(df$X), " + ", coef(fit), "X\n", sep = "")
## Y = -0.304 + 0.1685X
fit.difference <- fit # Save for end of chapter
This table will be ignored for lack of pedagogical value. The procedures have already been discussed, and it would be a tedious exercise to build the nice summary table as it is displayed. Instead, the various models have bee saved up to this point for any exploration required, save for the rho values, which are either implicit or explored during those sections. Below I show just a simple approach to obtaining some of that information. I leave the MSE column for the interested reader to figure out.
# First column of slope coefficients
c('Cochrane' = c(coef(fit.cochrane)[2]),
'Hildreth' = c(coef(fit.hildreth)[2]),
'Difference' = c(coef(fit.difference)),
'Ordinary' = c(coef(fit.ordinary)[2]))
## Cochrane.x Hildreth.x Difference.x Ordinary.X
## 0.1738 0.1605 0.1685 0.1763
# Second column of sloope coefficient standard errors
c('Cochrane' = round(summary(fit.cochrane)[['coefficients']][2, 2], 4),
'Hildreth' = round(summary(fit.hildreth)[['coefficients']][2, 2], 4),
'Difference' = round(summary(fit.difference)[['coefficients']][1, 2], 4),
'Ordinary' = round(summary(fit.ordinary)[['coefficients']][2, 2], 4))
## Cochrane Hildreth Difference Ordinary
## 0.0030 0.0068 0.0051 0.0014
This forecast was manually calculated from the Cochrane-Orcutt modified (but not transformed) regression model. The steps lack pedagogical value for using R since we're just using it as a calculator in that case.
This chapter was a very brief introduction to time series analysis with R. I suggest looking at Robert H. Shumway and David S. Stoffer, Time Series Analysis and Its Applicatioins. For details visit the text website.