pdf("F:\\Rmisc\\graphs\\ac_casino.pdf") ac_casino <- read.table("http://www.stat.ufl.edu/~winner/sta4210/ac_casino.dat",header=F, col.names=c("year","profit","slot","table")) attach(ac_casino) n <- length(year) casino.mod1 <- lm(profit ~ slot + table) summary(casino.mod1) anova(casino.mod1) SSE.mod1 <- deviance(casino.mod1) e.mod1 <- residuals(casino.mod1) plot(e.mod1,type="o",xlab="Time",ylab="Residual") #### Compute the Durbin-Watson Stat and Cochrane-Orcutt Autocorrelation estimate 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 ###### Transformed profit, slot, table Yt <- profit[2:n] - (COr.mod1 * profit[1:(n-1)]) X1t <- slot[2:n] - (COr.mod1 * slot[1:(n-1)]) X2t <- table[2:n] - (COr.mod1 * table[1:(n-1)]) CO.mod1t <- lm(Yt ~ X1t + X2t) (CO.mod1t.sum <- summary(CO.mod1t)) anova(CO.mod1t) e.CO1t <- residuals(CO.mod1t) SSE.CO1t <- deviance(CO.mod1t) ###### Compute the D-W Stat for transformed model 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) ###### Back-Transform the Coefficients (b0.CO1o <- coef(CO.mod1t.sum)[1,1]/(1-COr.mod1)) (b1.CO1o <- coef(CO.mod1t.sum)[2,1]) (b2.CO1o <- coef(CO.mod1t.sum)[3,1]) (s.b0.CO1o <- coef(CO.mod1t.sum)[1,2]/(1-COr.mod1)) (s.b1.CO1o <- coef(CO.mod1t.sum)[2,2]) (s.b2.CO1o <- coef(CO.mod1t.sum)[3,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 <- profit[2:n] - (r * profit[1:(n-1)]) X1t <- slot[2:n] - (r * slot[1:(n-1)]) X2t <- table[2:n] - (r * table[1:(n-1)]) SSE.r[r.count,1] <- r SSE.r[r.count,2] <- deviance(lm(Yt ~ X1t + X2t)) if (deviance(lm(Yt ~ X1t + X2t)) < min.SSE) { min.SSE <- deviance(lm(Yt ~ X1t + X2t)) HL.r <- r } } HL.r SSE.r Yt <- profit[2:n] - (HL.r * profit[1:(n-1)]) X1t <- slot[2:n] - (HL.r * slot[1:(n-1)]) X2t <- table[2:n] - (HL.r * table[1:(n-1)]) HL.mod1t <- lm(Yt ~ X1t + X2t) (HL.mod1t.sum <- summary(HL.mod1t)) anova(HL.mod1t) e.HL1t <- residuals(HL.mod1t) SSE.HL1t <- deviance(HL.mod1t) ##### Compute D-W Stat 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) ##### Back-transform the coefficients (b0.HL1o <- coef(HL.mod1t.sum)[1,1]/(1-HL.r)) (b1.HL1o <- coef(HL.mod1t.sum)[2,1]) (b2.HL1o <- coef(HL.mod1t.sum)[3,1]) (s.b0.HL1o <- coef(HL.mod1t.sum)[1,2]/(1-HL.r)) (s.b1.HL1o <- coef(HL.mod1t.sum)[2,2]) (s.b2.HL1o <- coef(HL.mod1t.sum)[3,2]) ############## First Differences Method Yt <- profit[2:n]-profit[1:(n-1)] X1t <- slot[2:n]-slot[1:(n-1)] X2t <- table[2:n]-table[1:(n-1)] FD.mod1t <- lm(Yt ~ -1 + X1t + X2t) (FD.mod1t.sum <- summary(FD.mod1t)) anova(FD.mod1t) (b0.FD1o <- mean(profit) - (coef(FD.mod1t.sum)[1,1]*mean(slot)) - (coef(FD.mod1t.sum)[2,1]*mean(table))) (b1.FD1o <- coef(FD.mod1t.sum)[1,1]) (b2.FD1o <- coef(FD.mod1t.sum)[2,1]) (s.b1.FD1o <- coef(FD.mod1t.sum)[1,2]) (s.b2.FD1o <- coef(FD.mod1t.sum)[2,2]) dev.off()