> attach(blaisdell)
The following object(s) are masked from 'blaisdell (position 3)':

    X, Y
> 
> n <- length(Y)
> 
> blais.mod1 <- lm(Y~X)
> summary(blais.mod1)

Call:
lm(formula = Y ~ X)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.149142 -0.054399 -0.000454  0.046425  0.163754 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.454750   0.214146  -6.793 2.31e-06 ***
X            0.176283   0.001445 122.017  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.08606 on 18 degrees of freedom
Multiple R-squared: 0.9988,     Adjusted R-squared: 0.9987 
F-statistic: 1.489e+04 on 1 and 18 DF,  p-value: < 2.2e-16 

> anova(blais.mod1)
Analysis of Variance Table

Response: Y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
X          1 110.257 110.257   14888 < 2.2e-16 ***
Residuals 18   0.133   0.007                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 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)
        2 
0.7347256 
> (COr.mod1 <- CO1.mod1/CO2.mod1)
        1 
0.6311636 
> 
> ################### 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))
             Yt       Xt
 [1,]  8.170812 49.65288
 [2,]  8.453100 50.64874
 [3,]  7.659648 45.64460
 [4,]  8.807360 53.32744
 [5,]  8.628248 51.89292
 [6,]  9.114717 54.66748
 [7,]  8.840280 53.67971
 [8,]  9.166670 55.36984
 [9,]  8.798958 53.46570
[10,]  9.385763 56.59193
[11,]  8.811246 52.79844
[12,]  9.662725 57.79765
[13,]  9.860911 58.29923
[14,] 10.176966 60.66886
[15,] 10.342529 61.41797
[16,] 10.491207 62.77202
[17,] 10.410379 61.96294
[18,] 10.706276 64.17931
[19,] 10.955941 65.22271
> 
> CO.mod1t <- lm(Yt ~ Xt)
> (CO.mod1t.sum <- summary(CO.mod1t))

Call:
lm(formula = Yt ~ Xt)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.097039 -0.056815  0.009902  0.034553  0.125048 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.394111   0.167230  -2.357   0.0307 *  
Xt           0.173758   0.002957  58.767   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.06715 on 17 degrees of freedom
Multiple R-squared: 0.9951,     Adjusted R-squared: 0.9948 
F-statistic:  3454 on 1 and 17 DF,  p-value: < 2.2e-16 

> anova(CO.mod1t)
Analysis of Variance Table

Response: Yt
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Xt         1 15.5748 15.5748  3453.6 < 2.2e-16 ***
Residuals 17  0.0767  0.0045                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 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)
       2 
1.650248 
> 
> (b0.CO1o <- coef(CO.mod1t.sum)[1,1]/(1-COr.mod1))
        1 
-1.068524 
> (b1.CO1o <- coef(CO.mod1t.sum)[2,1])
[1] 0.1737583
> 
> (s.b0.CO1o <- coef(CO.mod1t.sum)[1,2]/(1-COr.mod1))
        1 
0.4533986 
> (s.b1.CO1o <- coef(CO.mod1t.sum)[2,2])
[1] 0.002956708
> 
> 
> ##################### 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
[1] 0.96
> SSE.r
       [,1]       [,2]
  [1,] 0.01 0.13082501
  [2,] 0.02 0.12918794
  [3,] 0.03 0.12757684
  [4,] 0.04 0.12599172
  [5,] 0.05 0.12443255
  [6,] 0.06 0.12289932
  [7,] 0.07 0.12139200
  [8,] 0.08 0.11991059
  [9,] 0.09 0.11845505
 [10,] 0.10 0.11702537
 [11,] 0.11 0.11562152
 [12,] 0.12 0.11424349
 [13,] 0.13 0.11289124
 [14,] 0.14 0.11156476
 [15,] 0.15 0.11026401
 [16,] 0.16 0.10898897
 [17,] 0.17 0.10773960
 [18,] 0.18 0.10651588
 [19,] 0.19 0.10531777
 [20,] 0.20 0.10414523
 [21,] 0.21 0.10299823
 [22,] 0.22 0.10187672
 [23,] 0.23 0.10078067
 [24,] 0.24 0.09971003
 [25,] 0.25 0.09866476
 [26,] 0.26 0.09764479
 [27,] 0.27 0.09665008
 [28,] 0.28 0.09568058
 [29,] 0.29 0.09473622
 [30,] 0.30 0.09381693
 [31,] 0.31 0.09292266
 [32,] 0.32 0.09205332
 [33,] 0.33 0.09120884
 [34,] 0.34 0.09038914
 [35,] 0.35 0.08959413
 [36,] 0.36 0.08882372
 [37,] 0.37 0.08807780
 [38,] 0.38 0.08735627
 [39,] 0.39 0.08665901
 [40,] 0.40 0.08598591
 [41,] 0.41 0.08533682
 [42,] 0.42 0.08471162
 [43,] 0.43 0.08411014
 [44,] 0.44 0.08353223
 [45,] 0.45 0.08297771
 [46,] 0.46 0.08244639
 [47,] 0.47 0.08193807
 [48,] 0.48 0.08145253
 [49,] 0.49 0.08098953
 [50,] 0.50 0.08054883
 [51,] 0.51 0.08013015
 [52,] 0.52 0.07973319
 [53,] 0.53 0.07935763
 [54,] 0.54 0.07900313
 [55,] 0.55 0.07866930
 [56,] 0.56 0.07835574
 [57,] 0.57 0.07806201
 [58,] 0.58 0.07778762
 [59,] 0.59 0.07753205
 [60,] 0.60 0.07729474
 [61,] 0.61 0.07707506
 [62,] 0.62 0.07687235
 [63,] 0.63 0.07668586
 [64,] 0.64 0.07651480
 [65,] 0.65 0.07635831
 [66,] 0.66 0.07621543
 [67,] 0.67 0.07608513
 [68,] 0.68 0.07596629
 [69,] 0.69 0.07585769
 [70,] 0.70 0.07575800
 [71,] 0.71 0.07566579
 [72,] 0.72 0.07557949
 [73,] 0.73 0.07549742
 [74,] 0.74 0.07541778
 [75,] 0.75 0.07533862
 [76,] 0.76 0.07525789
 [77,] 0.77 0.07517339
 [78,] 0.78 0.07508285
 [79,] 0.79 0.07498390
 [80,] 0.80 0.07487415
 [81,] 0.81 0.07475122
 [82,] 0.82 0.07461282
 [83,] 0.83 0.07445687
 [84,] 0.84 0.07428163
 [85,] 0.85 0.07408590
 [86,] 0.86 0.07386921
 [87,] 0.87 0.07363215
 [88,] 0.88 0.07337669
 [89,] 0.89 0.07310654
 [90,] 0.90 0.07282763
 [91,] 0.91 0.07254850
 [92,] 0.92 0.07228073
 [93,] 0.93 0.07203920
 [94,] 0.94 0.07184225
 [95,] 0.95 0.07171147
 [96,] 0.96 0.07167118
 [97,] 0.97 0.07174743
 [98,] 0.98 0.07196656
 [99,] 0.99 0.07235338
