w9b: Static and Rolling Predictions
Nao Mimoto - Dept. of Statistics : The University of Akron
1. Ex: Static 50-step AR(2) prediction
Let’s use simulated data and fit AR(2). Then we’ll forecast for 10 steps.
set.seed(77224) # set seed
Y = arima.sim(list(ar = c(.6, .3) ), 200 ) + 22 # Simulate AR(2)
library(forecast) # install.packages("forecast")
Fit1 = Arima(Y, order=c(2,0,0), include.mean=TRUE) # Force to fit AR(2) w/o mean
Fit1## Series: Y
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 mean
## 0.5175 0.4038 21.6454
## s.e. 0.0647 0.0651 0.7748
##
## sigma^2 = 0.8633: log likelihood = -268.45
## AIC=544.9 AICc=545.11 BIC=558.09 Y.forec = predict(Fit1, n.ahead=10) # 10-step prediction of AR(2)
str(Y.forec) # See what's inside the forecast## List of 2
## $ pred: Time-Series [1:10] from 201 to 210: 22.9 23 22.9 22.8 22.8 ...
## $ se : Time-Series [1:10] from 201 to 210: 0.929 1.046 1.218 1.323 1.422 ... plot(forecast(Fit1, 50)) # plot forecast with prediction errors
abline(h=21.6454)
abline(h=Fit1$coef["intercept"])
# Short-term forecast and its forecasting error (PI)
forecast(Fit1, 5)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 22.90607 21.71536 24.09678 21.08503 24.72711
## 202 23.02815 21.68747 24.36883 20.97775 25.07855
## 203 22.87003 21.30896 24.43110 20.48258 25.25748
## 204 22.83751 21.14163 24.53339 20.24388 25.43113
## 205 22.75682 20.93491 24.57873 19.97045 25.54319 X = Y - Fit1$coef["intercept"]
X201_pred = Fit1$coef[1]*X[200] + Fit1$coef[2]*X[199]
Fit1$coef["intercept"]## intercept
## 21.64538 Y201 = Fit1$coef["intercept"] + X201_pred
Y201## intercept
## 22.90607 1.96*sqrt(Fit1$sigma2)## [1] 1.821073 # Long-term forecast and its forecasting error (PI)
forecast(Fit1, 50)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 22.90607 21.71536 24.09678 21.08503 24.72711
## 202 23.02815 21.68747 24.36883 20.97775 25.07855
## 203 22.87003 21.30896 24.43110 20.48258 25.25748
## 204 22.83751 21.14163 24.53339 20.24388 25.43113
## 205 22.75682 20.93491 24.57873 19.97045 25.54319
## 206 22.70194 20.77994 24.62393 19.76250 25.64137
## 207 22.64095 20.63171 24.65019 19.56808 25.71382
## 208 22.58723 20.50400 24.67045 19.40121 25.77325
## 209 22.53480 20.38736 24.68225 19.25057 25.81904
## 210 22.48598 20.28291 24.68905 19.11667 25.85528
## 211 22.43954 20.18792 24.69116 18.99598 25.88310
## 212 22.39579 20.10172 24.68987 18.88730 25.90429
## 213 22.35440 20.02307 24.68574 18.78893 25.91988
## 214 22.31532 19.95122 24.67942 18.69974 25.93090
## 215 22.27838 19.88540 24.67136 18.61864 25.93812
## 216 22.24348 19.82502 24.66194 18.54476 25.94220
## 217 22.21051 19.76951 24.65150 18.47733 25.94368
## 218 22.17935 19.71842 24.64028 18.41568 25.94301
## 219 22.14991 19.67131 24.62850 18.35923 25.94059
## 220 22.12209 19.62783 24.61635 18.30745 25.93673
## 221 22.09581 19.58765 24.60397 18.25991 25.93171
## 222 22.07097 19.55047 24.59148 18.21619 25.92576
## 223 22.04751 19.51603 24.57899 18.17594 25.91908
## 224 22.02534 19.48410 24.56658 18.13885 25.91182
## 225 22.00439 19.45448 24.55430 18.10463 25.90415
## 226 21.98460 19.42696 24.54223 18.07303 25.89616
## 227 21.96589 19.40138 24.53041 18.04381 25.88798
## 228 21.94822 19.37759 24.51886 18.01678 25.87967
## 229 21.93153 19.35544 24.50761 17.99174 25.87131
## 230 21.91575 19.33480 24.49670 17.96853 25.86297
## 231 21.90084 19.31557 24.48612 17.94700 25.85468
## 232 21.88676 19.29762 24.47589 17.92702 25.84650
## 233 21.87345 19.28087 24.46603 17.90845 25.83845
## 234 21.86088 19.26523 24.45652 17.89118 25.83057
## 235 21.84900 19.25061 24.44738 17.87511 25.82288
## 236 21.83777 19.23695 24.43859 17.86016 25.81538
## 237 21.82716 19.22416 24.43016 17.84622 25.80810
## 238 21.81714 19.21220 24.42208 17.83323 25.80105
## 239 21.80767 19.20100 24.41434 17.82111 25.79423
## 240 21.79872 19.19050 24.40694 17.80980 25.78765
## 241 21.79027 19.18067 24.39987 17.79923 25.78130
## 242 21.78228 19.17145 24.39311 17.78936 25.77520
## 243 21.77473 19.16281 24.38666 17.78014 25.76933
## 244 21.76760 19.15470 24.38051 17.77151 25.76370
## 245 21.76086 19.14708 24.37464 17.76343 25.75829
## 246 21.75450 19.13993 24.36906 17.75587 25.75312
## 247 21.74848 19.13322 24.36374 17.74879 25.74817
## 248 21.74280 19.12692 24.35867 17.74215 25.74344
## 249 21.73742 19.12099 24.35386 17.73593 25.73892
## 250 21.73235 19.11542 24.34928 17.73010 25.73460 Fit1$coef["intercept"] # mean estimate## intercept
## 21.64538 21.73235 - 1.96*sd(Y)## [1] 17.75652 21.73235 + 1.96*sd(Y)## [1] 25.708182. Ex: Static 10-step MA(2) prediction
set.seed(77224) # set seed
Y = arima.sim(list(ma = c(.6, .3) ), 200 ) + 8 # Simulate MA(2)
library(forecast) # install.packages("forecast")
Fit1 = Arima(Y, order=c(0,0,2), include.mean=TRUE) # Force to fit AR(2) w/o mean
Fit1## Series: Y
## ARIMA(0,0,2) with non-zero mean
##
## Coefficients:
## ma1 ma2 mean
## 0.6431 0.3133 7.8680
## s.e. 0.0666 0.0662 0.1273
##
## sigma^2 = 0.8659: log likelihood = -268.12
## AIC=544.23 AICc=544.44 BIC=557.43 Y.forec = predict(Fit1, n.ahead=10) # 10-step prediction of AR(2)
plot(forecast(Fit1, 10)) # plot forecast with prediction errors
abline(h=Fit1$coef["intercept"])
# Short-term forecast and its forecasting error (PI)
forecast(Fit1, 5)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 7.457452 6.264932 8.649972 5.633650 9.281254
## 202 7.534149 6.116309 8.951989 5.365750 9.702549
## 203 7.867954 6.401709 9.334200 5.625525 10.110383
## 204 7.867954 6.401709 9.334200 5.625525 10.110383
## 205 7.867954 6.401709 9.334200 5.625525 10.110383 res = Fit1$resid
X201_pred = Fit1$coef[1]*res[200] + Fit1$coef[2]*res[199]
Y201 = Fit1$coef["intercept"] + X201_pred
Y201## intercept
## 7.457452 1.96*sqrt(Fit1$sigma2)## [1] 1.823835 Y201 - 1.96*sqrt(Fit1$sigma2)## intercept
## 5.633617 Y201 + 1.96*sqrt(Fit1$sigma2)## intercept
## 9.281287 # Long-term forecast and its forecasting error (PI)
forecast(Fit1, 10)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 7.457452 6.264932 8.649972 5.633650 9.281254
## 202 7.534149 6.116309 8.951989 5.365750 9.702549
## 203 7.867954 6.401709 9.334200 5.625525 10.110383
## 204 7.867954 6.401709 9.334200 5.625525 10.110383
## 205 7.867954 6.401709 9.334200 5.625525 10.110383
## 206 7.867954 6.401709 9.334200 5.625525 10.110383
## 207 7.867954 6.401709 9.334200 5.625525 10.110383
## 208 7.867954 6.401709 9.334200 5.625525 10.110383
## 209 7.867954 6.401709 9.334200 5.625525 10.110383
## 210 7.867954 6.401709 9.334200 5.625525 10.110383 Fit1$coef["intercept"] # mean estimate## intercept
## 7.867954 Fit1$coef["intercept"] - 1.96*(1.144117)## intercept
## 5.625485 Fit1$coef["intercept"] + 1.96*(1.144117)## intercept
## 10.11042 # 3 ways to estimate sd of Y (1.144)
sd(Y) # sample SD## [1] 1.153121 sqrt(acf(Y, type="covariance")$acf[1]) # sample ACVF(0)
## [1] 1.150234 Fit1$sigma2* (1+Fit1$coef[1]^2 + Fit1$coef[2]^2) # use Tho ACVF of MA(2)## ma1
## 1.3090053. Ex: Rolling 1-step prediction
Suppose you have \(X_1, \ldots, X_{180}\) to begin. - Separate \(Y_t\) into two parts. Pretend that \(X_1,\ldots, X_{100}\) is all you have, and hide the second part \(\{X_{101}, \ldots, X_{180} \}\).
Use \(\{X_1, \ldots, X_{100}\}\) to find fitting model. Say it was AR(2). Produce 1-step predicton \(\{ \hat X_{101}\}\)
Observe \(X_{101}\)
Use \(\{X_2, \ldots, X_{101}\}\) to re-estimate parameters of AR(2). Then produce 1-step predicton \(\{\hat X_{102}\}\)
Observe \(X_{102}\)
Use \(\{X_3, \ldots, X_{102}\}\) to re-estimate parameters of AR(2). Then produce 1-step predicton \(\{\hat X_{103}\}\)
Observe \(X_{103}\)
keep going
Always using past 100 observation (Fixed Window Size) to predict the next obervation
Fixed Start Rolling 1-step prediction will use \(\{X_1, \ldots, X_{105}\}\) to predict \(X_{106}\).
Example: Rolling 1-step prediction (Fixed Window Size, Fixed Model)
set.seed(77224) # set seed
X = arima.sim(list(ar = c(.6, .3)), n=180, sd=1 ) # simulate the dataset
X = ts(X, start=c(432, 1), freq=7)
X1 = window(X, end=time(X)[100], freq=7) # cut TS to 1st 100
X2 = window(X, start=time(X)[100], freq=7) # last 29
plot(X, type="o", col="red") # Entire dataset.
lines(X1, type="o") # Pretend that black is all you have.
library(forecast)
Fit3 = auto.arima(X1, d=0, max.q=0) # Find best AR(p) by AIC
Fit3## Series: X1
## ARIMA(2,0,0) with zero mean
##
## Coefficients:
## ar1 ar2
## 0.5876 0.3392
## s.e. 0.0941 0.0962
##
## sigma^2 = 0.9046: log likelihood = -136.77
## AIC=279.55 AICc=279.8 BIC=287.37 #--- Suppose you are happy with AR(2). Perform rolling 1-step prediction with widnow of 100 obs.
#--- The model AR(2) is fixed as rolled, but parameters are re-estimated each time.
#--- Rolling 1-step Prediction w/ fixed window of 100
library(forecast)
Y = X # Entire data
window.size = 100 # Window size for estimation
order.in = c(2,0,0) # model order
Yhat = Yhat.CIu = Yhat.CIl = 0 # initialize what needs to be saved
for (i in 1:(length(Y)-window.size)) {
# Force to fit AR(2) each time on last 100 obs.
Fit00 = Arima( Y[i:(i+window.size-1)], order=order.in) # <--- AR(2)
Y.h = predict(Fit00, n.ahead=1) # one step prediction
Yhat[i] = Y.h$pred
Yhat.CIu[i] = Yhat[i]+1.96*Y.h$se
Yhat.CIl[i] = Yhat[i]-1.96*Y.h$se
}
# Yhat starts at window.size+1 up to length(Y)
Yhat = ts(Yhat, start=time(Y)[window.size+1], freq=frequency(Y))
Yhat.CIu = ts(Yhat.CIu, start=time(Y)[window.size+1], freq=frequency(Y))
Yhat.CIl = ts(Yhat.CIl, start=time(Y)[window.size+1], freq=frequency(Y))
Y1 = window(Y, end=time(Y)[window.size], freq=frequency(Y))
Y2 = window(Y, start=time(Y)[window.size+1], freq=frequency(Y))
#- Calculate prediction performance
Pred.error = Y2 - Yhat
Pred.rMSE = sqrt( mean( (Pred.error)^2 ) ) # prediction root Mean Squared Error
c( mean(Pred.error), Pred.rMSE ) # Av Pred Error, and pred rMSE## [1] 0.1184882 0.9430092 #- Plot the prediction result with original data
layout(matrix(c(1,1,1,1,1,1,2,3,4), 3, 3, byrow=TRUE)) # 3 plots at once
plot( Y, type="n")
lines( Y1, type="o", main="Rolling 1-step prediction (Red=Actual, Blue=Prediction)" )
lines( Y2, type="o", col="red")
lines( Yhat, type="o", col="blue")
lines( Yhat.CIu, type="l", col="gray30", lty=2)
lines( Yhat.CIl, type="l", col="gray30", lty=2)
# - Plot the prediction error
plot( (window.size+1):length(Y), Pred.error, type="o", main="Prediction Error (Red-Blue)" )
abline(h=0, lty=2)
acf(Pred.error)
layout(1) # reset the layout