# pdf("F:\\sta4210\\CH12TA02.pdf") blaisdell <- read.table("http://www.stat.ufl.edu/~winner/sta4210/data/CH12TA02.txt", header=F, col.names=c("Y","X")) attach(blaisdell) n <- length(Y) blais.mod1 <- lm(Y~X) summary(blais.mod1) anova(blais.mod1) SSE.mod1 <- deviance(blais.mod1) e.mod1 <- residuals(blais.mod1) plot(e.mod1,type="o",xlab="Time",ylab="Residual") DW1.mod1 <- 0 CO1.mod1 <- 0 CO2.mod1 <- 0 for (t in 2:n) { DW1.mod1 <- DW1.mod1 + (e.mod1[t] - e.mod1[t-1])^2 CO1.mod1 <- CO1.mod1 + e.mod1[t-1]*e.mod1[t] CO2.mod1 <- CO2.mod1 + (e.mod1[t-1])^2 } (DW.mod1 <- DW1.mod1/SSE.mod1) (COr.mod1 <- CO1.mod1/CO2.mod1) ################### Cochrane-Orcutt Method Yt <- Y[2:n] - (COr.mod1 * Y[1:(n-1)]) Xt <- X[2:n] - (COr.mod1 * X[1:(n-1)]) print(cbind(Yt,Xt)) CO.mod1t <- lm(Yt ~ Xt) (CO.mod1t.sum <- summary(CO.mod1t)) anova(CO.mod1t) e.CO1t <- residuals(CO.mod1t) SSE.CO1t <- deviance(CO.mod1t) DW1.CO1t <- 0 for (t in 2:(n-1)) DW1.CO1t <- DW1.CO1t + (e.CO1t[t] - e.CO1t[t-1])^2 (DW.CO1t <- DW1.CO1t/SSE.CO1t) (b0.CO1o <- coef(CO.mod1t.sum)[1,1]/(1-COr.mod1)) (b1.CO1o <- coef(CO.mod1t.sum)[2,1]) (s.b0.CO1o <- coef(CO.mod1t.sum)[1,2]/(1-COr.mod1)) (s.b1.CO1o <- coef(CO.mod1t.sum)[2,2]) ##################### Hildreth-Lu Method ######### SSE.r <- matrix(rep(0,2*length(seq(0.01,0.99,0.01))),ncol=2) r.count <- 0 HL.r <- 0 min.SSE <- 10000000 for (r in seq(0.01,0.99,0.01)) { r.count <- r.count + 1 Yt <- Y[2:n] - (r * Y[1:(n-1)]) Xt <- X[2:n] - (r * X[1:(n-1)]) SSE.r[r.count,1] <- r SSE.r[r.count,2] <- deviance(lm(Yt ~ Xt)) if (deviance(lm(Yt ~ Xt)) < min.SSE) { min.SSE <- deviance(lm(Yt ~ Xt)) HL.r <- r } } HL.r SSE.r Yt <- Y[2:n] - (HL.r * Y[1:(n-1)]) Xt <- X[2:n] - (HL.r * X[1:(n-1)]) print(cbind(Yt,Xt)) HL.mod1t <- lm(Yt ~ Xt) (HL.mod1t.sum <- summary(HL.mod1t)) anova(HL.mod1t) e.HL1t <- residuals(HL.mod1t) SSE.HL1t <- deviance(HL.mod1t) DW1.HL1t <- 0 for (t in 2:(n-1)) DW1.HL1t <- DW1.HL1t + (e.HL1t[t] - e.HL1t[t-1])^2 (DW.HL1t <- DW1.HL1t/SSE.HL1t) (b0.HL1o <- coef(HL.mod1t.sum)[1,1]/(1-HL.r)) (b1.HL1o <- coef(HL.mod1t.sum)[2,1]) (s.b0.HL1o <- coef(HL.mod1t.sum)[1,2]/(1-HL.r)) (s.b1.HL1o <- coef(HL.mod1t.sum)[2,2]) ############## First Differences Method Yt <- Y[2:n]-Y[1:(n-1)] Xt <- X[2:n]-X[1:(n-1)] FD.mod1t <- lm(Yt ~ -1 + Xt) (FD.mod1t.sum <- summary(FD.mod1t)) anova(FD.mod1t) (b0.FD1o <- mean(Y) - (coef(FD.mod1t.sum)[1,1]*mean(X))) (b1.FD1o <- coef(FD.mod1t.sum)[1,1]) (s.b1.FD1o <- coef(FD.mod1t.sum)[1,2]) ####### Generating next period Forecast (X_21=175.3) X_21 <- 175.3 ### Cochrane-Orthcutt Yhat_last.CO1 <- b0.CO1o + b1.CO1o*X[n] e_last.CO1 <- Y[n]-Yhat_last.CO1 (F21_CO <- b0.CO1o + b1.CO1o*X_21 + COr.mod1*e_last.CO1) Yt_CO <- Y[2:n] - COr.mod1*Y[1:(n-1)] Xt_CO <- X[2:n] - COr.mod1*X[1:(n-1)] Xt_21_CO <- X_21-COr.mod1*X[20] MSE_COt <- deviance(lm(Yt_CO ~ Xt_CO))/df.residual(lm(Yt_CO ~ Xt_CO)) n_CO <- length(Yt_CO) SS_XX_CO <- sum((Xt_CO-mean(Xt_CO))^2) s_pred_CO <- sqrt(MSE_COt * (1 + (1/n_CO) + (((Xt_21_CO-mean(Xt_CO))^2)/SS_XX_CO))) (Y21_PI_CO <- cbind(F21_CO-qt(0.975,df.residual(lm(Yt_CO ~ Xt_CO)))*s_pred_CO, F21_CO+qt(0.975,df.residual(lm(Yt_CO ~ Xt_CO)))*s_pred_CO)) ### Hildreth-Liu Yhat_last.HL1 <- b0.HL1o + b1.HL1o*X[n] e_last.HL1 <- Y[n]-Yhat_last.HL1 (F21_HL <- b0.HL1o + b1.HL1o*X_21 + HL.r*e_last.HL1) Yt_HL <- Y[2:n] - HL.r*Y[1:(n-1)] Xt_HL <- X[2:n] - HL.r*X[1:(n-1)] Xt_21_HL <- X_21-HL.r*X[20] MSE_HLt <- deviance(lm(Yt_HL ~ Xt_HL))/df.residual(lm(Yt_HL ~ Xt_HL)) n_HL <- length(Yt_HL) SS_XX_HL <- sum((Xt_HL-mean(Xt_HL))^2) s_pred_HL <- sqrt(MSE_HLt * (1 + (1/n_HL) + (((Xt_21_HL-mean(Xt_HL))^2)/SS_XX_HL))) (Y21_PI_HL <- cbind(F21_HL-qt(0.975,df.residual(lm(Yt_HL ~ Xt_HL)))*s_pred_HL, F21_HL+qt(0.975,df.residual(lm(Yt_HL ~ Xt_HL)))*s_pred_HL)) ### First Difference Yhat_last.FD1 <- b0.FD1o + b1.FD1o*X[n] e_last.FD1 <- Y[n]-Yhat_last.FD1 (F21_FD <- b0.FD1o + b1.FD1o*X_21 + 1.0*e_last.FD1) Yt_FD <- Y[2:n] - Y[1:(n-1)] Xt_FD <- X[2:n] - X[1:(n-1)] Xt_21_FD <- X_21-X[20] MSE_FDt <- deviance(lm(Yt_FD ~ -1 + Xt_FD))/df.residual(lm(Yt_FD ~ -1 + Xt_FD)) n_FD <- length(Yt_FD) SS_XX_FD <- sum((Xt_FD-mean(Xt_FD))^2) s_pred_FD <- sqrt(MSE_FDt * (1 + (1/n_FD) + (((Xt_21_FD-mean(Xt_FD))^2)/SS_XX_FD))) (Y21_PI_FD <- cbind(F21_FD-qt(0.975,df.residual(lm(Yt_FD ~ -1 + Xt_FD)))*s_pred_FD, F21_FD+qt(0.975,df.residual(lm(Yt_FD ~ Xt_FD)))*s_pred_FD)) # dev.off()