> 
> Yt <- Y[2:n] - (HL.r * Y[1:(n-1)])
> Xt <- X[2:n] - (HL.r * X[1:(n-1)])
> 
> print(cbind(Yt,Xt))
          Yt     Xt
 [1,] 1.2784  7.792
 [2,] 1.4160  7.900
 [3,] 0.4384  2.008
 [4,] 1.7308 10.776
 [5,] 1.2656  7.500
 [6,] 1.6304  9.584
 [7,] 1.1192  7.248
 [8,] 1.3864  8.412
 [9,] 0.8740  5.620
[10,] 1.4904  8.812
[11,] 0.7416  4.032
[12,] 1.6720  9.656
[13,] 1.6400  8.908
[14,] 1.7456 10.324
[15,] 1.6744  9.692
[16,] 1.6192  9.928
[17,] 1.3608  7.968
[18,] 1.5712  9.724
[19,] 1.6696  9.748
> 
> HL.mod1t <- lm(Yt ~ Xt)
> (HL.mod1t.sum <- summary(HL.mod1t))

Call:
lm(formula = Yt ~ Xt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.11494 -0.04399  0.01113  0.03968  0.13951 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.07117    0.05797   1.228    0.236    
Xt           0.16045    0.00684  23.458 2.18e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.06493 on 17 degrees of freedom
Multiple R-squared:  0.97,      Adjusted R-squared: 0.9683 
F-statistic: 550.3 on 1 and 17 DF,  p-value: 2.177e-14 

> anova(HL.mod1t)
Analysis of Variance Table

Response: Yt
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Xt         1 2.31988 2.31988  550.26 2.177e-14 ***
Residuals 17 0.07167 0.00422                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 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)
       2 
1.725439 
> 
> (b0.HL1o <- coef(HL.mod1t.sum)[1,1]/(1-HL.r))
[1] 1.77933
> (b1.HL1o <- coef(HL.mod1t.sum)[2,1])
[1] 0.1604536
> 
> (s.b0.HL1o <- coef(HL.mod1t.sum)[1,2]/(1-HL.r))
[1] 1.449373
> (s.b1.HL1o <- coef(HL.mod1t.sum)[2,2])
[1] 0.006840127
> 
> 
> ############## 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))

Call:
lm(formula = Yt ~ -1 + Xt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.08958 -0.03231  0.02412  0.05344  0.15139 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
Xt 0.168488   0.005096   33.06   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.06939 on 18 degrees of freedom
Multiple R-squared: 0.9838,     Adjusted R-squared: 0.9829 
F-statistic:  1093 on 1 and 18 DF,  p-value: < 2.2e-16 

> anova(FD.mod1t)
Analysis of Variance Table

Response: Yt
          Df Sum Sq Mean Sq F value    Pr(>F)    
Xt         1 5.2637  5.2637  1093.1 < 2.2e-16 ***
Residuals 18 0.0867  0.0048                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> (b0.FD1o <- mean(Y) - (coef(FD.mod1t.sum)[1,1]*mean(X)))
[1] -0.3040052
> (b1.FD1o <- coef(FD.mod1t.sum)[1,1])
[1] 0.1684878
> 
> 
> (s.b1.FD1o <- coef(FD.mod1t.sum)[1,2])
[1] 0.005096004
> 
> 
> 
> ####### 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)
       1 
29.40028 
> 
> 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))
      [,1]  [,2]
1 29.24056 29.56
> 
> 
> ### 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)
[1] 29.3796
> 
> 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))
         [,1]     [,2]
[1,] 29.23526 29.52394
> 
> 
> ### 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)
[1] 29.38656
> 
> 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))
        [,1]     [,2]
[1,] 29.2356 29.53815
> 
> 
>