Lab sessions: Day 3
Rob J Hyndman
28 September 2017
Lab Session 13
retaildata <- readxl::read_excel("retail.xlsx", skip = 1)
mytimeseries <- ts(retaildata[["A3349873A"]],
frequency=12, start=c(1982,4))
(etsmod <- ets(mytimeseries))## ETS(M,A,M)
##
## Call:
## ets(y = mytimeseries)
##
## Smoothing parameters:
## alpha = 0.5067
## beta = 1e-04
## gamma = 0.3222
##
## Initial states:
## l = 63.0202
## b = 0.7934
## s=0.9391 0.912 0.9294 1.5281 1.0577 0.9868
## 0.9604 0.941 0.9431 0.901 0.9661 0.9353
##
## sigma: 0.0488
##
## AIC AICc BIC
## 4018.808 4020.494 4085.835myfn <- function(y,h)
{
y %>% ets(model="MAM", damped=FALSE) %>% forecast(h=h)
}
e <- tsCV(mytimeseries, myfn)
sqrt(mean(e^2, na.rm=TRUE))## [1] 15.63519ggtsdisplay(e)## Warning: Removed 1 rows containing missing values (geom_point).
myfn <- function(y,h)
{
y %>% ets %>% forecast(h=h)
}
e2 <- tsCV(mytimeseries, myfn)
sqrt(mean(e2^2, na.rm=TRUE))## [1] 16.17141ggtsdisplay(e2)## Warning: Removed 1 rows containing missing values (geom_point).
Lab Session 14
lambda <- BoxCox.lambda(mytimeseries)
fit1 <- auto.arima(mytimeseries,
xreg=fourier(mytimeseries, 1),
seasonal=FALSE, lambda=lambda)
fit2 <- auto.arima(mytimeseries,
xreg=fourier(mytimeseries, 2),
seasonal=FALSE, lambda=lambda)
fit3 <- auto.arima(mytimeseries,
xreg=fourier(mytimeseries, 3),
seasonal=FALSE, lambda=lambda)
fit4 <- auto.arima(mytimeseries,
xreg=fourier(mytimeseries, 4),
seasonal=FALSE, lambda=lambda)
fit5 <- auto.arima(mytimeseries,
xreg=fourier(mytimeseries, 5),
seasonal=FALSE, lambda=lambda)
fit6 <- auto.arima(mytimeseries,
xreg=fourier(mytimeseries, 6),
seasonal=FALSE, lambda=lambda)
fit1$aicc## [1] -29.77232fit2$aicc## [1] -157.1155fit3$aicc## [1] -303.7034fit4$aicc## [1] -551.1655fit5$aicc## [1] -571.5915fit6$aicc## [1] -639.1607checkresiduals(fit6)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(0,1,1) errors
## Q* = 55.724, df = 11, p-value = 5.709e-08
##
## Model df: 13. Total lags used: 24fc <- forecast(fit6,
xreg=fourier(mytimeseries, 6,36))
autoplot(fc)
newretaildata <- readxl::read_excel("8501011.xls", sheet=2, skip=9)
mynewdata <- ts(newretaildata[["A3349873A"]],
frequency=12, start=c(1982,4))
mynewdata <- window(mynewdata, start=2014)
autoplot(fc) +
autolayer(mynewdata, series="Future")
Lab Session 15
library(hts)
plot(infantgts)
smatrix(infantgts)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 1 1 1 1 1 1 1 1 1 1 1 1
## [2,] 1 1 1 1 1 1 1 1 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 1 1 1 1 1
## [4,] 1 0 0 0 0 0 0 0 1 0 0 0 0
## [5,] 0 1 0 0 0 0 0 0 0 1 0 0 0
## [6,] 0 0 1 0 0 0 0 0 0 0 1 0 0
## [7,] 0 0 0 1 0 0 0 0 0 0 0 1 0
## [8,] 0 0 0 0 1 0 0 0 0 0 0 0 1
## [9,] 0 0 0 0 0 1 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 1 0 0 0 0 0 0
## [11,] 0 0 0 0 0 0 0 1 0 0 0 0 0
## [12,] 1 0 0 0 0 0 0 0 0 0 0 0 0
## [13,] 0 1 0 0 0 0 0 0 0 0 0 0 0
## [14,] 0 0 1 0 0 0 0 0 0 0 0 0 0
## [15,] 0 0 0 1 0 0 0 0 0 0 0 0 0
## [16,] 0 0 0 0 1 0 0 0 0 0 0 0 0
## [17,] 0 0 0 0 0 1 0 0 0 0 0 0 0
## [18,] 0 0 0 0 0 0 1 0 0 0 0 0 0
## [19,] 0 0 0 0 0 0 0 1 0 0 0 0 0
## [20,] 0 0 0 0 0 0 0 0 1 0 0 0 0
## [21,] 0 0 0 0 0 0 0 0 0 1 0 0 0
## [22,] 0 0 0 0 0 0 0 0 0 0 1 0 0
## [23,] 0 0 0 0 0 0 0 0 0 0 0 1 0
## [24,] 0 0 0 0 0 0 0 0 0 0 0 0 1
## [25,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [26,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [27,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16]
## [1,] 1 1 1
## [2,] 0 0 0
## [3,] 1 1 1
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
## [7,] 0 0 0
## [8,] 0 0 0
## [9,] 1 0 0
## [10,] 0 1 0
## [11,] 0 0 1
## [12,] 0 0 0
## [13,] 0 0 0
## [14,] 0 0 0
## [15,] 0 0 0
## [16,] 0 0 0
## [17,] 0 0 0
## [18,] 0 0 0
## [19,] 0 0 0
## [20,] 0 0 0
## [21,] 0 0 0
## [22,] 0 0 0
## [23,] 0 0 0
## [24,] 0 0 0
## [25,] 1 0 0
## [26,] 0 1 0
## [27,] 0 0 1# Forecast 10-step-ahead and reconcile the forecasts
infantforecast <- forecast(infantgts, h=10)
# Plot the forecasts including the last ten historical years
plot(infantforecast, include=10)
# Create a matrix of all aggregated time series
allts_infant <- aggts(infantgts)
# Forecast all series using ARIMA models
allf <- matrix(, nrow=10, ncol=ncol(allts_infant))
for(i in 1:ncol(allts_infant))
allf[,i] <- forecast(auto.arima(allts_infant[,i]), h=10)$mean
allf <- ts(allf, start=2004)
# combine the forecasts with the group matrix to get a gts object
y.f <- combinef(allf, groups = infantgts$groups)# set up training and testing samples
data <- window(infantgts, end=1993)
test <- window(infantgts, start=1994)
# Compute forecasts on training data
forecast <- forecast(data, h=10)
# calculate ME, RMSE, MAE, MAPE, MPE and MASE
accuracy(forecast, test)## Total Sex/female Sex/male State/NSW State/VIC State/QLD
## ME -176.6273007 -94.236398 -82.3909023 -18.9655562 -47.8419698 -54.76347
## RMSE 187.5697408 101.665490 90.2361037 43.9579759 52.2161348 59.56858
## MAE 176.6273007 94.236398 82.3909023 31.0069140 47.8419698 54.76347
## MAPE 13.4073094 16.397044 11.1439738 7.4636054 15.8447985 20.04418
## MPE -13.4073094 -16.397044 -11.1439738 -4.8394048 -15.8447985 -20.04418
## MASE 0.9161167 1.009676 0.7554178 0.3478711 0.6513543 1.18152
## State/SA State/WA State/NT State/ACT State/TAS NSW female
## ME -26.1252165 -31.325056 -1.154467 -5.8985681 9.4470017 18.2371228
## RMSE 30.3308774 35.870321 3.276476 8.7802047 14.4095410 27.8699565
## MAE 27.2237241 32.179220 1.969059 6.4454241 11.6983720 22.5470092
## MAPE 34.5193106 28.757408 5.446349 39.6241455 26.1601892 11.7054120
## MPE -33.5804152 -28.206335 -3.504192 -37.9152203 19.0087895 9.0774325
## MASE 0.8867663 1.112826 0.159438 0.9967151 0.5153468 0.5179252
## VIC female QLD female SA female WA female NT female ACT female
## ME -30.5623558 -39.747346 -16.9049390 -13.7193250 -14.598216 1.764803
## RMSE 35.9775320 43.707182 18.1595354 17.3057700 16.377837 4.862153
## MAE 30.5623558 39.747346 16.9049390 15.6654718 14.598216 4.236210
## MAPE 24.0753097 35.001144 49.3217301 32.6891019 100.193203 52.773149
## MPE -24.0753097 -35.001144 -49.3217301 -29.8685992 -100.193203 -10.653638
## MASE 0.8477769 1.361988 0.9595992 0.8960232 2.051272 1.008621
## TAS female NSW male VIC male QLD male SA male WA male
## ME 1.2938582 -37.2026790 -17.2796140 -15.0161231 -9.2202775 -17.605731
## RMSE 4.7119423 44.2538294 21.9312363 19.6645834 14.4344083 21.320034
## MAE 4.2501755 37.2026790 18.5483228 15.6371538 11.8662186 19.660420
## MAPE 21.2464883 16.3662946 10.9862151 10.3865131 29.1346925 31.873229
## MPE 1.9326002 -16.3662946 -10.3376522 -10.0512901 -25.6531910 -29.797786
## MASE 0.3744648 0.6675122 0.4244467 0.6488446 0.6207263 1.203699
## NT male ACT male TAS male
## ME 13.443750 -7.663371 8.1531435
## RMSE 15.687685 8.453458 10.5715649
## MAE 14.296623 7.663371 8.6168884
## MAPE 61.481614 73.895074 38.8128901
## MPE 54.921048 -73.895074 36.0849789
## MASE 2.071974 1.522524 0.7591972# Overall forecast accuracy
rmse <- accuracy(forecast, test)["RMSE",]
sqrt(sum(rmse^2))## [1] 270.5795forecast_bu <- forecast(data, h=10, method="bu")
sqrt(sum(accuracy(forecast_bu, test)["RMSE",]^2))## [1] 362.0634Lab Session 16
library(thief)
lambda <- BoxCox.lambda(mytimeseries)
fcarima <- forecast(auto.arima(mytimeseries, lambda=lambda), h=36)
fcthief <- thief(mytimeseries, usemodel='arima', lambda=lambda, h=36)
autoplot(mytimeseries) +
autolayer(fcthief, series="Thief") +
autolayer(fcarima, series="ARIMA", PI=FALSE) +
autolayer(mynewdata, series="New data")
autoplot(window(mytimeseries, start=2012)) +
autolayer(fcthief, series="Thief") +
autolayer(fcarima, series="ARIMA", PI=FALSE) +
autolayer(mynewdata, series="New data")
Lab Session 17
fit <- tbats(mytimeseries)
checkresiduals(fit)
##
## Ljung-Box test
##
## data: Residuals from TBATS
## Q* = 71.339, df = 3, p-value = 2.22e-15
##
## Model df: 27. Total lags used: 30forecast(fit) %>% autoplot
fit <- tbats(gasoline)
checkresiduals(fit)
##
## Ljung-Box test
##
## data: Residuals from TBATS
## Q* = 160.5, df = 82.357, p-value = 5.739e-07
##
## Model df: 22. Total lags used: 104.357142857143forecast(fit) %>% autoplot
Lab Session 18
fit <- nnetar(mytimeseries)
checkresiduals(fit)
forecast(fit, h=36, PI=TRUE) %>% autoplot
Lab Session 19
lambda <- BoxCox.lambda(mytimeseries)
e1 <- tsCV(mytimeseries, snaive, h=12)
e2 <- tsCV(mytimeseries, hw, h=12,
seasonal='mult')
e3 <- tsCV(mytimeseries, thief, h=12)
e4 <- tsCV(mytimeseries, stlf, lambda=lambda, h=12)f0 <- function(y,h, ...) {
forecast(ets(y, ...), h=h)
}
f1 <- function(y,h, ...) {
forecast(auto.arima(y, ...), h=h)
}
f2 <- function(y,h, ...) {
forecast(tbats(y, ...), h=h)
}
f3 <- function(y,h, ...) {
forecast(nnetar(y, ...), h=h)
}
e5 <- tsCV(mytimeseries, f0, h=12)
e6 <- tsCV(mytimeseries, f1, lambda=lambda, h=12)
e7 <- tsCV(mytimeseries, f2, h=12)
e8 <- tsCV(mytimeseries, f3, lambda=lambda, h=12)c(
SNAIVE=mean(e1^2, na.rm=TRUE),
HW=mean(e2^2, na.rm=TRUE),
THIEF=mean(e3^2, na.rm=TRUE),
STLF=mean(e4^2, na.rm=TRUE),
ETS=mean(e5^2, na.rm=TRUE),
ARIMA=mean(e6^2, na.rm=TRUE),
TBATS=mean(e7^2, na.rm=TRUE),
NNETAR=mean(e8^2, na.rm=TRUE)
)## SNAIVE HW THIEF STLF ETS ARIMA TBATS
## 676.6104 639.7708 652.2999 652.8351 628.3876 629.0767 1138.1412
## NNETAR
## 3446.4710ETS does the best on this example. TBATS and NNETAR are not even competitive.