Notes from Udemy Forecasting Models
Stephen Ewing 8/3/18
Data setup
Data
library(readr)
data <- read_csv("./Forecasting-Models-Data.txt")Packages Used
library("forecast")
library("tseries")
library("smooth")create a time series. 21 days in the average business month
spy <- ts(data[,2],frequency=21)Delimit training range. Ending after 84 months on the last day in the month
spyt <- window(spy,end=c(84,21))Delimit testing range. Starting on the first day of the 85th month
spyf <- window(spy,start=c(85,1))Training and testing ranges chart
plot(spy,main="SPY 2007-2015",ylab="Price",xlab="Month")
lines(spyt,col="blue")
lines(spyf,col="green")
legend("bottomright",col=c("blue","green"),lty=1,legend=c("Training","Testing"))
Simple Forecasting Methods
Arithmetic Mean Method
A constant forecast equal to the arithmetic mean of previos periods historical data
mean <- meanf(spyt,h=502) #forecasting for 502 periods, length of testing range
plot(mean,main="Arithmetic Mean Method",ylab="Price",xlab="Month",ylim=c(min(spy),max(spy)))
lines(spy,lty=3)
Random Walk Method
Naive or random walk method is a forecast equal to previous period data
Multi-Steps Forecast Forecasting the full testing range
rw1 <- rwf(spyt,h=502) #forecasting for 502 periods, length of testing range
plot(rw1,main="Random Walk Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation Updating the forecast one step at a time but without the resetimation of the best fitting model
rw2 <- lag(spyf,k=-1) # uses testing range spyf, k=-1 means lag one day
plot(spyt,main="Random Walk Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(rw2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="rw2")
Seasonal Random Walk Method
Forecasts equal to previous season data
Multi-Steps Forecast
srw1 <- snaive(spyt,h=502)
plot(srw1,main="Seasonal Random Walk Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
srw2 <- lag(spyf,k=-21) #k=number of periods to lag, so this is one month
plot(spyt,main="Seasonal Random Walk Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(srw2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="srw2")
Random Walk with Drift Method
Forecast equal to previous period data plus the arithmetic mean of previous periods historical data differences
Basically it incorporates drift based on the rate of change observed in the previous periods
Multi-Steps Forecast
rwd1 <- rwf(spyt,drift=TRUE,h=502)
plot(rwd1,main="Random Walk with Drift Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
rwd2 <- lag(spyf,k=-1)+mean(na.exclude(diff(spyt,lag=1))) #the addition of the mean gives this drift
plot(spyt,main="Random Walk with Drift Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(rwd2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="rwd2")
Forecasting Accuracy
main scale-dependent metrics include mean absolute error and root mean squared error mean absolute error: mean of absolute value of residuals root mean squared error: square root of mean of squared residuals
main scale-independent metrics include mean absolute percentage error mean absolute percentage error: mean of absolute value of residuals as a percentage of actual data
Multi-Steps Forecast
accuracy(mean,spyf)[2,1:6] ME RMSE MAE MPE MAPE MASE
76.96225 77.68007 76.96225 40.40627 40.40627 17.78549 accuracy(rw1,spyf)[2,1:6] ME RMSE MAE MPE MAPE MASE
19.092649 21.806753 19.338467 9.783360 9.932039 4.468997 accuracy(srw1,spyf)[2,1:6] ME RMSE MAE MPE MAPE MASE
21.121117 23.683972 21.314148 10.858796 10.975791 4.925564 accuracy(rwd1,spyf)[2,1:6] ME RMSE MAE MPE MAPE MASE
10.999963 13.119287 11.368150 5.618607 5.837789 2.627107 One-Step Forecast without Re-Estimation
accuracy(rw2,spyf)[1,1:5] ME RMSE MAE MPE MAPE
0.05541398 1.61096422 1.17425250 0.02655375 0.62091878 accuracy(srw2,spyf)[1,1:5] ME RMSE MAE MPE MAPE
1.2892505 5.8805092 4.3931968 0.6551872 2.3143711 accuracy(rwd2,spyf)[1,1:5] ME RMSE MAE MPE MAPE
0.023236303 1.610178538 1.171930665 0.009530051 0.619652447 Moving Averages and Exponential Smoothing Methods
Exponential smoothing methods consist of flattening time series data
Simple moving averages consist of forecast based on previous periods data with equal weighted influence.
Weighted moving averages consist of forecasts based on previous periods data with linearly decaying influence the older the become.
Exponential moving averages consist of forecasts based on previous periods data with exponentially decaying influence the older the become.
Exponential smoothing methods are a special case of exponential moving averages with notation ETS (error, trend, seasonality) where each can be none (N), additive (A), additive damped (Ad), multiplicative (M) or multiplicative damped (Md)
Simple Moving Average SMA
Forecast equal to the arithmetic mean of previous periods data
Multi-Steps Forecast
sma1 <- sma(spyt,21,h=502)
plot(spyt,main="Simple Moving Average SMA Method 1",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(sma1$forecast,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="sma1")
One-Step Forecast without Re-Estimation
sma2 <- sma(spyf,21)$fitted
plot(spyt,main="Simple Moving Average SMA Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(sma2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="sma2")
Brown Simple Exponential Smoothing ETS(A,N,N)
No trend or seasonal smoothing
Multi-Steps Forecast
brown1 <- ses(spyt,h=502)
plot(brown1,main="Brown Simple Exponential Smoothing ETS(A,N,N) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
brown2t <- ets(spyt,model="ANN",damped=F)
brown2t## ETS(A,N,N)
##
## Call:
## ets(y = spyt, model = "ANN", damped = F)
##
## Smoothing parameters:
## alpha = 0.903
##
## Initial states:
## l = 113.7761
##
## sigma: 1.3852
##
## AIC AICc BIC
## 14339.99 14340.00 14356.41brown2 <- ets(spyf,model="ANN",damped=F,alpha=brown2t$par[1])
brown2 <- brown2$fitted
plot(spyt,main="Brown Simple Exponential Smoothing ETS(A,N,N) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(brown2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="brown2")
Holt Linear Trend Method ETS(A,A,N)
Forecast with a linear trend pattern
Multi-Steps Forecast
holt1 <- holt(spyt,h=502)
plot(holt1,main="Holt Linear Trend ETS(A,A,N) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
holt2t <- ets(spyt,model="AAN",damped=F)
holt2t## ETS(A,A,N)
##
## Call:
## ets(y = spyt, model = "AAN", damped = F)
##
## Smoothing parameters:
## alpha = 0.902
## beta = 1e-04
##
## Initial states:
## l = 113.231
## b = 0.032
##
## sigma: 1.3856
##
## AIC AICc BIC
## 14343.13 14343.17 14370.51holt2 <- ets(spyf,model="AAN",damped=F,alpha=holt2t$par[1],beta=holt2t$par[2])
holt2 <- holt2$fitted
plot(spyt,main="Holt Linear Trend ETS(A,A,N) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(holt2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="holt2")
Exponential Trend Method ETS(A,M,N)
Forecast with an exponential trend pattern
Multi-Steps Forecast
exp1 <- holt(spyt,h=502,exponential=T)
plot(exp1,main="Exponential Trend ETS(A,M,N) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)One-Step Forecast without Re-Estimation
exp2t <- es(spyt,model="AMN")
exp2t## Time elapsed: 1.1 seconds
## Model estimated: ETS(AMN)
## Persistence vector g:
## alpha beta
## [1,] 0.9 0.002
## Initial values were optimised.
## 5 parameters were estimated in the process
## Residuals standard deviation: 1.386
## Cost function type: MSE; Cost function value: 2
##
## Information criteria:
## AIC AICc BIC
## 6161.483 6161.517 6188.860exp2 <- es(spyf,model="AMN",alpha=exp2t$persistence[1],beta=exp2t$persistence[2]) 
exp2 <- exp2$fitted
plot(spyt,main="Exponential Trend ETS(A,M,N) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(exp2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="exp2")
Gardner Additive Damped Trend Method ETS(A,Ad,N)
Multi-Steps Forecast
gardner1 <- holt(spyt,h=502,damped=T)
plot(gardner1,main="Gardner Additive Damped Trend ETS(A,Ad,N) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
gardner2t <- ets(spyt,model="AAN",damped=T)
gardner2t## ETS(A,Ad,N)
##
## Call:
## ets(y = spyt, model = "AAN", damped = T)
##
## Smoothing parameters:
## alpha = 0.8999
## beta = 0.0014
## phi = 0.98
##
## Initial states:
## l = 113.1943
## b = 0.1365
##
## sigma: 1.3863
##
## AIC AICc BIC
## 14345.80 14345.85 14378.65gardner2 <- ets(spyf,model="AAN",damped=T,alpha=gardner2t$par[1],beta=gardner2t$par[2],phi=gardner2t$par[3])
gardner2 <- gardner2$fitted
plot(spyt,main="Gardner Additive Damped Trend ETS(A,Ad,N) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(gardner2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="gardner2")
Taylor Multiplicative Damped Trend Method ETS(A,Md,N)
Multi-Steps Forecast
taylor1 <- holt(spyt,h=502,exponential=T,damped=T)
plot(taylor1,main="Taylor Multiplicative Damped Trend ETS(A,Md,N) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
taylor2t <- es(spyt,model="AMdN")
taylor2t## Time elapsed: 1.43 seconds
## Model estimated: ETS(AMdN)
## Persistence vector g:
## alpha beta
## [1,] 0.898 0.398
## Damping parameter: 0.013
## Initial values were optimised.
## 6 parameters were estimated in the process
## Residuals standard deviation: 1.387
## Cost function type: MSE; Cost function value: 2
##
## Information criteria:
## AIC AICc BIC
## 6165.496 6165.544 6198.348taylor2 <- es(spyf,model="AMdN",alpha=taylor2t$persistence[1],beta=taylor2t$persistence[2],phi=taylor2t$phi) 
taylor2 <- taylor2$fitted
plot(spyt,main="Taylor Multiplicative Damped Trend ETS(A,Md,N) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(taylor2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="taylor2")
Holt-Winters Additive Method ETS(A,A,A)
Multi-Steps Forecast
hwa1 <- hw(spyt,h=502,seasonal="additive")
plot(hwa1,main="Holt-Winters Additive ETS(A,A,A) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
hwa2t <- ets(spyt,model="AAA",damped=F)
hwa2t## ETS(A,A,A)
##
## Call:
## ets(y = spyt, model = "AAA", damped = F)
##
## Smoothing parameters:
## alpha = 0.8632
## beta = 0.0019
## gamma = 3e-04
##
## Initial states:
## l = 114.8671
## b = -0.0153
## s = 0.1028 0.1121 -0.1641 -0.2105 -0.4711 -0.1812
## -0.0644 -0.0081 0.2191 0.091 0.1776 0.05 0.0073 0.1239 -0.0575 0.0749 0.2443 0.3591 -0.0332 -0.2352 -0.1368
##
## sigma: 1.3944
##
## AIC AICc BIC
## 14386.14 14386.94 14528.50hwa2 <- ets(spyf,model="AAA",damped=F,alpha=hwa2t$par[1],beta=hwa2t$par[2],gamma=hwa2t$par[3])
hwa2 <- hwa2$fitted
plot(spyt,main="Holt-Winters Additive ETS(A,A,A) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(hwa2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="hwa2")
Holt-Winters Multiplicative Method ETS(A,A,M)
Multi-Steps Forecast
hwm1 <- hw(spyt,h=502,seasonal="multiplicative")
plot(hwm1,main="Holt-Winters Multiplicative ETS(A,A,M) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
hwm2t <- es(spyt,model="AAM")
hwm2t## Time elapsed: 1.13 seconds
## Model estimated: ETS(AAM)
## Persistence vector g:
## alpha beta gamma
## [1,] 0.874 0.006 0.013
## Initial values were optimised.
## 27 parameters were estimated in the process
## Residuals standard deviation: 1.405
## Cost function type: MSE; Cost function value: 2
##
## Information criteria:
## AIC AICc BIC
## 6231.300 6232.171 6379.135hwm2 <- es(spyf,model="AAM",alpha=hwm2t$persistence[1],beta=hwm2t$persistence[2],gamma=hwm2t$persistance[3]) 
hwm2 <- hwm2$fitted
plot(spyt,main="Holt-Winters Multiplicative ETS(A,A,M) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(hwm2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="hwm2")
Holt-Winters Damped Method ETS(A,Ad,M)
Multi-Steps Forecast
hwd1 <- hw(spyt,h=502,seasonal="multiplicative",damped=T)
plot(hwd1,main="Holt-Winters Damped ETS(A,Ad,M) Method 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
hwd2t <- es(spyt,model="AAdM")
hwd2t## Time elapsed: 1.02 seconds
## Model estimated: ETS(AAdM)
## Persistence vector g:
## alpha beta gamma
## [1,] 0.853 0.051 0.014
## Damping parameter: 0.931
## Initial values were optimised.
## 28 parameters were estimated in the process
## Residuals standard deviation: 1.415
## Cost function type: MSE; Cost function value: 2
##
## Information criteria:
## AIC AICc BIC
## 6258.061 6258.997 6411.371hwd2 <- es(spyf,model="AAdM",alpha=hwd2t$persistence[1],beta=hwd2t$persistence[2],gamma=hwd2t$persistance[3],phi=hwd2t$phi) 
hwd2 <- hwd2$fitted
plot(spyt,main="Holt-Winters Damped ETS(A,Ad,M) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(hwd2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="hwd2")
Method Selection
Multi-Steps Forecast
sbest1t <- ets(spyt)
sbest1t## ETS(A,N,N)
##
## Call:
## ets(y = spyt)
##
## Smoothing parameters:
## alpha = 0.903
##
## Initial states:
## l = 113.7761
##
## sigma: 1.3852
##
## AIC AICc BIC
## 14339.99 14340.00 14356.41sbest1 <- forecast(sbest1t,h=502)
plot(sbest1,main="Brown Simple Exponential Smoothing ETS(A,N,N) Method 1",ylab="Price",xlab="Month",ylim=c(min(spy),max(sbest1$upper)))
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
sbest2 <- ets(spyf,model="ANN",damped=F,alpha=sbest1t$par[1])
sbest2 <- sbest2$fitted
plot(spyt,main="Brown Simple Exponential Smoothing ETS(A,N,N) Method 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(sbest2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="sbest2")
Method Forecasting Accuracy
Multi-Steps Forecast
accuracy(sma1$forecast,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 20.10903 22.69469 20.29319 10.32171 10.43336accuracy(brown1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.075064 21.791358 19.322143 9.774055 9.923488 4.465224accuracy(holt1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 10.829943 12.950342 11.208678 5.530797 5.755959 2.590254#accuracy(exp1,spyf)[2,1:6]
accuracy(gardner1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 18.546415 21.287171 18.808521 9.497056 9.655471 4.346530accuracy(taylor1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.061782 21.779130 19.309372 9.767066 9.916806 4.462273accuracy(hwa1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## -2.919532 7.654723 5.096729 -1.544089 2.675100 1.177822accuracy(hwm1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## -34.510508 41.162978 34.510508 -17.814962 17.814962 7.975159accuracy(hwd1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 17.977088 20.762987 18.256869 9.197784 9.366738 4.219046accuracy(sbest1,spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.075142 21.791426 19.322215 9.774096 9.923526 4.465241One-Step Forecast without Re-Estimation
accuracy(sma2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.6150055 3.6465211 2.7800489 0.3054230 1.4753627accuracy(brown2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.06148384 1.62273900 1.17609739 0.02942017 0.62214308accuracy(holt2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## -0.0006216599 1.6218572332 1.1727153895 -0.0034269723 0.6203194683accuracy(exp2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.008638834 1.608678645 1.169324923 0.001890390 0.618264094accuracy(gardner2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.04487804 1.62375935 1.17389510 0.02005779 0.62086929accuracy(taylor2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.05442528 1.60857149 1.17214020 0.02607466 0.61976984accuracy(hwa2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## -0.02414681 1.60065195 1.16780883 -0.01590531 0.61768071accuracy(hwm2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## -0.02534185 1.63474421 1.20669416 -0.01603413 0.63814909accuracy(hwd2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.019762077 1.666849577 1.238133085 0.008224394 0.653851616accuracy(sbest2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.06148384 1.62273900 1.17609739 0.02942017 0.62214308Auto Regressive Integrated Moving Average Models
First Order Trend Stationarity
Unit Root Tests
Augmented Dickey-Fuller Test ADF
adf.test(spyt,alternative="stationary")##
## Augmented Dickey-Fuller Test
##
## data: spyt
## Dickey-Fuller = -0.9315, Lag order = 12, p-value = 0.9498
## alternative hypothesis: stationaryPhillips-Perron Test
pp.test(spyt,alternative="stationary")##
## Phillips-Perron Unit Root Test
##
## data: spyt
## Dickey-Fuller Z(alpha) = -2.4147, Truncation lag parameter = 8,
## p-value = 0.9568
## alternative hypothesis: stationaryTime Series Differencing
plot(spy, type="l",main="SPY 2007-2015",ylab="Price",xlab="Month")
plot(diff(spy), type="l",main="SPY 2007-2015",ylab="Price Differences",xlab="Month")
Augmented Dickey-Fuller Test ADF
adf.test(diff(spyt),alternative="stationary")##
## Augmented Dickey-Fuller Test
##
## data: diff(spyt)
## Dickey-Fuller = -11.188, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationaryPhillips-Perron Test PP
pp.test(diff(spyt),alternative="stationary")##
## Phillips-Perron Unit Root Test
##
## data: diff(spyt)
## Dickey-Fuller Z(alpha) = -1782.4, Truncation lag parameter = 8,
## p-value = 0.01
## alternative hypothesis: stationaryARIMA Models Specification
Normal and Partial Autocorrelation Functions ACF & PACF
acf(diff(spyt))
pacf(diff(spyt))
Random Walk Model ARIMA(0,1,0) without constant
Multi-Steps Forecast
arw1 <- Arima(spyt,order=c(0,1,0))
arw1## Series: spyt
## ARIMA(0,1,0)
##
## sigma^2 estimated as 1.934: log likelihood=-3083.01
## AIC=6168.02 AICc=6168.03 BIC=6173.5plot(forecast(arw1,h=502),main="Random Walk ARIMA(0,1,0) Without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
arw2 <- Arima(spyf,model=arw1)$fitted
plot(spyt,main="Random Walk ARIMA(0,1,0) Without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(arw2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="arw2")
Random Walk with Drift Model ARIMA(0,1,0) with constant
Multi-Steps Forecast
arwd1 <- Arima(spyt,order=c(0,1,0),include.constant=T)
arwd1## Series: spyt
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## 0.0322
## s.e. 0.0331
##
## sigma^2 estimated as 1.934: log likelihood=-3082.54
## AIC=6169.08 AICc=6169.09 BIC=6180.03plot(forecast(arwd1,h=502),main="Random Walk ARIMA(0,1,0) With Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
arwd2 <- Arima(spyf,model=arwd1)$fitted
plot(spyt,main="Random Walk ARIMA(0,1,0) With Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(arwd2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="arwd2")
Differentiated First Order Autoregressive ARIMA(1,1,0) with constant
Multi-Steps Forecast
dar1 <- Arima(spyt,order=c(1,1,0),include.constant=T)
dar1## Series: spyt
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## -0.0877 0.0322
## s.e. 0.0237 0.0303
##
## sigma^2 estimated as 1.92: log likelihood=-3075.74
## AIC=6157.47 AICc=6157.48 BIC=6173.9plot(forecast(dar1,h=502),main="Differentiated First Order Autoregressive ARIMA(1,1,0) With Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
dar2 <- Arima(spyf,model=dar1)$fitted
plot(spyt,main="Differentiated First Order Autoregressive ARIMA(1,1,0) With Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(dar2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="dar2")
Brown Simple Exponential Smoothing ARIMA(0,1,1) without constant
Multi-Steps Forecast
abrown1 <- Arima(spyt,order=c(0,1,1))
abrown1## Series: spyt
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.0971
## s.e. 0.0249
##
## sigma^2 estimated as 1.919: log likelihood=-3075.51
## AIC=6155.03 AICc=6155.03 BIC=6165.98plot(forecast(abrown1,h=502),main="Brown Simple Exponential Smoothing ARIMA(0,1,1) Without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
abrown2 <- Arima(spyf,model=abrown1)$fitted
plot(spyt,main="Brown Simple Exponential Smoothing ARIMA(0,1,1) Without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(abrown2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="abrown2")
Simple Exponential Smoothing with Growth ARIMA(0,1,1) with constant
Multi-Steps Forecast
abrowng1 <- Arima(spyt,order=c(0,1,1),include.constant=T)
abrowng1## Series: spyt
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## -0.0979 0.0322
## s.e. 0.0249 0.0297
##
## sigma^2 estimated as 1.919: log likelihood=-3074.93
## AIC=6155.86 AICc=6155.87 BIC=6172.28plot(forecast(abrowng1,h=502),main="Simple Exponential Smoothing with Growth ARIMA(0,1,1) With Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
abrowng2 <- Arima(spyf,model=abrowng1)$fitted
plot(spyt,main="Simple Exponential Smoothing with Growth ARIMA(0,1,1) With Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(abrowng2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="abrowng2")
Holt Linear Trend ARIMA (0,2,2) without constant
Multi-Steps Forecast
aholt1 <- Arima(spyt,order=c(0,2,2))
aholt1## Series: spyt
## ARIMA(0,2,2)
##
## Coefficients:
## ma1 ma2
## -1.0961 0.0982
## s.e. 0.0249 0.0250
##
## sigma^2 estimated as 1.921: log likelihood=-3076.91
## AIC=6159.82 AICc=6159.83 BIC=6176.24plot(forecast(aholt1,h=502),main="Holt Linear Trend ARIMA(0,2,2) without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
aholt2 <- Arima(spyf,model=aholt1)$fitted
plot(spyt,main="Holt Linear Trend ARIMA(0,2,2) without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(aholt2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="aholt2")
Gardner Additive Damped Trend ARIMA(1,1,2) without constant
Multi-Steps Forecast
agardner1 <- Arima(spyt,order=c(1,1,2))
agardner1## Series: spyt
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.1346 0.0424 -0.0614
## s.e. 0.3998 0.3991 0.0406
##
## sigma^2 estimated as 1.916: log likelihood=-3073.42
## AIC=6154.84 AICc=6154.87 BIC=6176.74plot(forecast(agardner1,h=502),main="Gardner Additive Damped Trend ARIMA(1,1,2) without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
agardner2 <- Arima(spyf,model=agardner1)$fitted
plot(spyt,main="Gardner Additive Damped Trend ARIMA(1,1,2) without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(agardner2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="agardner2")
Seasonal Random Walk with Drift ARIMA(0,0,0)x(0,1,0)[m] with constant
Multi-Steps Forecast
srwd1 <- Arima(spyt,order=c(0,0,0),seasonal=c(0,1,0),include.constant=T)
srwd1## Series: spyt
## ARIMA(0,0,0)(0,1,0)[21] with drift
##
## Coefficients:
## drift
## 0.0309
## s.e. 0.0062
##
## sigma^2 estimated as 29.86: log likelihood=-5432.8
## AIC=10869.6 AICc=10869.61 BIC=10880.53plot(forecast(srwd1,h=502),main="Seasonal Random Walk with Drift ARIMA(0,0,0)x(0,1,0)[21] With Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
srwd2 <- Arima(spyf,model=srwd1)$fitted
plot(spyt,main="Seasonal Random Walk with Drift ARIMA(0,0,0)x(0,1,0)[21] With Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(srwd2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="srwd2")
Seasonal Random Trend ARIMA(0,1,0)x(0,1,0)[m] without constant
Multi-Steps Forecast
srt1 <- Arima(spyt,order=c(0,1,0),seasonal=c(0,1,0))
srt1## Series: spyt
## ARIMA(0,1,0)(0,1,0)[21]
##
## sigma^2 estimated as 4.123: log likelihood=-3705.6
## AIC=7413.19 AICc=7413.2 BIC=7418.66plot(forecast(srt1,h=502),main="Seasonal Random Trend ARIMA(0,1,0)x(0,1,0)[21] Without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
srt2 <- Arima(spyf,model=srt1)$fitted
plot(spyt,main="Seasonal Random Trend ARIMA(0,1,0)x(0,1,0)[21] Without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(srt2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="srt2")
General Seasonal Model ARIMA(0,1,1)x(0,1,1)[m] without constant
Multi-Steps Forecast
gseas1 <- Arima(spyt,order=c(0,1,1),seasonal=c(0,1,1))
gseas1## Series: spyt
## ARIMA(0,1,1)(0,1,1)[21]
##
## Coefficients:
## ma1 sma1
## -0.098 -1.000
## s.e. 0.025 0.012
##
## sigma^2 estimated as 1.926: log likelihood=-3087.93
## AIC=6181.85 AICc=6181.87 BIC=6198.24plot(forecast(gseas1,h=502),main="General Seasonal Model ARIMA(0,1,1)x(0,1,1)[21] Without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
gseas2 <- Arima(spyf,model=gseas1)$fitted
plot(spyt,main="General Seasonal Model ARIMA(0,1,1)x(0,1,1)[21] Without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(gseas2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="gseas2")
General First Order Autoregressive Seasonal Model ARIMA(1,0,1)x(0,1,1)[m] with constant
Multi-Steps Forecast
gsar1 <- Arima(spyt,order=c(1,0,1),seasonal=c(0,1,1),include.constant=T)
gsar1## Series: spyt
## ARIMA(1,0,1)(0,1,1)[21] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## 0.9993 -0.0984 -1.0000 0.0317
## s.e. 0.0014 0.0251 0.0097 0.0226
##
## sigma^2 estimated as 1.926: log likelihood=-3088.25
## AIC=6186.49 AICc=6186.53 BIC=6213.81plot(forecast(gsar1,h=502),main="General First Order Autoregressive Seasonal Model ARIMA(1,0,1)x(0,1,1)[21] With Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
gsar2 <- Arima(spyf,model=gsar1)$fitted
plot(spyt,main="General First Order Autoregressive Seasonal Model ARIMA(1,0,1)x(0,1,1)[21] With Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(gsar2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="gsar2")
Seasonally Differentiated First Order Autoregressive ARIMA(1,0,0)x(0,1,0)[m] with constant
Multi-Steps Forecast
sdar1 <- Arima(spyt,order=c(1,0,0),seasonal=c(0,1,0),include.constant=T)
sdar1## Series: spyt
## ARIMA(1,0,0)(0,1,0)[21] with drift
##
## Coefficients:
## ar1 drift
## 0.9305 0.0309
## s.e. 0.0087 0.0325
##
## sigma^2 estimated as 3.984: log likelihood=-3677.73
## AIC=7361.47 AICc=7361.48 BIC=7377.86plot(forecast(sdar1,h=502),main="Seasonally Differentiated First Order Autoregressive ARIMA(1,0,0)x(0,1,0)[21] With Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
sdar2 <- Arima(spyf,model=sdar1)$fitted
plot(spyt,main="Seasonally Differentiated First Order Autoregressive ARIMA(1,0,0)x(0,1,0)[21] With Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(sdar2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="sdar2")
Holt-Winters Additive Seasonality ARIMA (0,1,m+1)X(0,1,0)[m] without constant
Multi-Steps Forecast
ahwa1 <- Arima(spyt,order=c(0,1,22),seasonal=c(0,1,0))
ahwa1## Series: spyt
## ARIMA(0,1,22)(0,1,0)[21]
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## -0.0976 0.0030 0.0018 0.0066 -0.0038 0.0041 0.0029 0.0056
## s.e. 0.0259 0.0097 0.0090 0.0098 0.0089 0.0096 0.0091 0.0091
## ma9 ma10 ma11 ma12 ma13 ma14 ma15 ma16
## -0.0034 0.0060 -0.0014 0.0070 -0.0010 0.0002 0.0023 0.0059
## s.e. 0.0087 0.0091 0.0087 0.0087 0.0088 0.0091 0.0092 0.0088
## ma17 ma18 ma19 ma20 ma21 ma22
## -0.0006 0.0005 0.0046 0.0087 -0.9972 0.0904
## s.e. 0.0096 0.0088 0.0096 0.0091 0.0106 0.0258
##
## sigma^2 estimated as 1.946: log likelihood=-3084.96
## AIC=6215.93 AICc=6216.57 BIC=6341.57plot(forecast(ahwa1,h=502),main="Holt-Winters Additive Seasonality ARIMA (0,1,22)X(0,1,0)[21] Without Constant Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
ahwa2 <- Arima(spyf,model=ahwa1)$fitted
plot(spyt,main="Holt-Winters Additive Seasonality ARIMA (0,1,22)X(0,1,0)[21] Without Constant Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(ahwa2,col="green")
lines(spyf,lty=3)
legend("bottomright",col="green",lty=1,legend="ahwa2")
Model Selection
Multi-Steps Forecast
abest1 <- auto.arima(spyt,max.d=1)
abest1## Series: spyt
## ARIMA(1,1,1)(0,0,1)[21]
##
## Coefficients:
## ar1 ma1 sma1
## 0.4548 -0.5443 -0.0538
## s.e. 0.1829 0.1724 0.0243
##
## sigma^2 estimated as 1.911: log likelihood=-3071.14
## AIC=6150.28 AICc=6150.3 BIC=6172.18plot(forecast(abest1,h=502),main="Best ARIMA (1,1,1)X(0,0,1)[21] Model 1",ylab="Price",xlab="Month")
lines(spy,lty=3)
One-Step Forecast without Re-Estimation
abest2 <- Arima(spyf,model=abest1)$fitted
plot(spyt,main="Best ARIMA (1,1,1)X(0,0,1)[21] Model 2",ylab="Price",xlab="Month",xlim=c(1,108),ylim=c(min(spy),max(spy)))
lines(abest2,col="green")
lines(spyf,lty=3)
legend("bottom",col="green",lty=1,legend="abest2")
Model Forecasting Accuracy
Multi-Steps Forecast
accuracy(forecast(arw1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.092649 21.806753 19.338467 9.783360 9.932039 4.468997accuracy(forecast(arwd1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 10.999963 13.119287 11.368150 5.618607 5.837789 2.627107accuracy(forecast(dar1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 10.999972 13.120168 11.368527 5.618541 5.837944 2.627194accuracy(forecast(abrown1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.075109 21.791398 19.322185 9.774079 9.923510 4.465234accuracy(forecast(abrowng1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 10.978235 13.100997 11.348987 5.607115 5.827804 2.622679accuracy(forecast(aholt1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## -4.507895 8.833436 5.746543 -2.362641 3.010410 1.327990accuracy(forecast(agardner1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.001951 21.727394 19.254289 9.735365 9.887947 4.449544accuracy(forecast(srwd1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 13.038919 15.049819 13.328243 6.694731 6.867039 3.080072accuracy(forecast(srt1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## -26.604272 32.533819 26.646404 -13.730014 13.754567 6.157815accuracy(forecast(gseas1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 11.148177 13.245968 11.496770 5.696880 5.904604 2.656830accuracy(forecast(gsar1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 13.120317 15.285071 13.400982 6.712831 6.882116 3.096882accuracy(forecast(sdar1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 11.017799 13.259949 11.480642 5.625770 5.899679 2.653103accuracy(forecast(ahwa1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## -4.008874 8.448890 5.511739 -2.104960 2.889244 1.273728accuracy(forecast(abest1,h=502),spyf)[2,1:6]## ME RMSE MAE MPE MAPE MASE
## 19.201945 21.906451 19.440915 9.840941 9.985498 4.492672One-Step Forecast without Re-Estimation
accuracy(arw2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.05564223 1.60937675 1.17225198 0.02670005 0.61988109accuracy(arwd2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.023415512 1.608581909 1.169821638 0.009643717 0.618550730accuracy(dar2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.02581641 1.62061446 1.17310645 0.01063506 0.62050134accuracy(abrown2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.06203304 1.62277740 1.17626341 0.02974240 0.62224018accuracy(abrowng2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.02638518 1.62197367 1.17394464 0.01087345 0.62098689accuracy(aholt2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## -0.02523900 1.62374731 1.17461137 -0.01602416 0.62138787accuracy(agardner2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.06505048 1.61933361 1.17464928 0.03120562 0.62147586accuracy(srwd2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.6182851 5.6510653 4.0818960 0.3032664 2.1517638accuracy(srt2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.0013631413 2.1408497959 1.5879475256 0.0005683963 0.8339597397accuracy(gseas2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## -0.006667440 1.589746708 1.142638461 -0.004962661 0.601294346accuracy(gsar2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.04360977 1.58929200 1.14820903 0.02092702 0.60445442accuracy(sdar2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.04539233 2.10474248 1.56159623 0.02224295 0.82096440accuracy(ahwa2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## -0.004890197 1.595816380 1.150836381 -0.004040726 0.605539032accuracy(abest2,spyf)[1,1:5]## ME RMSE MAE MPE MAPE
## 0.07172688 1.62335603 1.17936377 0.03456213 0.62393556Residuals White Noise
Best fitting model
abest1## Series: spyt
## ARIMA(1,1,1)(0,0,1)[21]
##
## Coefficients:
## ar1 ma1 sma1
## 0.4548 -0.5443 -0.0538
## s.e. 0.1829 0.1724 0.0243
##
## sigma^2 estimated as 1.911: log likelihood=-3071.14
## AIC=6150.28 AICc=6150.3 BIC=6172.18Normal and Partial Autocorrelation Functions ACF & PACF
acf(residuals(abest1))
pacf(residuals(abest1))
Ljung.Box Autocorrelation Test
Box.test(residuals(abest1),lag=21,type="Ljung")##
## Box-Ljung test
##
## data: residuals(abest1)
## X-squared = 33.773, df = 21, p-value = 0.03832