9.7 Exercises #2
#Loading required packages
setwd("C:/Users/Public")
library(fpp2)
## Warning: package 'fpp2' was built under R version 3.5.1
## Loading required package: ggplot2
## Loading required package: forecast
## Warning: package 'forecast' was built under R version 3.5.1
## Loading required package: fma
## Warning: package 'fma' was built under R version 3.5.1
## Loading required package: expsmooth
## Warning: package 'expsmooth' was built under R version 3.5.1
#Input data
data(LakeHuron)
plot(LakeHuron)
#Piecewise Linear Model
huron.lm<-tslm(LakeHuron~trend)
summary(huron.lm)
##
## Call:
## tslm(formula = LakeHuron ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.50997 -0.72726 0.00083 0.74402 2.53565
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 580.202037 0.230111 2521.398 < 2e-16 ***
## trend -0.024201 0.004036 -5.996 3.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.13 on 96 degrees of freedom
## Multiple R-squared: 0.2725, Adjusted R-squared: 0.2649
## F-statistic: 35.95 on 1 and 96 DF, p-value: 3.545e-08
t <- time(huron)
tb1 <- ts(pmax(0, t - 1920), start = 1875)
huron.pw <- tslm(huron ~ t + tb1)
summary(huron.pw)
##
## Call:
## tslm(formula = huron ~ t + tb1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5571 -0.6261 -0.1377 0.8337 2.3938
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 124.42259 17.19863 7.234 1.19e-10 ***
## t -0.06048 0.00904 -6.690 1.54e-09 ***
## tb1 0.06555 0.01490 4.398 2.85e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.036 on 95 degrees of freedom
## Multiple R-squared: 0.3956, Adjusted R-squared: 0.3828
## F-statistic: 31.08 on 2 and 95 DF, p-value: 4.117e-11
#Forecast for next 30 years
fc<-forecast(huron.lm, h=30)
plot(fc)
9.7 Exercises #6
library(ggplot2)
library(fpp2)
#Data
retail <- read.csv("retail.csv")
mytimeseries <- ts(retail[,4], frequency=12, start=c(1982,4))
train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)
lambda <- 0
# Find Fourier terms with min AIC
bestmodel <- list(aicc=Inf)
for(k in seq(6)) {
fit <- auto.arima(train, lambda=lambda,
xreg=fourier(train, K=k), seasonal=FALSE)
if(fit$aicc < bestmodel$aicc) {
bestmodel <- fit
bestK <- k
}
}
fc <- forecast(bestmodel,
xreg=fourier(train, bestK, h=length(test)))
#Check residuals
checkresiduals(fc)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(0,1,5) errors
## Q* = 76.258, df = 17, p-value = 1.757e-09
##
## Model df: 7. Total lags used: 24
#ggplot
autoplot(fc)
autoplot(fc) + autolayer(test)
accuracy(fc, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 6.93591 55.71739 29.34976 -72.369765 92.39859 0.748925
## Test set 21.61329 50.60509 40.84163 6.497492 22.08678 1.042166
## ACF1 Theil's U
## Training set 0.09311753 NA
## Test set 0.53445609 1.279399
10.8 Exercises #4
library(fpp2)
library(hts)
## Warning: package 'hts' was built under R version 3.5.1
library(ggplot2)
library(forecast)
head(visnights)
## NSWMetro NSWNthCo NSWSthCo NSWSthIn NSWNthIn QLDMetro QLDCntrl
## 1998 Q1 9.047095 8.565678 5.818029 2.679538 2.977507 12.106052 2.748374
## 1998 Q2 6.962126 7.124468 2.466437 3.010732 3.477703 7.786687 4.040915
## 1998 Q3 6.871963 4.716893 1.928053 3.328869 3.014770 11.380024 5.343964
## 1998 Q4 7.147293 6.269299 2.797556 2.417772 3.757972 9.311460 4.260419
## 1999 Q1 7.956923 9.493901 4.853681 3.224285 3.790760 12.671942 4.186113
## 1999 Q2 6.542243 5.401201 2.759843 2.428489 3.395284 9.582965 4.237806
## QLDNthCo SAUMetro SAUCoast SAUInner VICMetro VICWstCo VICEstCo
## 1998 Q1 2.137234 2.881372 2.591997 0.8948773 7.490382 2.4420048 3.381972
## 1998 Q2 2.269596 2.124736 1.375780 0.9792509 5.198178 0.9605047 1.827940
## 1998 Q3 4.890227 2.284870 1.079542 0.9803289 5.244217 0.7559744 1.351952
## 1998 Q4 2.621548 1.785889 1.497664 1.5094343 6.274246 1.2716040 1.493415
## 1999 Q1 2.483203 2.293873 2.247684 0.9635227 9.187422 2.3850583 2.896929
## 1999 Q2 3.377830 2.197418 1.672802 0.9968803 4.992303 1.3288638 1.547901
## VICInner WAUMetro WAUCoast WAUInner OTHMetro OTHNoMet
## 1998 Q1 5.326655 3.075779 3.066555 0.6949954 3.437924 2.073469
## 1998 Q2 4.441119 2.154929 3.334405 0.5576796 2.677081 1.787939
## 1998 Q3 3.815645 2.787286 4.365844 1.0061844 3.793743 2.345021
## 1998 Q4 3.859567 2.752910 4.521996 1.1725514 3.304231 1.943689
## 1999 Q1 4.588755 3.519564 3.579347 0.3981829 3.510819 2.165838
## 1999 Q2 4.070401 3.160430 3.408533 0.5960182 2.871867 1.803940
autoplot(visnights)
tourism.hts <- hts(visnights, characters = c(3,5))
#Compare with Question 2
fc_visnights_arima <- forecast(
tourism.hts, h = 8, method = "bu", fmethod = "arima")
# Plot
plot(fc_visnights_arima, levels = 0, color_lab = TRUE)
title(main = "Total visnights")
plot(fc_visnights_arima, levels = 1, color_lab = TRUE)
title(main = "Grouped by State")
plot(fc_visnights_arima, levels = 2, color_lab = TRUE)
title(main = "Grouped by zone in each state")
#Question 4
# Forecast 8-steps-ahead
visforecast <- forecast(visnights, h=8)
plot(visforecast, include=8)
# aggregate visnights time series matrix.
visnights.total <- rowSums(visnights)
visnights.total.ts <- ts(visnights.total,
start = 1998,
frequency = 4)
str(visnights.total.ts)
## Time-Series [1:76] from 1998 to 2017: 83.4 64.6 71.3 70 86.4 ...
# Auto arima
visnights.autoarima <- auto.arima(visnights.total.ts)
autoplot(visnights.total.ts) +
autolayer(visnights.autoarima$fitted)
# the model has similar trend
fc_visnights_total_autoarima <- forecast(
visnights.autoarima, h = 8
)
fc_visnights_total_autoarima$model
## Series: visnights.total.ts
## ARIMA(0,1,1)(0,1,1)[4]
##
## Coefficients:
## ma1 sma1
## -0.5608 -0.8616
## s.e. 0.0991 0.1401
##
## sigma^2 estimated as 11.11: log likelihood=-188.2
## AIC=382.39 AICc=382.75 BIC=389.18
#based on the result, arima(0, 1, 1)(0, 1, 1)is the best fit.
# get ts object of bottom-up forecasts.
visnights.bu.ts <- ts(
rowSums(fc_visnights_arima$bts),
start = 2017,
frequency = 4)
autoplot(fc_visnights_total_autoarima) +
autolayer(visnights.bu.ts)
#the plot shows the forecast
10.8 Exercises #5
hts.train <- window(tourism.hts, end=c(2014,4))
hts.test <- window(tourism.hts, start=2015)
forecast.visnights.ets = forecast(
hts.train, h = 8,
method = "comb", weights = "wls", fmethod = "ets"
)
forecast.visnights.ets = forecast(
hts.train, h = 8,
method = "bu", fmethod = "ets"
)
#compare accuracy based on mse
vis <- matrix(NA, ncol = 4, nrow = 4)
rownames(vis) <- c("Total", "State", "Bottom", "All series")
colnames(vis) <- c("Bottom-up MAPE", "Bottom-up MASE", "Optimal MAPE", "Optimal MASE")
vis[1,] <- c(
accuracy.gts(
forecast.visnights.ets,
hts.test,
levels = 0
)[c("MAPE","MASE"),"Total"],
accuracy.gts(
forecast.visnights.ets,
hts.test,
levels = 0
)[c("MAPE","MASE"),"Total"]
)
j=2
for(i in c(1:2)){
vis[j,] <- c(
mean(accuracy.gts(
forecast.visnights.ets,
hts.test,
levels = i)["MAPE",]),
mean(accuracy.gts(
forecast.visnights.ets,
hts.test,
levels = i)["MASE",]),
mean(accuracy.gts(
forecast.visnights.ets,
hts.test,
levels = i)["MAPE",]),
mean(accuracy.gts(
forecast.visnights.ets,
hts.test,
levels = i)["MASE",])
)
j=j+1
}
vis[4,] <- c(
mean(accuracy.gts(
forecast.visnights.ets,
hts.test
)["MAPE",]),
mean(accuracy.gts(
forecast.visnights.ets,
hts.test
)["MASE",]),
mean(accuracy.gts(
forecast.visnights.ets,
hts.test
)["MAPE",]),
mean(accuracy.gts(
forecast.visnights.ets,
hts.test
)["MASE",])
)
library(knitr)
## Warning: package 'knitr' was built under R version 3.5.1
knitr::kable(vis, digits=2, booktabs=TRUE)
| Total |
8.59 |
2.14 |
8.59 |
2.14 |
| State |
10.48 |
1.44 |
10.48 |
1.44 |
| Bottom |
12.57 |
1.15 |
12.57 |
1.15 |
| All series |
11.95 |
1.25 |
11.95 |
1.25 |
#the optimally reconciled forecasts are better than the other two based on the mse