Hyndman Chapters 7 & 8 Exercises
Sam Neil
7/28/2021
library(fpp2)## Warning: package 'fpp2' was built under R version 4.0.5## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## -- Attaching packages ---------------------------------------------- fpp2 2.4 --## v ggplot2 3.3.3 v fma 2.4
## v forecast 8.15 v expsmooth 2.3## Warning: package 'forecast' was built under R version 4.0.5## Warning: package 'fma' was built under R version 4.0.5## Warning: package 'expsmooth' was built under R version 4.0.5## library(forecast)
library(seasonal)## Warning: package 'seasonal' was built under R version 4.0.5library(psych)## Warning: package 'psych' was built under R version 4.0.5##
## Attaching package: 'psych'## The following object is masked from 'package:seasonal':
##
## outlier## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.0.5##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)###CHAPTER 7###
#Exercise 1# Consider the pigs series — the number of pigs slaughtered in Victoria each month.
- Use the ses() function in R to find the optimal values of α and ℓ0, and generate forecasts for the next four months.
pigses = ses(pigs, h = 4)
summary(pigses)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = pigs, h = 4)
##
## Smoothing parameters:
## alpha = 0.2971
##
## Initial states:
## l = 77260.0561
##
## sigma: 10308.58
##
## AIC AICc BIC
## 4462.955 4463.086 4472.665
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249 0.01282239
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995 98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995 98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995 98816.41 83958.37 113674.4 76092.99 121539.8- Compute a 95% prediction interval for the first forecast using ^y±1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
sd = sd(pigses$residuals)
pigPF = pigses$mean[1]
upper95 = pigPF + 1.96*sd
lower95 = pigPF - 1.96*sdupper95## [1] 118952.8lower95## [1] 78679.97intervals from R
pigses$upper[1, "95%"]## 95%
## 119020.8pigses$lower[1, "95%"]## 95%
## 78611.97The intervals are not exactly the same, but they are very similar. Both the upper 95% and lower 95% is off by 68 when comparing the self-calculated values to those generated by R.
Exercise 2 Write your own function to implement simple exponential smoothing. The function should take arguments y (the time series), alpha (the smoothing parameter α) and level (the initial level ℓ0).It should return the forecast of the next observation in the series. Does it give the same forecast as ses()?
ses_homebrew = function(y, alpha, level){
yhat = level
for (i in 1:length(y)){
yhat = alpha*y[i] + (1-alpha)*yhat
}
return(yhat)
}alpha = pigses$model$par[1]
level = pigses$model$par[2]
ses_homebrew(pigs, alpha = alpha, level=level)## alpha
## 98816.41This function yields identical results to the regular ses() function, which also predicted 98816.41.
Exercise 3 Modify your function from the previous exercise to return the sum of squared errors rather than the forecast of the next observation. Then use the optim() function to find the optimal values of α and ℓ0. Do you get the same values as the ses() function?
ses_homebrew2 = function(y, pars=c(alpha, level)){
yhat = level
error = 0
SSE = 0
alpha = pars[1]
level = pars[2]
for (i in 1:length(y)){
yhat = alpha*y[i] + (1-alpha)*yhat
error = y[i]-yhat
SSE = SSE + error^2
}
return(SSE)
}optimizedPig = optim(par = c(alpha, pigs[1]), y = pigs, fn = ses_homebrew2)
optimizedPig$par[1]## alpha
## 0.999957optimizedPig$par[2]##
## 76381.26pigses$model$par[1]## alpha
## 0.2971488pigses$model$par[2]## l
## 77260.06The alpha values I got from my function were significantly different from those in the actual model. I tried several different methods, but they all yielded a value of approximately 1. The variation between these different function attempts all kept a level around 76,000-77,000 though, which is within a reasonable range of the actual value.
Exercise 4 Combine your previous two functions to produce a function which both finds the optimal values of
α and ℓ0, and produces a forecast of the next observation in the series.
ses_homebrew3 = function(y, pars=c(alpha, level)){
fc = 0
SSE = function(y, pars=c(alpha, level)){
yhat = level
error = 0
SSE = 0
alpha = pars[1]
level = pars[2]
for (i in 1:length(y)){
error = y[i]-yhat
SSE = SSE + error^2
yhat = alpha*y[i] + (1-alpha)*yhat
}
fc <<- yhat
return(SSE)
}
optimizedPig = optim(par = c(alpha, pigs[1]), y = pigs, fn = SSE)
return(fc)
}
ses_homebrew3(pigs, c(alpha, pigs[1]))## alpha
## 98816.29pigses$mean[1]## [1] 98816.41This function took a very long time to refine, but produced incredible results. The function is only 0.05 off from the actual value. Its initial version was off by a factor of 10^26.
Exercise 5 Data set books contains the daily sales of paperback and hardcover books at the same store. The task is to forecast the next four days’ sales for paperback and hardcover books.
a.Plot the series and discuss the main features of the data.
autoplot(books)
summary(books)## Paperback Hardcover
## Min. :111.0 Min. :128.0
## 1st Qu.:167.2 1st Qu.:170.5
## Median :189.0 Median :200.5
## Mean :186.4 Mean :198.8
## 3rd Qu.:207.2 3rd Qu.:222.0
## Max. :247.0 Max. :283.0This plot covers a time series for both paperback and hardcover books. While the plots are different, both appear to have an upward rend and some form of seasonality to them. Hardcover initially has a smaller magnitude than paperback, but has a larger one by the end of the time series.
- Use the ses() function to forecast each series, and plot the forecasts.
pb <- books[,1]
hc <- books[,2]
sesPB <- ses(pb)
sesHC <- ses(hc)summary(sesPB)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = pb)
##
## Smoothing parameters:
## alpha = 0.1685
##
## Initial states:
## l = 170.8271
##
## sigma: 34.8183
##
## AIC AICc BIC
## 318.9747 319.8978 323.1783
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303 -0.2117522
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 207.1097 162.4882 251.7311 138.8670 275.3523
## 32 207.1097 161.8589 252.3604 137.9046 276.3147
## 33 207.1097 161.2382 252.9811 136.9554 277.2639
## 34 207.1097 160.6259 253.5935 136.0188 278.2005
## 35 207.1097 160.0215 254.1979 135.0945 279.1249
## 36 207.1097 159.4247 254.7946 134.1818 280.0375
## 37 207.1097 158.8353 255.3840 133.2804 280.9389
## 38 207.1097 158.2531 255.9663 132.3899 281.8294
## 39 207.1097 157.6777 256.5417 131.5099 282.7094
## 40 207.1097 157.1089 257.1105 130.6400 283.5793summary(sesHC)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = hc)
##
## Smoothing parameters:
## alpha = 0.3283
##
## Initial states:
## l = 149.2861
##
## sigma: 33.0517
##
## AIC AICc BIC
## 315.8506 316.7737 320.0542
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887 -0.1417763
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 239.5601 197.2026 281.9176 174.7799 304.3403
## 32 239.5601 194.9788 284.1414 171.3788 307.7414
## 33 239.5601 192.8607 286.2595 168.1396 310.9806
## 34 239.5601 190.8347 288.2855 165.0410 314.0792
## 35 239.5601 188.8895 290.2306 162.0662 317.0540
## 36 239.5601 187.0164 292.1038 159.2014 319.9188
## 37 239.5601 185.2077 293.9124 156.4353 322.6848
## 38 239.5601 183.4574 295.6628 153.7584 325.3618
## 39 239.5601 181.7600 297.3602 151.1625 327.9577
## 40 239.5601 180.1111 299.0091 148.6406 330.4795autoplot(sesPB)
autoplot(sesHC)
c. Compute the RMSE values for the training data in each case.
rmsePB <- sqrt(mean(sesPB$residuals^2))
rmseHC <- sqrt(mean(sesHC$residuals^2))rmsePB## [1] 33.63769rmseHC## [1] 31.93101Exercise 6 We will continue with the daily sales of paperback and hardcover books in data set books.
- Apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
holtPB <- holt(pb, h=4)
holtHC <- holt(hc, h=4)summary(holtPB)##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = pb, h = 4)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
##
## Initial states:
## l = 170.699
## b = 1.2621
##
## sigma: 33.4464
##
## AIC AICc BIC
## 318.3396 320.8396 325.3456
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.717178 31.13692 26.18083 -5.508526 15.58354 0.6602122
## ACF1
## Training set -0.1750792
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 209.4668 166.6035 252.3301 143.9130 275.0205
## 32 210.7177 167.8544 253.5811 145.1640 276.2715
## 33 211.9687 169.1054 254.8320 146.4149 277.5225
## 34 213.2197 170.3564 256.0830 147.6659 278.7735autoplot(holtPB)
summary(holtHC)##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = hc, h = 4)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
##
## Initial states:
## l = 147.7935
## b = 3.303
##
## sigma: 29.2106
##
## AIC AICc BIC
## 310.2148 312.7148 317.2208
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1357882 27.19358 23.15557 -2.114792 12.1626 0.6908555
## ACF1
## Training set -0.03245186
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 250.1739 212.7390 287.6087 192.9222 307.4256
## 32 253.4765 216.0416 290.9113 196.2248 310.7282
## 33 256.7791 219.3442 294.2140 199.5274 314.0308
## 34 260.0817 222.6468 297.5166 202.8300 317.3334autoplot(holtHC)
b. Compare the RMSE measures of Holt’s method for the two series to those of simple exponential smoothing in the previous question. (Remember that Holt’s method is using one more parameter than SES.) Discuss the merits of the two forecasting methods for these data sets.
rmsePB.holt <- sqrt(mean(holtPB$residuals^2))
rmseHC.holt <- sqrt(mean(holtHC$residuals^2))rmsePB.holt## [1] 31.13692rmsePB## [1] 33.63769rmseHC.holt## [1] 27.19358rmseHC## [1] 31.93101The RMSE for both paperback and hardcovers is lower in the Holt model. This is noteworthy as the holt model used one more parameter than the SES model. SES is better for forecasting data with no clear trend or seasonal pattern. As there is a clear trend in the data, and somewhat of a seasonal pattern, Holt is a better method to use.
- The Holt method appears to be better based on the prior discussion of RMSE values for both paperback and hardcover. The difference in these models can also be shown in the graphs below.
sesPB <- ses(pb)
holtPB <- holt(pb)
autoplot(books[,1]) + autolayer(sesPB, series = "SES Forecast") + autolayer(holtPB, series = "Holt Forecast", alpha=0.5)+ggtitle("Paperback")
sesHC <- ses(hc)
holtHC <- holt(hc)
autoplot(books[,2]) + autolayer(sesHC, series = "SES Forecast") + autolayer(holtHC, series = "Holt Forecast", alpha=0.5) +ggtitle("Hardcover")
In addition to the RMSE values, the plots for both series show that Holt Forecasts do a better job of capturing the upward trend in each dataset.
- Calculate a 95% prediction interval for the first forecast for each series, using the RMSE values and assuming normal errors. Compare your intervals with those produced using ses and holt.
sesPB self-calculated
sd = sd(sesPB$residuals)
sesPF = sesPB$mean[1]
upper95sesPB = sesPF + 1.96*sd
lower95sesPB = sesPF - 1.96*sdupper95sesPB## [1] 272.623lower95sesPB## [1] 141.5964sesPB$upper[1, "95%"]## 95%
## 275.3523sesPB$lower[1, "95%"]## 95%
## 138.867holtPB self-calculated
sd = sd(holtPB$residuals)
holtPF = holtPB$mean[1]
upper95holtPB = holtPF + 1.96*sd
lower95holtPB = holtPF - 1.96*sdupper95holtPB## [1] 271.0945lower95holtPB## [1] 147.839holtPB$upper[1, "95%"]## 95%
## 275.0205holtPB$lower[1, "95%"]## 95%
## 143.913sesHC self-calculated
sd = sd(sesHC$residuals)
sesPF = sesHC$mean[1]
upper95sesHC = sesPF + 1.96*sd
lower95sesHC = sesPF - 1.96*sdupper95sesHC## [1] 300.5354lower95sesHC## [1] 178.5848sesHC$upper[1, "95%"]## 95%
## 304.3403sesHC$lower[1, "95%"]## 95%
## 174.7799holtHC self-calculated
sd = sd(holtHC$residuals)
holtPF = holtHC$mean[1]
upper95holtHC = holtPF + 1.96*sd
lower95holtHC = holtPF - 1.96*sdupper95holtHC## [1] 304.3838lower95holtHC## [1] 195.964holtHC$upper[1, "95%"]## 95%
## 307.4256holtHC$lower[1, "95%"]## 95%
## 192.9222The self-calculated values are all off by approximately 3-4. This could be a potential result of rounding values differently.
Exercise 7. For this exercise use data set eggs, the price of a dozen eggs in the United States from 1900–1993. Experiment with the various options in the holt() function to see how much the forecasts change with damped trend, or with a Box-Cox transformation. Try to develop an intuition of what each argument is doing to the forecasts.
summary(eggs)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 62.27 148.87 209.15 206.15 276.71 358.78autoplot(eggs)
There is a clear downward trend, and data appears to be yearly. I will start with the most basic Holt model and expand upon it as needed.
holt1 <- holt(eggs, h=100)
autoplot(holt1)
This is highly unrealistic as egg prices would not approach negative numbers. The second model type in the prompt will likely address this issue.
holt2 <- holt(eggs, damped=TRUE, h=100)
autoplot(holt2)
This forecast appears to be much stronger, but the lower boundaries still dip into negative numbers a considerable amount. I will test the addition of a Box-Cox transformation to see if it further restricts this.
holt3 <- holt(eggs, damped = TRUE, lambda = BoxCox.lambda(eggs),h=100)
autoplot(holt3)
The addition of a Box-Cox transformation did an excellent job of restricting the boundaries to be non-negative.
sqrt(mean(holt1$residuals^2))## [1] 26.58219sqrt(mean(holt2$residuals^2))## [1] 26.54019sqrt(mean(holt3$residuals^2))## [1] 1.039187The residuals from model 3 are by far the best out of the three options.
Exercise 8 Recall your retail time series data (from Exercise 3 in Section 2.10).
- Why is multiplicative seasonality necessary for this series?
setwd("C:\\Users\\samne\\OneDrive\\Desktop\\Master\\School\\Grad\\1st year\\Summer\\Forecasting\\Week 4\\HW 2")
retail <-read.csv("retail.csv", stringsAsFactors= T, skip=1)
retailts <- ts(retail[,"A3349873A"],frequency=12, start=c(1982,4))
autoplot(retailts)
The plot of the model clearly shows that the magnitude of the seasonal trend was increasing over time. In the early years the magnitude is incredibly small, but it quickly grows throughout the time period. As the magnitude is changing over time, a multiplicative model best captures these changes.
- Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
retail.holt <- hw(retailts, seasonal = "multiplicative")
retail.damp <- hw(retailts, seasonal = "multiplicative", damped=TRUE)
autoplot(retailts) + autolayer(retail.holt, series = "Holt Forecast") + autolayer(retail.damp, series = "Holt Dampened Forecast", alpha=0.4) +ggtitle("Multiplicative Holt vs. Dampened Multiplicative Holt")
While it is somewhat hard to see, the peaks of the dampened model are much lower than those of the non-dampened model.
- Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
rmse.holt <- sqrt(mean(retail.holt$residuals^2))
rmse.holt## [1] 0.05020746rmse.damp <- sqrt(mean(retail.damp$residuals^2))
rmse.damp## [1] 0.05150921The residuals are very similar, but the non-dampened model is actually slightly better for the time period that is looked at. In a much longer term scenario, the dampened model may be preferable.
- Check that the residuals from the best method look like white noise.
checkresiduals(retail.holt)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 40.405, df = 8, p-value = 2.692e-06
##
## Model df: 16. Total lags used: 24The residuals from the model look largely random and are fairly well contained between small values. The distribution of residuals is largely normal, but it is notable that the bin just above zero is somewhat small compared to the bins around it. Virtually all of the lag on the ACF plot is not significantly different than 0, but there are a handful of points that are. Overall it appears to largely be white noise.
- Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 8 in Section 3.7?
retail.train <- window(retailts, end=c(2010, 12))
retail.test <- window(retailts, start=2011)
retail.snaive <- snaive(retail.train)
accuracy(retail.snaive, retail.test)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.772973 20.24576 15.95676 4.702754 8.109777 1.000000 0.7385090
## Test set 55.300000 71.44309 55.78333 14.900996 15.082019 3.495907 0.5315239
## Theil's U
## Training set NA
## Test set 1.297866retail.holt <- hw(retail.train, seasonal = "multiplicative")
accuracy(retail.holt, retail.test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.03021223 9.107356 6.553533 0.001995484 3.293399 0.4107058
## Test set 55.33771444 70.116586 55.337714 15.173207794 15.173208 3.4679801
## ACF1 Theil's U
## Training set 0.02752875 NA
## Test set 0.35287516 1.30502The two models have a very similar performance. The Holt-Winters model is better using the RMSE metric that has been used throughout this question. The error terms between the two models largely have negligible differences though.
Exercise 9 For the same retail data, try an STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data. How does that compare with your best previous forecasts on the test set?
retail.stl <- stlm(retail.train, s.window=7, lambda= BoxCox.lambda(retail.train), method="ets")
stl.fc <- forecast(retail.stl, lambda=BoxCox.lambda(retail.train))## Warning in InvBoxCox(fcast$mean, lambda, biasadj, fcast): biasadj information
## not found, defaulting to FALSE.autoplot(stl.fc)
accuracy(stl.fc, retail.test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.5906471 7.222165 5.186348 -0.2680566 2.532699 0.3250252
## Test set 55.3633691 70.727261 55.443323 15.0848589 15.114804 3.4745985
## ACF1 Theil's U
## Training set 0.02087455 NA
## Test set 0.40502917 1.306429accuracy(retail.holt, retail.test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.03021223 9.107356 6.553533 0.001995484 3.293399 0.4107058
## Test set 55.33771444 70.116586 55.337714 15.173207794 15.173208 3.4679801
## ACF1 Theil's U
## Training set 0.02752875 NA
## Test set 0.35287516 1.30502This model performed incredibly similarly to the Holt model. The RMSE for this model was slightly worse than the final Holt-Winters model though.
Exercise 10 For this exercise use data set ukcars, the quarterly UK passenger vehicle production data from 1977Q1–2005Q1.
- Plot the data and describe the main features of the series.
autoplot(ukcars)
There appears to be a considerable seasonal trend in the data, and a general trend line that decreases into the mid-1980s and increases after that (with a notable dip after 2000). The seasonal trend appears to decrease in magnitude a considerable amount.
- Decompose the series using STL and obtain the seasonally adjusted data.
uk.stl <- stl(ukcars, s.window = "periodic", robust=TRUE)
uk.stl.sa <- seasadj(uk.stl)
autoplot(uk.stl.sa)
c. Forecast the next two years of the series using an additive damped trend method applied to the seasonally adjusted data. (This can be done in one step using stlf() with arguments etsmodel=“AAN”, damped=TRUE.)
uk.stl.ets <- stlf(uk.stl.sa, etsmodel="AAN", damped=TRUE, h=8)
autoplot(uk.stl.ets)
d. Forecast the next two years of the series using Holt’s linear method applied to the seasonally adjusted data (as before but with damped=FALSE).
uk.holt <- holt(uk.stl.sa, h=8)
autoplot(uk.holt)
e. Now use ets() to choose a seasonal model for the data.
uk.ets <- ets(ukcars)
autoplot(uk.ets)
ETS automatic selection chose an (A, N, A) model.
uk.ets.fc <- forecast(uk.ets, h=8)
autoplot(uk.ets.fc)
- Compare the RMSE of the ETS model with the RMSE of the models you obtained using STL decompositions. Which gives the better in-sample fits?
accuracy(uk.ets)## ME RMSE MAE MPE MAPE MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
## ACF1
## Training set 0.02573334accuracy(uk.stl.ets)## ME RMSE MAE MPE MAPE MASE
## Training set 1.040088 22.81499 17.73846 -0.08623281 5.73168 0.5780878
## ACF1
## Training set 0.01590553accuracy(uk.holt)## ME RMSE MAE MPE MAPE MASE
## Training set -0.2697968 25.32878 20.05443 -0.6328861 6.474306 0.653564
## ACF1
## Training set 0.0276835The RMSE had the lowest value with the STL ETS model. The Holt model and normal ETS models had very similar RMSE values.
- Compare the forecasts from the three approaches? Which seems most reasonable?
autoplot(ukcars) + autolayer(uk.stl.ets, series="STL ETS", PI=FALSE) + autolayer(uk.holt, series="Holt", alpha=0.7, PI=FALSE) + autolayer(uk.ets.fc, series="STL", alpha = 0.7, PI = FALSE)
The STL ETS forecast seems the most reasonable. The Holt forecast is clearly the worst as it does not capture the seasonal trends. The STL forecast is not terrible, but it fails to account for the decrease in magnitude of the seasonal trends. The STL ETS is conservative while still capturing the seasonal fluctuations.
- Check the residuals of your preferred model.
checkresiduals(uk.stl.ets)## Warning in checkresiduals(uk.stl.ets): The fitted degrees of freedom is based on
## the model used for the seasonally adjusted data.
##
## Ljung-Box test
##
## data: Residuals from STL + ETS(A,Ad,N)
## Q* = 23.825, df = 3, p-value = 2.717e-05
##
## Model df: 5. Total lags used: 8The Ljung-Box test indicated that there is serial correlation, and he residuals appear to have some pattern to them. There are some significant points on the ACF plot, and the residuals are somewhat normally distributed.
Exercise 11 For this exercise use data set visitors, the monthly Australian short-term overseas visitors data, May 1985–April 2005.
a.Make a time plot of your data and describe the main features of the series.
autoplot(visitors)
There is a clear upward trend in the data and clear seasonal trends that are increasing in magnitude.
- Split your data into a training set and a test set comprising the last two years of available data. Forecast the test set using Holt-Winters’ multiplicative method.
summary(visitors)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 75.4 189.2 303.1 288.2 378.7 593.1str(visitors) #1:240## Time-Series [1:240] from 1985 to 2005: 75.7 75.4 83.1 82.9 77.3 ...visitors.train <- visitors[1:216]
visitors.test <- visitors[217:240]
visitors.ts <- ts(visitors, frequency = 12, start = c(1985, 5))
visitors.train.ts <- ts(visitors.train, frequency = 12, start = c(1985, 5))
visitors.test.ts <- ts(visitors.test, frequency = 12, start = c(2003, 5))
visitors.holt <- hw(visitors.train.ts, seasonal = "multiplicative", h=24)
summary(visitors.holt)##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = visitors.train.ts, h = 24, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha = 0.4379
## beta = 0.0164
## gamma = 1e-04
##
## Initial states:
## l = 90.9387
## b = 2.584
## s = 0.945 1.0553 1.0882 0.9681 1.3072 1.0711
## 1.0264 0.9095 0.938 1.0401 0.8509 0.8001
##
## sigma: 0.0571
##
## AIC AICc BIC
## 2326.608 2329.699 2383.988
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.9749466 14.06539 10.35763 -0.5792169 4.223204 0.3970304
## ACF1
## Training set 0.1356528
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2003 293.5512 272.0790 315.0235 260.7123 326.3902
## Jun 2003 311.3817 286.3466 336.4168 273.0939 349.6695
## Jul 2003 379.6865 346.4414 412.9316 328.8425 430.5305
## Aug 2003 341.5672 309.2321 373.9023 292.1150 391.0195
## Sep 2003 330.3820 296.7647 363.9994 278.9688 381.7953
## Oct 2003 371.9441 331.4615 412.4267 310.0313 433.8569
## Nov 2003 387.1607 342.2705 432.0509 318.5070 455.8144
## Dec 2003 471.3185 413.3054 529.3316 382.5951 560.0419
## Jan 2004 348.1741 302.8173 393.5309 278.8068 417.5414
## Feb 2004 390.3863 336.7055 444.0672 308.2885 472.4841
## Mar 2004 377.6150 322.9349 432.2951 293.9890 461.2410
## Apr 2004 337.3049 285.9782 388.6316 258.8075 415.8024
## May 2004 284.8762 239.4088 330.3437 215.3398 354.4127
## Jun 2004 302.1571 251.6621 352.6521 224.9316 379.3825
## Jul 2004 368.4105 304.0470 432.7739 269.9751 466.8459
## Aug 2004 331.3981 270.9579 391.8384 238.9627 423.8335
## Sep 2004 320.5215 259.5777 381.4653 227.3160 413.7270
## Oct 2004 360.8154 289.3783 432.2525 251.5618 470.0690
## Nov 2004 375.5478 298.2123 452.8832 257.2733 493.8222
## Dec 2004 457.1458 359.3352 554.9564 307.5574 606.7342
## Jan 2005 337.6781 262.6845 412.6717 222.9853 452.3709
## Feb 2005 378.5881 291.3958 465.7805 245.2390 511.9373
## Mar 2005 366.1740 278.7938 453.5542 232.5375 499.8105
## Apr 2005 327.0594 246.2591 407.8596 203.4861 450.6326autoplot(visitors.holt) +autolayer(visitors.test.ts)
c. Why is multiplicative seasonality necessary here? Multiplicative seasonality is needed because the seasonality changes over time.
- Forecast the two-year test set using each of the following methods:
- an ETS model;
visitors.ets <- ets(visitors.train.ts)
visitors.ets.fc <- forecast(visitors.ets, h =24)
autoplot(visitors.ets.fc)+autolayer(visitors.test.ts)
II .an additive ETS model applied to a Box-Cox transformed series;
visitors.ets.add <- ets(visitors.train.ts, additive.only = TRUE, lambda = BoxCox.lambda(visitors.train.ts))
visitors.ets.add.fc <- forecast(visitors.ets.add, h =24)
autoplot(visitors.ets.add.fc)+autolayer(visitors.test.ts)
III. a seasonal naïve method;
visitors.sn <- snaive(visitors.train.ts)
autoplot(visitors.sn)+autolayer(visitors.test.ts)
IV. an STL decomposition applied to the Box-Cox transformed data followed by an ETS model applied to the seasonally adjusted (transformed) data.
visitors.stl <- stlf(visitors.train.ts, s.window=7, robust=TRUE, method="ets",lambda = BoxCox.lambda(visitors.train.ts))
visitors.stl.fc <- forecast(visitors.stl, lambda = BoxCox.lambda(visitors.train.ts), h=24)
autoplot(visitors.stl.fc)+autolayer(visitors.test.ts)
e. Which method gives the best forecasts? Does it pass the residual tests?
accuracy(visitors.ets.fc, visitors.test.ts)## ME RMSE MAE MPE MAPE MASE
## Training set 0.7640074 14.53480 10.57657 0.1048224 3.994788 0.405423
## Test set 72.1992664 80.23124 74.55285 15.9202832 16.822384 2.857773
## ACF1 Theil's U
## Training set -0.05311217 NA
## Test set 0.58716982 1.127269accuracy(visitors.ets.add.fc, visitors.test.ts)## ME RMSE MAE MPE MAPE MASE
## Training set 1.001363 14.97096 10.82396 0.1609336 3.974215 0.4149057
## Test set 69.458843 78.61032 72.41589 15.1662261 16.273089 2.7758586
## ACF1 Theil's U
## Training set -0.02535299 NA
## Test set 0.67684148 1.086953accuracy(visitors.sn, visitors.test.ts)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 17.29363 31.15613 26.08775 7.192445 10.285961 1.000000 0.6327669
## Test set 32.87083 50.30097 42.24583 6.640781 9.962647 1.619375 0.5725430
## Theil's U
## Training set NA
## Test set 0.6594016accuracy(visitors.stl.fc, visitors.test.ts)## ME RMSE MAE MPE MAPE MASE
## Training set 0.5244387 13.88325 9.060303 0.1194425 3.30638 0.3473011
## Test set 81.7275178 89.34210 82.896081 18.1670388 18.61494 3.1775871
## ACF1 Theil's U
## Training set -0.04096429 NA
## Test set 0.62843225 1.252356The seasonal naive model appears to do the best based on every measurement of error on the test dataset, while the final STL model appears to do the best on the training data. The seasonal naive model will be selected as the test data errors are what truly matters.
checkresiduals(visitors.sn)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 295.02, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24There is a clear issue of autocorrelation with the seasonal naive data, but the residuals appear to be somewhat normally distributed. The ACF plot shows a clear trend and many points are outside of the confidence interval. It does not appear to pass the tests.
- Compare the same four methods using time series cross-validation with the tsCV() function instead of using a training and test set. Do you come to the same conclusions?
fets <- function(y, h) {
forecast(ets(y), h = h)
}
fets.bc <- function(y, h) {
forecast(ets(y, lambda = BoxCox.lambda(y)), h = h)
}
fsn <- function(y) {
forecast(snaive(y), h = 1)
}
fstl <- function(y, h) {
forecast(stlm(y, lambda = BoxCox.lambda(y), s.window = frequency(y) + 1, robust = TRUE, method = "ets"), h = h)
}sqrt(mean(tsCV(visitors, fets, h = 1)^2, na.rm=TRUE))## [1] 18.52985sqrt(mean(tsCV(visitors, fets.bc, h = 1)^2, na.rm=TRUE))## [1] 18.86439sqrt(mean(tsCV(visitors, fstl, h = 1)^2, na.rm = TRUE))## [1] 17.49642sqrt(mean(tsCV(visitors, snaive, h = 1)^2, na.rm = TRUE))## [1] 32.78675The above results indicate that fstl is the best, but all models are relatively similar besides the snaive model. This is the exact opposite conclusion of the prior section as snaive performed best on the test dataset.
Exercise 12 The fets() function below returns ETS forecasts.
- Apply tsCV() for a forecast horizon of
h=4, for both ETS and seasonal naïve methods to the qcement data, (Hint: use the newly created fets() and the existing snaive() functions as your forecast function arguments.)
autoplot(qcement)
fets <- function(y, h) {
forecast(ets(y), h = h)
}qcement.ets <- tsCV(qcement, forecastfunction = fets, h=4)
qcement.sn <- tsCV(qcement, forecastfunction = snaive, h=4)
summary(qcement.ets)## h=1 h=2 h=3 h=4
## Min. :-0.284697 Min. :-0.505359 Min. :-0.467304 Min. :-0.45008
## 1st Qu.:-0.040635 1st Qu.:-0.043311 1st Qu.:-0.052818 1st Qu.:-0.05847
## Median : 0.003049 Median : 0.005511 Median : 0.017342 Median : 0.01804
## Mean : 0.001308 Mean : 0.003296 Mean : 0.002916 Mean : 0.00374
## 3rd Qu.: 0.046562 3rd Qu.: 0.065395 3rd Qu.: 0.076212 3rd Qu.: 0.07748
## Max. : 0.273002 Max. : 0.307940 Max. : 0.316248 Max. : 0.34043
## NA's :1 NA's :2 NA's :3 NA's :4summary(qcement.sn)## h=1 h=2 h=3 h=4
## Min. :-0.41700 Min. :-0.41700 Min. :-0.41700 Min. :-0.41700
## 1st Qu.:-0.01575 1st Qu.:-0.01725 1st Qu.:-0.01825 1st Qu.:-0.01875
## Median : 0.04900 Median : 0.04700 Median : 0.04700 Median : 0.04750
## Mean : 0.03491 Mean : 0.03392 Mean : 0.03394 Mean : 0.03386
## 3rd Qu.: 0.10425 3rd Qu.: 0.10375 3rd Qu.: 0.10425 3rd Qu.: 0.10475
## Max. : 0.35100 Max. : 0.35100 Max. : 0.35100 Max. : 0.35100
## NA's :1 NA's :3 NA's :5 NA's :7- Compute the MSE of the resulting 4-step-ahead errors. (Hint: make sure you remove missing values.) Why are there missing values? Comment on which forecasts are more accurate. Is this what you expected?
ets.mse <- mean(qcement.ets^2, na.rm=TRUE)
ets.mse## [1] 0.01250954sn.mse <- mean(qcement.sn^2, na.rm=TRUE)
sn.mse## [1] 0.01792147Both models have very strong MSE, but the ETS model is slightly stronger. As there is a strong seasonal component that increases in magnitude and a very strong trend in the data, it makes sense that the ETS model would be strongest. The missing values relate to the size of the forecast. Each additional unit of time in the forecast means that there are more unknowns that are estimated.
Exercise 13 Compare ets(), snaive() and stlf() on the following six time series. For stlf(), you might need to use a Box-Cox transformation. Use a test set of three years to decide what gives the best forecasts. ausbeer, bricksq, dole, a10, h02, usmelec.
ausbeer
autoplot(ausbeer)
str(ausbeer) #1:218 quarterly data## Time-Series [1:218] from 1956 to 2010: 284 213 227 308 262 228 236 320 272 233 ...aus.train <- ts(ausbeer[1:206], frequency = 4, start = c(1956,1))
aus.test <- ts(ausbeer[207:218], frequency = 4, start = c(2007,3))
aus.ets <- ets(aus.train)
aus.ets.fc <- forecast(aus.ets, h=12)
aus.sn <- snaive(aus.train, h=12)
aus.stl <- stlf(aus.train, h = 12, lambda = BoxCox.lambda(aus.train))
autoplot(aus.test) + autolayer(aus.ets.fc, series = "ETS", PI=FALSE) + autolayer(aus.sn, series = "Seasonal Naive",PI=FALSE) + autolayer(aus.stl, series = "STL",PI=FALSE)
bricksq
autoplot(bricksq)
str(bricksq) #1:155 quarterly data## Time-Series [1:155] from 1956 to 1994: 189 204 208 197 187 214 227 223 199 229 ...brick.train <- ts(bricksq[1:143], frequency = 4, start = c(1956,1))
brick.test <- ts(bricksq[144:155], frequency = 4, start = c(1991,4))
brick.ets <- ets(brick.train)
brick.ets.fc <- forecast(brick.ets, h=12)
brick.sn <- snaive(brick.train, h=12)
brick.stl <- stlf(brick.train, h = 12, lambda = BoxCox.lambda(brick.train))
autoplot(brick.test) + autolayer(brick.ets.fc, series = "ETS", PI=FALSE) + autolayer(brick.sn, series = "Seasonal Naive",PI=FALSE) + autolayer(brick.stl, series = "STL",PI=FALSE)
dole
autoplot(dole)
str(dole) #1:439 monthly data## Time-Series [1:439] from 1956 to 1992: 4742 6128 6494 5379 6011 ...dole.train <- ts(dole[1:403], frequency = 12, start = c(1956,1))
dole.test <- ts(dole[404:439], frequency = 12, start = c(1989,8))
dole.ets <- ets(dole.train)
dole.ets.fc <- forecast(dole.ets, h=36)
dole.sn <- snaive(dole.train, h=36)
dole.stl <- stlf(dole.train, h = 36, lambda = BoxCox.lambda(dole.train))
autoplot(dole.test) + autolayer(dole.ets.fc, series = "ETS", PI=FALSE) + autolayer(dole.sn, series = "Seasonal Naive",PI=FALSE) + autolayer(dole.stl, series = "STL",PI=FALSE)
a10
autoplot(a10)
str(a10) #204 monthly data## Time-Series [1:204] from 1992 to 2008: 3.53 3.18 3.25 3.61 3.57 ...a10.train <- ts(a10[1:169], frequency = 12, start = c(1992,1))
a10.test <- ts(a10[170:204], frequency = 12, start = c(2006,2))
a10.ets <- ets(a10.train)
a10.ets.fc <- forecast(a10.ets, h=36)
a10.sn <- snaive(a10.train, h=36)
a10.stl <- stlf(a10.train, h = 36, lambda = BoxCox.lambda(a10.train))
autoplot(a10.test) + autolayer(a10.ets.fc, series = "ETS", PI=FALSE) + autolayer(a10.sn, series = "Seasonal Naive",PI=FALSE) + autolayer(a10.stl, series = "STL",PI=FALSE)
h02
autoplot(h02)
str(h02) #204 monthly data, same parameters as a10## Time-Series [1:204] from 1992 to 2008: 0.43 0.401 0.432 0.493 0.502 ...h02.train <- ts(h02[1:169], frequency = 12, start = c(1992,1))
h02.test <- ts(h02[170:204], frequency = 12, start = c(2006,2))
h02.ets <- ets(h02.train)
h02.ets.fc <- forecast(h02.ets, h=36)
h02.sn <- snaive(h02.train, h=36)
h02.stl <- stlf(h02.train, h = 36, lambda = BoxCox.lambda(h02.train))
autoplot(h02.test) + autolayer(h02.ets.fc, series = "ETS", PI=FALSE) + autolayer(h02.sn, series = "Seasonal Naive",PI=FALSE) + autolayer(h02.stl, series = "STL",PI=FALSE)
usmelec
autoplot(usmelec)
str(usmelec) #486 monthly data## Time-Series [1:486] from 1973 to 2013: 160 144 148 140 147 ...usmelec.train <- ts(usmelec[1:450], frequency = 12, start = c(1973,1))
usmelec.test <- ts(usmelec[451:486], frequency = 12, start = c(2010,7))
usmelec.ets <- ets(usmelec.train)
usmelec.ets.fc <- forecast(usmelec.ets, h=36)
usmelec.sn <- snaive(usmelec.train, h=36)
usmelec.stl <- stlf(usmelec.train, h = 36, lambda = BoxCox.lambda(usmelec.train))
autoplot(usmelec.test) + autolayer(usmelec.ets.fc, series = "ETS", PI=FALSE) + autolayer(usmelec.sn, series = "Seasonal Naive",PI=FALSE) + autolayer(usmelec.stl, series = "STL",PI=FALSE)
Exercise 14 a. Use ets() on the following series: bicoal, chicken, dole, usdeaths, lynx, ibmclose, eggs. Does it always give good forecasts?
bicoal
autoplot(forecast(ets(bicoal)))
chicken
autoplot(forecast(ets(chicken)))
dole
autoplot(forecast(ets(dole)))
usdeaths
autoplot(forecast(ets(usdeaths)))
lynx
autoplot(forecast(ets(lynx)))
ibmclose
autoplot(forecast(ets(ibmclose)))
eggs
autoplot(forecast(ets(eggs)))
The bicoal, chicken, ibmclose, and eggs models are not great. These are all straight line predictions that are very similar to a naive model. The lynx plot is similarly flat, but appears to match the data better. Only dole and usdeaths appear to have a forecast that accounts for seasonality.
- Find an example where it does not work well. Can you figure out why?
str(bicoal)## Time-Series [1:49] from 1920 to 1968: 569 416 422 565 484 520 573 518 501 505 ...The bicoal dataset covers yearly data for 49 years. There is no clear seasonality in the data as the data is yearly, and while there appears to be some autocorrelation, there are no clear cyclical trends. These factors mean that there is no clear seasonality to the data and the trend line also is largely unclear.
Exercise 15 Show that the point forecasts from an ETS(M,A,M) model are the same as those obtained using Holt-Winters’ multiplicative method.
I am using the retail data as we have already run a HW multiplicative model
retail.holt <- hw(retailts, seasonal = "multiplicative")
summary(retail.holt)##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = retailts, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha = 0.504
## beta = 1e-04
## gamma = 0.4578
##
## Initial states:
## l = 62.8715
## b = 0.8152
## s = 0.9514 0.886 0.9114 1.5529 1.0184 0.9813
## 0.9589 0.9898 0.9593 0.8883 0.9094 0.9929
##
## sigma: 0.0513
##
## AIC AICc BIC
## 4040.084 4041.770 4107.112
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1170648 13.29378 8.991856 -0.1217735 3.918351 0.4748948
## ACF1
## Training set 0.08635577
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2014 390.3784 364.7154 416.0413 351.1303 429.6264
## Feb 2014 391.1995 362.4039 419.9951 347.1605 435.2386
## Mar 2014 427.9732 393.4376 462.5088 375.1555 480.7909
## Apr 2014 394.1500 359.7834 428.5167 341.5908 446.7093
## May 2014 403.4598 365.8492 441.0704 345.9394 460.9802
## Jun 2014 392.3988 353.6036 431.1940 333.0667 451.7309
## Jul 2014 410.9940 368.1710 453.8169 345.5019 476.4860
## Aug 2014 405.6186 361.3056 449.9315 337.8478 473.3893
## Sep 2014 416.5669 369.0509 464.0828 343.8975 489.2362
## Oct 2014 437.9753 385.9982 489.9524 358.4832 517.4674
## Nov 2014 585.8096 513.6953 657.9240 475.5203 696.0990
## Dec 2014 577.7851 504.1964 651.3737 465.2409 690.3292
## Jan 2015 399.6599 342.8992 456.4206 312.8519 486.4679
## Feb 2015 400.4831 342.1250 458.8412 311.2321 489.7341
## Mar 2015 438.1104 372.6939 503.5270 338.0644 538.1564
## Apr 2015 403.4687 341.8115 465.1258 309.1722 497.7652
## May 2015 412.9807 348.4595 477.5019 314.3041 511.6574
## Jun 2015 401.6414 337.5529 465.7300 303.6264 499.6565
## Jul 2015 420.6566 352.1637 489.1496 315.9057 525.4076
## Aug 2015 415.1371 346.2205 484.0538 309.7383 520.5360
## Sep 2015 426.3243 354.2214 498.4272 316.0524 536.5961
## Oct 2015 448.2152 371.0413 525.3891 330.1879 566.2425
## Nov 2015 599.4807 494.4676 704.4937 438.8771 760.0842
## Dec 2015 591.2440 485.9383 696.5497 430.1928 752.2952The above gives two years worth of data, so I will make an ETS with 2 years of forecasts.
retail.ets <- ets(retailts, model = "MAM")
retail.ets.fc <- forecast(retail.ets, h=24)
summary(retail.ets.fc)##
## Forecast method: ETS(M,Ad,M)
##
## Model Information:
## ETS(M,Ad,M)
##
## Call:
## ets(y = retailts, model = "MAM")
##
## Smoothing parameters:
## alpha = 0.5449
## beta = 0.0101
## gamma = 0.273
## phi = 0.98
##
## Initial states:
## l = 63.1031
## b = 0.3194
## s = 0.9511 0.9049 0.9219 1.5382 1.0807 0.9941
## 0.9575 0.9387 0.9314 0.8936 0.9659 0.9219
##
## sigma: 0.0506
##
## AIC AICc BIC
## 4027.141 4029.030 4098.111
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.8553146 13.9672 8.722404 0.3510144 3.697244 0.460664 -0.06028984
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2014 399.0019 373.1260 424.8777 359.4282 438.5756
## Feb 2014 393.1791 364.0290 422.3291 348.5978 437.7603
## Mar 2014 420.2126 385.4916 454.9335 367.1115 473.3137
## Apr 2014 388.1142 352.9749 423.2536 334.3733 441.8552
## May 2014 402.1398 362.7251 441.5545 341.8602 462.4194
## Jun 2014 394.2989 352.8420 435.7559 330.8960 457.7018
## Jul 2014 411.0768 365.0407 457.1129 340.6706 481.4829
## Aug 2014 412.2801 363.3819 461.1783 337.4967 487.0634
## Sep 2014 422.6129 369.7766 475.4492 341.8068 503.4190
## Oct 2014 440.3068 382.5066 498.1070 351.9091 528.7046
## Nov 2014 563.7567 486.3104 641.2030 445.3128 682.2006
## Dec 2014 592.2239 507.3290 677.1189 462.3883 722.0596
## Jan 2015 412.2604 348.3305 476.1902 314.4881 510.0326
## Feb 2015 405.9474 340.7413 471.1535 306.2233 505.6715
## Mar 2015 433.5499 361.5341 505.5658 323.4112 543.6887
## Apr 2015 400.1549 331.5207 468.7892 295.1879 505.1219
## May 2015 414.3350 341.0524 487.6176 302.2590 526.4110
## Jun 2015 405.9883 332.0341 479.9424 292.8851 519.0914
## Jul 2015 422.9911 343.7248 502.2575 301.7637 544.2186
## Aug 2015 423.9630 342.3158 505.6103 299.0944 548.8317
## Sep 2015 434.3226 348.4482 520.1970 302.9890 565.6562
## Oct 2015 452.2365 360.5159 543.9571 311.9620 592.5111
## Nov 2015 578.6937 458.4011 698.9864 394.7220 762.6655
## Dec 2015 607.5696 478.2277 736.9115 409.7582 805.3809autoplot(retail.test) + autolayer(retail.ets.fc, PI=FALSE, series = "ETS") + autolayer(retail.holt, PI=FALSE, series = "Holt")
The two series are nearly identical apart from minor differences that could be a result of rounding errors.
Exercise 16 Show that the forecast variance for an ETS(A,N,N) model is given by σ2[1+α2(h−1)].
forecast variance = σ2h
This math is outside of my area of expertise, but the basic idea from Hyndman’s other book is that σ2h = vn+h|n for an (A, N, N) model = σ2[1+α2(h−1)]
therefore σ2h=σ2[1+α2(h−1)] for an (A, N, N) model.
Exercise 17. Write down 95% prediction intervals for an ETS(A,N,N) model as a function of
ℓT, α, h and σ, assuming normally distributed errors.
I took the function from pigs and the predictions from that exercise and combined them into a single function to produce confidence intervals using the required metrics.
confint95 = function(y, pars=c(alpha, level), sigma){
fc = 0
SSE = function(y, pars=c(alpha, level)){
yhat = level
error = 0
SSE = 0
alpha = pars[1]
level = pars[2]
for (i in 1:length(y)){
error = y[i]-yhat
SSE = SSE + error^2
yhat = alpha*y[i] + (1-alpha)*yhat
}
fc <<- yhat
return(SSE)
}
optimizedsig = optim(par = c(alpha, pigs[1]), y = sigma, fn = SSE)
return(fc)
sd = sd(optimizedsig$residuals)
sigPF = optimizedsig$mean[1]
upper95 = sigPF + 1.96*sd
lower95 = sigPF - 1.96*sd
return(upper95)
return(lower95)
}###Chapter 8###
Exercise 1 Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.
- Explain the differences among these figures. Do they all indicate that the data are white noise?
While these figures have varying levels of ACF on the y-axis, all of them are within the bounds of the blue lines. These lines indicate that the points are not statistically significantly different than 0, and thus can be considered to be white noise.
b.Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?
The three metrics have vastly different sample sizes, thus the margins of the blue lines vary considerably. Model 1 only has 36 random numbers, model 2 has 360, and model 3 has 1000. The more numbers there are in a sample, the tighter the bounds can be. It is harder to account for white noise in a smaller sample as it is more difficult to tell what actually matters.
Exercise 2 A classic example of a non-stationary series is the daily closing IBM stock price series (data set ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.
autoplot(ibmclose)
checkresiduals(ibmclose)## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
The ACF plot shows strong autocorrelation and all values are outside of the blue lines, meaning that they are significantly different from 0. The value of these points decrease as lag level increases, but all of those present in the graph are significant.
autoplot(pacf(ibmclose))
While most of the values in PACF are not significantly different than 0, the lag of 1 approaches the value 1. This is incredibly significant. Thus both plots indicate autocorrelation and require differencing to remedy.
Exercise 3 For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
- usnetelec
autoplot(usnetelec)
library(urca)## Warning: package 'urca' was built under R version 4.0.5ndiffs(usnetelec)## [1] 1First order differencing is needed.
autoplot(diff(usnetelec))
- usgdp
autoplot(usgdp)
ndiffs(usgdp)## [1] 2autoplot(diff(diff(usgdp)))
c. mcopper
autoplot(mcopper)
ndiffs(mcopper)## [1] 1autoplot(diff(mcopper))
Seems box-cox is still needed
mc.bc <- BoxCox.lambda(mcopper)
autoplot(diff(BoxCox(mcopper, mc.bc)))
- enplanements
autoplot(enplanements)
ndiffs(enplanements)## [1] 1ep.bc <- BoxCox.lambda(enplanements)
autoplot(diff(BoxCox(enplanements, ep.bc)))
seems like a seasonal trend is still present.
nsdiffs(enplanements)## [1] 1autoplot(diff(diff(BoxCox(enplanements, ep.bc),lag = 12)))
This data appears much better than the prior plots for enplacements.
- visitors
autoplot(visitors)
ndiffs(visitors)## [1] 1nsdiffs(visitors)## [1] 1vis.bc <- BoxCox.lambda(visitors)
autoplot(diff(diff(BoxCox(visitors, vis.bc),lag = 12)))
This looks much better than the initial plot.
Exercise 4 For the enplanements data, write down the differences you chose above using backshift operator notation.
I used a first difference and a first seasonal difference. First difference is (1-B)yt first seasonal difference is (1-Bm)yt
Final equation is (1-B)(1.B^m)yt
Exercise 5 For your retail data (from Exercise 3 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data
autoplot(retailts)
ndiffs(retailts)## [1] 1nsdiffs(retailts)## [1] 1retail.bc <- BoxCox.lambda(retailts)
autoplot(diff(diff(BoxCox(retailts, retail.bc),lag = 12)))
Exercise 6 Use R to simulate and plot some data from simple ARIMA models.
- Use the following R code to generate data from an AR(1) model with
ϕ1=0.6 and σ2=1. The process starts with y1=0.
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <- 0.6*y[i-1] + e[i]
}- Produce a time plot for the series. How does the plot change as you change ϕ1?
autoplot(y)
y1 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y1[i] <- 0.06*y1[i-1] + e[i]
}
autoplot(y1)
The data appears to be much tighter and volatile when phi is changed.
c.Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1.
thetaboiz <- function(theta1){
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <- theta1*y[i-1] + e[i]
}
return(y)
}- Produce a time plot for the series. How does the plot change as you change θ1?
autoplot(thetaboiz(0)) + autolayer(thetaboiz(0.25))+ autolayer(thetaboiz(0.5)) + autolayer(thetaboiz(0.75)) + autolayer(thetaboiz(1))
As theta1 increases in value, the variation in value for y increases.
- Generate data from an ARMA(1,1) model with
ϕ1=0.6,θ1=0.6 and σ2=1.
yARMA <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
yARMA[i] <- 0.6*yARMA[i-1] + 0.6*e[i-1] + e[i]
}
plot(yARMA)
f. Generate data from an AR(2) model with ϕ1=−0.8, ϕ2=0.3 and σ2 = 1.(Note that these parameters will give a non-stationary series.)
yAR2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
yAR2[i] <- -0.8*yAR2[i-1] + 0.3*yAR2[i-1] + e[i]
}
autoplot(yAR2)
g. Graph the latter two series and compare them.
autoplot(yARMA)
autoplot(yAR2) The second model is much more volatile than the first, though the magnitudes in the first model are greater.
Exercise 7 Consider wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.
autoplot(wmurders)
ndiffs(wmurders)## [1] 2autoplot(diff(diff(wmurders)))
ur.kpss(diff(diff(wmurders)))##
## #######################################
## # KPSS Unit Root / Cointegration Test #
## #######################################
##
## The value of the test statistic is: 0.0458Unit root test on twice differenced data indicates that it is now stationary and that d should be 2
wm2d <- diff(diff(wmurders))
checkresiduals(wm2d)## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
autoplot(pacf(wm2d))
An AR(1) model seems to work based on the PACF. Thus and ARIMA(1,2,0) model will be selected
Should you include a constant in the model? Explain. I would not include a constant in the model as the data is twice differenced. Such a constant would relate to a quadratic trend in the actual data and could change the model a great deal.
Write this model in terms of the backshift operator.
((1−B)^2)yt = C + 1 − φ1B
d.Fit the model using R and examine the residuals. Is the model satisfactory?
wmAR120 <- Arima(wmurders, order=c(1,2,0))
checkresiduals(wmAR120)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,2,0)
## Q* = 16.452, df = 9, p-value = 0.05802
##
## Model df: 1. Total lags used: 10The residuals appear to be white noise, but the distribution of residuals is somewhat abnormal. This may just be a factor of the bin size though. The residual plot appears to have considerable correlation, but the p-value of the Ljung-Box plot is not enough to reject the null hypothesis (though it is very close). The model is not amazing, but it is satisfactory.
- Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.
wmAR120.fc <- forecast(wmAR120, h=3)
wmAR120.fc## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 2.474461 2.174701 2.774221 2.016018 2.932904
## 2006 2.387811 1.889464 2.886158 1.625655 3.149966
## 2007 2.282165 1.477486 3.086843 1.051515 3.512814wmAR120.fc$mean[1]## [1] 2.474461The point forecast for the first point matches the mean for h=1.
- Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
autoplot(wmAR120.fc)
g. Does auto.arima() give the same model you have chosen? If not, which model do you think is better?
wmaa <- auto.arima(wmurders)auto.arima yields a (1,2,1) model.
wmaa.fc <- forecast(wmaa, h=3)
autoplot(wmaa.fc) + autolayer(wmAR120.fc, series = "(1,2,0) Model", alpha = 0.4)
accuracy(wmaa.fc)## ME RMSE MAE MPE MAPE MASE
## Training set -0.01065956 0.2072523 0.1528734 -0.2149476 4.335214 0.9400996
## ACF1
## Training set 0.02176343accuracy(wmAR120.fc)## ME RMSE MAE MPE MAPE MASE
## Training set -0.001376898 0.2274352 0.1777919 0.001486708 5.080777 1.093337
## ACF1
## Training set -0.1593845Both models are very similar, but the auto.arima model had slightly better RMSE and MASE. The MASE on the auto.arima model is less than 0, while the one for the model I selected is slightly above 1. Thus the auto.arima model is slightly better.
Exercise 8 Consider austa, the total international visitors to Australia (in millions) for the period 1980-2015. a. Use auto.arima() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
autoaus <- auto.arima(austa)
autoaus## Series: austa
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.3006 0.1735
## s.e. 0.1647 0.0390
##
## sigma^2 estimated as 0.03376: log likelihood=10.62
## AIC=-15.24 AICc=-14.46 BIC=-10.57(0,1,1) with drift was selected.
checkresiduals(autoaus)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1) with drift
## Q* = 2.297, df = 5, p-value = 0.8067
##
## Model df: 2. Total lags used: 7The Ljung-Box test indicates that there is no significant proof of autocorrelation and the residual plots indicates a largely normal distribution of residuals. The ACF plot indicates that the residuals are white noise.
autoaus.fc <- forecast(auto.arima(austa), h=10)
autoplot(autoaus.fc)
b. Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part a. Remove the MA term and plot again.
arimawithnobrim <- forecast(Arima(austa, order=c(0,1,1)), h=10)
autoplot(arimawithnobrim)
With the removal of drift, the arima (0,1,1) model is a flat line that is effectively a naive model (though it is slightly higher than the last point).
ari.fc <- forecast(Arima(austa, order=c(0,1,0)), h=10)
autoplot(ari.fc)
This model is similar to the prior one, but appears to just take the last known value like a naive model.
- Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.
arima213d <- forecast(Arima(austa, c(2,1,3), include.drift = TRUE), h=10)
autoplot(arima213d)
#arima213 <- forecast(Arima(austa, c(2,1,3), include.drift = FALSE), h=10)“Error in stats::arima(x = x, order = order, seasonal = seasonal, include.mean = include.mean, : non-stationary AR part from CSS” If you attempt to remove the constant, you get the above error. The constant is needed for the model.
- Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.
arima001c <- forecast(Arima(austa, c(0,0,1), include.constant = TRUE), h=10)
autoplot(arima001c)
arima000c <- forecast(Arima(austa, c(0,0,0), include.constant = TRUE), h=10)
autoplot(arima000c)
e. Plot forecasts from an ARIMA(0,2,1) model with no constant.
arima021 <- forecast(Arima(austa, c(0,2,1), include.constant = FALSE), h=10)
autoplot(arima021)
Exercise 9 For the usgdp series: a. if necessary, find a suitable Box-Cox transformation for the data;
autoplot(usgdp) Strong upward trend
autoplot(BoxCox(usgdp, BoxCox.lambda(usgdp)))
Box-Cox is slightly smoother, but also on a much more manageable scale.
usBox <- BoxCox.lambda(usgdp)- fit a suitable ARIMA model to the transformed data using auto.arima();
usAuto <- auto.arima(usgdp, lambda = usBox)
usAuto## Series: usgdp
## ARIMA(2,1,0) with drift
## Box Cox transformation: lambda= 0.366352
##
## Coefficients:
## ar1 ar2 drift
## 0.2795 0.1208 0.1829
## s.e. 0.0647 0.0648 0.0202
##
## sigma^2 estimated as 0.03518: log likelihood=61.56
## AIC=-115.11 AICc=-114.94 BIC=-101.26autoplot(usAuto$fitted) + autolayer(usgdp, alpha = 0.7)
The model fits the data incredibly well.
- Try some other plausible models by experimenting with the orders chosen; The model chosen was (2,1,0) with drift. I will try (2,1,1)d, (1,1,0)d, and (2,1,0)
us211d <- Arima(usgdp, order = c(2, 1, 1), lambda = usBox, include.drift = TRUE)
autoplot(us211d$fitted) + autolayer(usgdp, alpha = 0.7)
us110d <- Arima(usgdp, order = c(1, 1, 0), lambda = usBox, include.drift = TRUE)
autoplot(us110d$fitted) + autolayer(usgdp, alpha = 0.7)
us210 <- Arima(usgdp, order = c(2, 1, 0), lambda = usBox, include.drift = FALSE)
autoplot(us210$fitted) + autolayer(usgdp, alpha = 0.7)
All of the models have incredibly similar performance.
- choose what you think is the best model and check the residual diagnostics;
accuracy(usAuto)## ME RMSE MAE MPE MAPE MASE
## Training set 1.195275 39.2224 29.29521 -0.01363259 0.6863491 0.1655687
## ACF1
## Training set -0.03824844accuracy(us211d)## ME RMSE MAE MPE MAPE MASE
## Training set 1.224559 39.08125 29.27782 -0.01431716 0.683324 0.1654704
## ACF1
## Training set -0.03850778accuracy(us210)## ME RMSE MAE MPE MAPE MASE
## Training set 11.06831 42.21158 32.00633 0.2244843 0.7393035 0.1808912
## ACF1
## Training set -0.1210488accuracy(us110d)## ME RMSE MAE MPE MAPE MASE
## Training set 1.315796 39.90012 29.5802 -0.01678591 0.6834509 0.1671794
## ACF1
## Training set -0.08544569All models have similar errors. The auto model has the lowest mean error, but the (2,1,1) with drift model has the lowest metrics related to variance. It also has the lowest MASE (by a very small margin). I will select the (2,1, 1) with drift model as it has the best metrics for forecasting.
checkresiduals(us211d)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1) with drift
## Q* = 5.4997, df = 4, p-value = 0.2398
##
## Model df: 4. Total lags used: 8The Ljung-Box test indicates no serial correlation and the distribution of means is largely normal. The ACF plot indicates that most of the residuals are white noise, save for lag 12. This may mean that there is some seasonality that is not being accounted for.
- produce forecasts of your fitted model. Do the forecasts look reasonable?
us211d.fc <- forecast(us211d)
autoplot(us211d.fc)
The forecast looks reasonable and appears to follow the current upward trend well.
- compare the results with what you would obtain using ets() (with no transformation).
usETS.fc <- forecast(ets(usgdp))
autoplot(usETS.fc)
The ETS model appears to do a better job of capturing the upward trend line seen in the data.
Exercise 10 Consider austourists, the quarterly visitor nights (in millions) spent by international tourists to Australia for the period 1999–2015.
- Describe the time plot.
autoplot(austourists)
There is a clear seasonal pattern with an upward trend. The magnitude of the data also appears to be increasing with time.
- What can you learn from the ACF graph?
autoplot(Acf(austourists))
The data is decreasing with significant spikes at each yearly mark. There is an indication that there is strong seasonal trends. Data for the first 12 lag levels are significant and lag level 16 is significant.
- What can you learn from the PACF graph?
autoplot(Pacf(austourists))
This chart also indicates a strong yearly pattern. There are no significant spikes after lag 8, and most of the spikes are within the first 5 levels (lag 3 does not have a significant spike). This again indicates strong seasonality.
- Produce plots of the seasonally differenced data (1−B4)Yt. What model do these graphs suggest?
As the data is quarterly, we are effectively doing a yearly seasonal lag.
aus.season <-diff(austourists, lag=4, differences=1)
autoplot(aus.season)
ur.kpss(aus.season)##
## #######################################
## # KPSS Unit Root / Cointegration Test #
## #######################################
##
## The value of the test statistic is: 0.2575Data needs to be differenced again.
aus.season.dif <- diff(aus.season)
ur.kpss(aus.season.dif)##
## #######################################
## # KPSS Unit Root / Cointegration Test #
## #######################################
##
## The value of the test statistic is: 0.0359autoplot(aus.season.dif)
The graphs and tests suggest that the arima model will need to factor in first differencing and have a seasonal component with a lag of 4.
- Does auto.arima() give the same model that you chose? If not, which model do you think is better?
autoaus<- auto.arima(austourists)
autoaus## Series: austourists
## ARIMA(1,0,0)(1,1,0)[4] with drift
##
## Coefficients:
## ar1 sar1 drift
## 0.4705 -0.5305 0.5489
## s.e. 0.1154 0.1122 0.0864
##
## sigma^2 estimated as 5.15: log likelihood=-142.48
## AIC=292.97 AICc=293.65 BIC=301.6The auto.arima model matches the specifications that I had from part d.
- Write the model in terms of the backshift operator, then without using the backshift operator.
autoaus$model## $phi
## [1] 0.4704867 0.0000000 0.0000000 -0.5304552 0.2495721
##
## $theta
## [1] 0 0 0 0
##
## $Delta
## [1] 0 0 0 1
##
## $Z
## [1] 1 0 0 0 0 0 0 0 1
##
## $a
## [1] 4.0257726 -3.0452679 0.1557664 -1.5879471 1.0336022 24.3229131 11.4706473
## [8] 37.5799294 24.7062553
##
## $P
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000000e+00 -3.011281e-17 -1.410195e-18 5.724459e-18 4.322798e-19
## [2,] -3.011281e-17 1.348336e-17 1.222393e-18 1.275158e-18 -1.572649e-18
## [3,] -1.410195e-18 1.222393e-18 -8.190646e-19 1.807235e-18 -7.596270e-19
## [4,] 5.724459e-18 1.275158e-18 1.807235e-18 -2.897583e-18 6.816372e-19
## [5,] 4.322798e-19 -1.572649e-18 -7.596270e-19 6.816372e-19 -5.356201e-38
## [6,] 1.732084e-18 -9.848056e-22 3.630797e-23 2.419025e-19 -1.138095e-19
## [7,] 2.661857e-17 -1.001249e-21 5.248740e-18 -8.274611e-18 2.731223e-18
## [8,] 5.092621e-17 -1.851678e-17 7.069789e-18 -5.097169e-18 2.761385e-18
## [9,] -1.241618e-17 3.011131e-17 1.410210e-18 -5.724445e-18 -4.321738e-19
## [,6] [,7] [,8] [,9]
## [1,] 1.730204e-18 2.661599e-17 5.092091e-17 -1.241698e-17
## [2,] 2.798843e-35 3.265529e-34 -1.851485e-17 3.011281e-17
## [3,] -1.993554e-38 5.248928e-18 7.069203e-18 1.410195e-18
## [4,] 2.419614e-19 -8.274619e-18 -5.096903e-18 -5.724459e-18
## [5,] -1.138396e-19 2.731207e-18 2.761448e-18 -4.322798e-19
## [6,] -4.561392e-19 1.094356e-17 1.106473e-17 -1.732084e-18
## [7,] -4.174965e-35 -9.895140e-18 3.095056e-18 -2.661857e-17
## [8,] -8.890727e-35 -7.460735e-34 3.490371e-17 -5.092621e-17
## [9,] 5.875338e-34 -6.066331e-35 9.017924e-35 1.241618e-17
##
## $T
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.4704867 1 0 0 0 0 0 0 0
## [2,] 0.0000000 0 1 0 0 0 0 0 0
## [3,] 0.0000000 0 0 1 0 0 0 0 0
## [4,] -0.5304552 0 0 0 1 0 0 0 0
## [5,] 0.2495721 0 0 0 0 0 0 0 0
## [6,] 1.0000000 0 0 0 0 0 0 0 1
## [7,] 0.0000000 0 0 0 0 1 0 0 0
## [8,] 0.0000000 0 0 0 0 0 1 0 0
## [9,] 0.0000000 0 0 0 0 0 0 1 0
##
## $V
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 0
##
## $h
## [1] 0
##
## $Pn
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000e+00 -1.572148e-17 -1.420446e-18 3.700301e-18 2.008451e-19
## [2,] -1.572148e-17 1.348336e-17 1.222393e-18 1.275158e-18 -1.572649e-18
## [3,] -1.420446e-18 1.222393e-18 -8.190646e-19 1.807235e-18 -7.596270e-19
## [4,] 3.700301e-18 1.275158e-18 1.807235e-18 -2.897583e-18 6.816372e-19
## [5,] 2.008451e-19 -1.572649e-18 -7.596270e-19 6.816372e-19 0.000000e+00
## [6,] -2.127279e-19 -9.848056e-22 3.630797e-23 2.419025e-19 -1.138095e-19
## [7,] 5.151377e-18 -1.001249e-21 5.248740e-18 -8.274611e-18 2.731223e-18
## [8,] 5.211112e-18 -1.851678e-17 7.069789e-18 -5.097169e-18 2.761385e-18
## [9,] -7.400349e-18 3.011131e-17 1.410210e-18 -5.724445e-18 -4.321738e-19
## [,6] [,7] [,8] [,9]
## [1,] -2.146074e-19 5.148800e-18 5.205808e-18 -7.401150e-18
## [2,] 0.000000e+00 1.761151e-35 -1.851485e-17 3.011281e-17
## [3,] 0.000000e+00 5.248928e-18 7.069203e-18 1.410195e-18
## [4,] 2.419614e-19 -8.274619e-18 -5.096903e-18 -5.724459e-18
## [5,] -1.138396e-19 2.731207e-18 2.761448e-18 -4.322798e-19
## [6,] -4.561392e-19 1.094356e-17 1.106473e-17 -1.732084e-18
## [7,] 0.000000e+00 -9.895140e-18 3.095056e-18 -2.661857e-17
## [8,] 0.000000e+00 2.353014e-34 3.490371e-17 -5.092621e-17
## [9,] 5.777790e-34 -1.683391e-34 -1.391200e-34 1.241618e-17(1-B)(1−B4)Yt = c+ϕ1y(t−1) +εt
Exercise 11 Consider usmelec, the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 – June 2013). In general there are two peaks per year: in mid-summer and mid-winter.
- Examine the 12-month moving average of this series to see what kind of trend is involved.
usmelec12ma <- ma(usmelec, order = 12)
autoplot(usmelec12ma)
b. Do the data need transforming? If so, find a suitable transformation.
The magnitude changes slightly in recent years, but not by much.
usmelec.bc <- BoxCox.lambda(usmelec)
autoplot(BoxCox(usmelec, usmelec.bc))
c. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
ndiffs(usmelec)## [1] 1nsdiffs(usmelec)## [1] 1Data needs season and regular differences.
usmelec.stat <- diff(diff(BoxCox(usmelec, usmelec.bc),lag = 12))
autoplot(usmelec.stat)
ur.kpss(usmelec.stat)##
## #######################################
## # KPSS Unit Root / Cointegration Test #
## #######################################
##
## The value of the test statistic is: 0.0167- Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
auto.arima(usmelec, lambda = usmelec.bc)## Series: usmelec
## ARIMA(1,1,3)(2,1,1)[12]
## Box Cox transformation: lambda= -0.5738331
##
## Coefficients:
## ar1 ma1 ma2 ma3 sar1 sar2 sma1
## 0.3952 -0.8194 -0.0468 0.0414 0.0403 -0.0934 -0.8462
## s.e. 0.3717 0.3734 0.1707 0.1079 0.0560 0.0533 0.0343
##
## sigma^2 estimated as 1.278e-06: log likelihood=2547.75
## AIC=-5079.5 AICc=-5079.19 BIC=-5046.23auto.arima(usmelec)## Series: usmelec
## ARIMA(1,0,2)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 ma1 ma2 sma1 drift
## 0.9717 -0.4374 -0.2774 -0.7061 0.3834
## s.e. 0.0163 0.0483 0.0493 0.0310 0.0868
##
## sigma^2 estimated as 57.67: log likelihood=-1635.13
## AIC=3282.26 AICc=3282.44 BIC=3307.22arima2 <- Arima(usmelec, lambda = usmelec.bc, order = c(1, 0, 2), seasonal = c(0, 1, 1))
arima2## Series: usmelec
## ARIMA(1,0,2)(0,1,1)[12]
## Box Cox transformation: lambda= -0.5738331
##
## Coefficients:
## ar1 ma1 ma2 sma1
## 0.9987 -0.4311 -0.2545 -0.8528
## s.e. 0.0016 0.0440 0.0441 0.0265
##
## sigma^2 estimated as 1.28e-06: log likelihood=2550.34
## AIC=-5090.69 AICc=-5090.56 BIC=-5069.88arima3 <- Arima(usmelec, lambda = usmelec.bc, order = c(1, 1, 2), seasonal = c(0, 1, 1))
arima3## Series: usmelec
## ARIMA(1,1,2)(0,1,1)[12]
## Box Cox transformation: lambda= -0.5738331
##
## Coefficients:
## ar1 ma1 ma2 sma1
## 0.2411 -0.6625 -0.1151 -0.8577
## s.e. 0.1856 0.1887 0.1261 0.0258
##
## sigma^2 estimated as 1.282e-06: log likelihood=2545.46
## AIC=-5080.93 AICc=-5080.8 BIC=-5060.13arima4 <- Arima(usmelec, lambda = usmelec.bc, order = c(0, 1, 2), seasonal = c(0, 1, 1))
arima4## Series: usmelec
## ARIMA(0,1,2)(0,1,1)[12]
## Box Cox transformation: lambda= -0.5738331
##
## Coefficients:
## ma1 ma2 sma1
## -0.4317 -0.2552 -0.8536
## s.e. 0.0439 0.0440 0.0261
##
## sigma^2 estimated as 1.284e-06: log likelihood=2544.75
## AIC=-5081.51 AICc=-5081.42 BIC=-5064.87The auto.arima model is the best in terms of AIC of the models tested above.
- Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
checkresiduals(auto.arima(usmelec))
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,1)[12] with drift
## Q* = 42.725, df = 19, p-value = 0.001413
##
## Model df: 5. Total lags used: 24There is significant proof of autocorrelation in the model and the ACF plot indicates that there are a few points of lag that are significantly different than 0.
I will use the last model I made above to see if it has better results.
checkresiduals(arima4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 32.525, df = 21, p-value = 0.05176
##
## Model df: 3. Total lags used: 24The Ljung-Box test indicates that the data is not autocorrelated, but it is very close to being significantly correlated. There is only one significant spike in the ACF plot, and the distribution of residuals appears to largely be normal. Overall, this is a better plot.
- Forecast the next 15 years of electricity generation by the U.S. electric industry. Get the latest figures from the EIA to check the accuracy of your forecasts.
arima4.fc <- forecast(arima4, h = 180)
autoplot(arima4.fc)
The forecast is not very accurate when comparing it to the charts on the EIC website. This makes sense though as forecasting 180 points out is highly excessive.
- Eventually, the prediction intervals are so wide that the forecasts are not particularly useful. How many years of forecasts do you think are sufficiently accurate to be usable?
The accuracy of a forecast will depend on a great deal of factors. The volatility of the data and the varying influence of exogenous factors on models means that there is not one clear answer to this question. In general forecasts are more accurate the more data there is to build them and the shorter they are out from the actual data. There are some situations where minor fluctuations are hard to calculate, but the general trend is clear. For instance, it would be safe to assume that the price of the S&P 500 will be higher in 50 years than it is today.
Exercise 12 For the mcopper data:
- if necessary, find a suitable Box-Cox transformation for the data;
autoplot(mcopper)
mcopper.bc <- BoxCox.lambda(mcopper)
autoplot(BoxCox(mcopper, lambda = mcopper.bc))
The variation in the data is somewhat cleaner after the box-cox transformation.
- fit a suitable ARIMA model to the transformed data using auto.arima();
mcopper.auto<-auto.arima(mcopper, lambda = mcopper.bc)
mcopper.auto## Series: mcopper
## ARIMA(0,1,1)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ma1
## 0.3720
## s.e. 0.0388
##
## sigma^2 estimated as 0.04997: log likelihood=45.05
## AIC=-86.1 AICc=-86.08 BIC=-77.43- try some other plausible models by experimenting with the orders chosen;
Arima(mcopper, order=c(0,1,2), lambda=mcopper.bc)## Series: mcopper
## ARIMA(0,1,2)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ma1 ma2
## 0.3704 -0.0037
## s.e. 0.0424 0.0406
##
## sigma^2 estimated as 0.05006: log likelihood=45.05
## AIC=-84.11 AICc=-84.06 BIC=-71.11Arima(mcopper, order=c(1,1,1), lambda=mcopper.bc)## Series: mcopper
## ARIMA(1,1,1)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ar1 ma1
## -0.0092 0.3797
## s.e. 0.1053 0.0961
##
## sigma^2 estimated as 0.05006: log likelihood=45.05
## AIC=-84.1 AICc=-84.06 BIC=-71.1Arima(mcopper, order=c(1,1,2), lambda=mcopper.bc)## Series: mcopper
## ARIMA(1,1,2)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ar1 ma1 ma2
## -0.8974 1.2828 0.3609
## s.e. 0.1839 0.1787 0.0642
##
## sigma^2 estimated as 0.05002: log likelihood=45.76
## AIC=-83.52 AICc=-83.45 BIC=-66.19Arima(mcopper, order=c(0,1,0), lambda=mcopper.bc)## Series: mcopper
## ARIMA(0,1,0)
## Box Cox transformation: lambda= 0.1919047
##
## sigma^2 estimated as 0.05674: log likelihood=8.85
## AIC=-15.7 AICc=-15.69 BIC=-11.37The auto.arima model is the strongest out of the models tested when comparing AIC values.
- choose what you think is the best model and check the residual diagnostics;
checkresiduals(mcopper.auto)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)
## Q* = 22.913, df = 23, p-value = 0.4659
##
## Model df: 1. Total lags used: 24There is no proof of autocorrelation, the residuals are mostly normally distributed, and the ACF plot shows no significant spikes.
- produce forecasts of your fitted model. Do the forecasts look reasonable?
mcopper.auto.fc <- forecast(mcopper.auto)
autoplot(mcopper.auto.fc)
The forecasts do not look very reasonable.
- compare the results with what you would obtain using ets() (with no transformation).
mcopper.ets<-ets(mcopper)
autoplot(mcopper.ets)
mcopper.ets.fc<-forecast(mcopper.ets)
autoplot(mcopper.ets.fc)
The forecast seems better than the arima model as the arima model looks almost like a naive model.
Exercise 13 Choose one of the following seasonal time series: hsales, auscafe, qauselec, qcement, qgas.
- Do the data need transforming? If so, find a suitable transformation.
autoplot(qcement)
qcement.bc <- BoxCox.lambda(qcement)
autoplot(BoxCox(qcement, lambda=qcement.bc))
BoxCox helped bring variance closer together.
- Are the data stationary? If not, find an appropriate differencing which yields stationary data.
ndiffs(qcement)## [1] 1nsdiffs(qcement)## [1] 1autoplot()
qcement.sa <- diff(diff(BoxCox(qcement, qcement.bc),lag = 12))
autoplot(qcement.sa)
ur.kpss(qcement.sa)##
## #######################################
## # KPSS Unit Root / Cointegration Test #
## #######################################
##
## The value of the test statistic is: 0.0302c.Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
qcement.auto<-auto.arima(qcement, lambda = qcement.bc)
qcement.auto## Series: qcement
## ARIMA(1,1,0)(1,1,2)[4]
## Box Cox transformation: lambda= -0.4319725
##
## Coefficients:
## ar1 sar1 sma1 sma2
## -0.2174 -0.5432 -0.2417 -0.5970
## s.e. 0.0660 0.2202 0.1998 0.1719
##
## sigma^2 estimated as 0.001643: log likelihood=406.17
## AIC=-802.35 AICc=-802.08 BIC=-785.2qcement.arima <- Arima(qcement, lambda = qcement.bc, order=c(1,1,1))
qcement.arima## Series: qcement
## ARIMA(1,1,1)
## Box Cox transformation: lambda= -0.4319725
##
## Coefficients:
## ar1 ma1
## 0.1075 -0.6287
## s.e. 0.0950 0.0647
##
## sigma^2 estimated as 0.006341: log likelihood=258.66
## AIC=-511.32 AICc=-511.22 BIC=-500.98Arima(qcement, lambda = qcement.bc, order=c(0,1,0))## Series: qcement
## ARIMA(0,1,0)
## Box Cox transformation: lambda= -0.4319725
##
## sigma^2 estimated as 0.007957: log likelihood=231.51
## AIC=-461.03 AICc=-461.01 BIC=-457.58The ARIMA model (1,1,0)(1,1,2)[4] is the best according to AIC
- Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
checkresiduals(qcement.auto)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(1,1,2)[4]
## Q* = 9.6635, df = 4, p-value = 0.04649
##
## Model df: 4. Total lags used: 8The plots indicate significant autocorrelation and lag spikes at the 17 month level.
checkresiduals(qcement.arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)
## Q* = 274.7, df = 6, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 8This plot has much stronger autocorrelation and many more significant lag levels. The auto appears to be better.
- Forecast the next 24 months of data using your preferred model.
qcement.auto.fc <- forecast(qcement.auto, h=24)
autoplot(qcement.auto.fc)
f. Compare the forecasts obtained using ets().
qcement.ets <- ets(qcement)
autoplot(qcement.ets)
qcement.ets.fc<- forecast(qcement.ets, h=24)
autoplot(qcement.ets.fc) + autolayer(qcement.auto.fc, series="ARIMA", PI=FALSE)
The ETS plot is very similar to the ARIMA plot, but the ETS plot is more conservative.
Exercise 14 For the same time series you used in the previous exercise, try using a non-seasonal model applied to the seasonally adjusted data obtained from STL. The stlf() function will make the calculations easy (with method=“arima”). Compare the forecasts with those obtained in the previous exercise. Which do you think is the best approach?
qcement.stl <- stlf(qcement, s.window = 13, robust=TRUE, method="arima", h=24)
autoplot(qcement.stl)
autoplot(qcement.stl) + autolayer(qcement.auto.fc, series="ARIMA", PI=FALSE) +autolayer(qcement.ets.fc, series="ETS", PI=FALSE)
The ETS and STL models appear to be very similar, while the ARIMA model appears to be much larger. I believe that the ETS or STL are better approaches than the pure ARIMA as they better capture the trends in the data.
Exercise 15 For your retail time series (Exercise 5 above): a.develop an appropriate seasonal ARIMA model;
retail.arima<-auto.arima(retailts)
retail.arima.fc<-forecast(retail.arima)
autoplot(retail.arima.fc) 
autoplot(retail.arima.fc) +autolayer(forecast(snaive(retailts)), series = "Seasonal naive", PI=FALSE)
The plot compared to prior sections plot of seasonal naive data appears to be very similar, but somewhat better at capturing the data.
c,Obtain up-to-date retail data from the ABS website (Cat 8501.0, Table 11), and compare your forecasts with the actual numbers. How good were the forecasts from the various models?
I had difficulty downloading the data, so I decided to make a test dataset out of the current data instead.
train <- subset(retailts, end = length(retailts) - 12)
retail.arima1<-auto.arima(train)
retail.arima1.fc<-forecast(retail.arima1, h=12)
retail.snaive1 <-forecast(snaive(train), h=12)
autoplot(retailts) +autolayer(retail.snaive1, series = "Seasonal naive", PI=FALSE) + autolayer(retail.arima1.fc, series = "ARIMA", PI=FALSE)
The ARIMA forecast is much closer to the actual data.
Exercise 16 Consider sheep, the sheep population of England and Wales from 1867–1939.
- Produce a time plot of the time series.
autoplot(sheep)
b. Assume you decide to fit the following model: yt=yt−1+ϕ1(yt−1−yt−2)+ϕ2(yt−2−yt−3)+ϕ3(yt−3−yt−4)+εt, where εt is a white noise series. What sort of ARIMA model is this (i.e., what are p ,d , and q)?
The three ϕ mean it is a (3,x,x) model
The single Yt-1 outside of ϕ terms means it is a (x,1,x) model
Model should be (3,1,0)
c.By examining the ACF and PACF of the differenced data, explain why this model is appropriate.
checkresiduals(diff(sheep))## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
pacf(diff(sheep))
Pacf implies that there are no signifcant lag amounts after 3. This means a (3,x,x) model is warranted.
d.The last five values of the series are given below:
The estimated parameters are
ϕ1=0.42,ϕ2=−0.20, and ϕ3=−0.30. Without using the forecast function, calculate forecasts for the next three years (1940–1942).
sheep1940 <- 1797 + 0.42*(1797-1791) - 0.2*(1791-1627) -0.3*(1627-1665)
sheep1941 <- sheep1940 + 0.42*(sheep1940-1797) -0.20*(1797-1791) - 0.30*(1791-1627)
sheep1942 <- sheep1941 + 0.42*(sheep1941-sheep1940) - 0.20*(sheep1940-1797) - 0.30*(1797-1791)
sheep1940## [1] 1778.12sheep1941## [1] 1719.79sheep1942## [1] 1697.268- Now fit the model in R and obtain the forecasts using forecast. How are they different from yours? Why?
forecast(Arima(sheep, order=c(3,1,0)), h=3)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1940 1777.996 1687.454 1868.537 1639.524 1916.467
## 1941 1718.869 1561.544 1876.193 1478.261 1959.476
## 1942 1695.985 1494.151 1897.819 1387.306 2004.664The predictions are very similar, but are slightly off. This is likely a result of rounding in the ϕ values.
Exercise 17 The annual bituminous coal production in the United States from 1920 to 1968 is in data set bicoal.
autoplot(bicoal)
b. You decide to fit the following model to the series: yt=c+ϕ1yt−1+ϕ2yt−2+ϕ3yt−3+ϕ4yt−4+εt where yt is the coal production in year t and εt is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?
There are four ϕ terms in the model meaning it is a (4,x,x) model There are zero Yt-1 terms outside of ϕ terms means it is a (x,0,x) model The model is a (4,0,0) model.
- Explain why this model was chosen using the ACF and PACF.
checkresiduals(bicoal)## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
pacf(bicoal)
There is significant lag up to lag level 4, and no significant lag after. Thus a (4, x, x) model is warranted.
- The last five values of the series are given below. The estimated parameters are
c=162.00, ϕ1=0.83,ϕ2=−0.34, ϕ3=0.55, and ϕ4=−0.38. Without using the forecast function, calculate forecasts for the next three years (1969–1971).
c=162
phi1=0.83
phi2=-0.34
phi3=0.55
phi4=-0.38
bicoal1969 <- c + phi1*545 + phi2*552 + phi3*534 + phi4*512
bicoal1970 <- c + phi1*bicoal1969 + phi2*545 + phi3*552 + phi4*534
bicoal1971 <- c + phi1*bicoal1970 + phi2*bicoal1969 + phi3*545 + phi4*552
bicoal1969## [1] 525.81bicoal1970## [1] 513.8023bicoal1971## [1] 499.6705e.Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?
bicoal.fc<-forecast(Arima(bicoal, order=c(4,0,0)), h=3)
bicoal.fc## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1969 527.6291 459.8804 595.3779 424.0164 631.2419
## 1970 517.1923 429.0014 605.3832 382.3160 652.0686
## 1971 503.8051 412.4786 595.1315 364.1334 643.4768The values I calculated are off by 2-4 under. This likely is a result of phi variables being rounded to two decimal points.
Exercise 18
#install.packages("Quandl")
library(Quandl)## Warning: package 'Quandl' was built under R version 4.0.5## Loading required package: xts## Warning: package 'xts' was built under R version 4.0.5## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.0.5##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, lastcarbon <-Quandl("BP/C02_EMMISSIONS_EUR")
summary(carbon)## Date Value
## Min. :1965-12-31 Min. :3428
## 1st Qu.:1979-09-30 1st Qu.:4396
## Median :1993-07-01 Median :4722
## Mean :1993-07-01 Mean :4639
## 3rd Qu.:2007-04-01 3rd Qu.:4919
## Max. :2020-12-31 Max. :5510I am using the “Carbon Dioxide (CO2) Emmissions - Total Europe” dataset
head(carbon)## Date Value
## 1 2020-12-31 3592.885
## 2 2019-12-31 4087.100
## 3 2018-12-31 4246.544
## 4 2017-12-31 4295.847
## 5 2016-12-31 4255.108
## 6 2015-12-31 4206.920- Plot graphs of the data, and try to identify an appropriate ARIMA model.
carbon$Date <-rev(carbon$Date)
carbon$Value<-rev(carbon$Value)
tail(carbon)## Date Value
## 51 2015-12-31 4206.920
## 52 2016-12-31 4255.108
## 53 2017-12-31 4295.847
## 54 2018-12-31 4246.544
## 55 2019-12-31 4087.100
## 56 2020-12-31 3592.885carbon.ts <- ts(carbon$Value, frequency = 1, start = c(1965, 1))
autoplot(carbon.ts)
auto.arima(carbon.ts)## Series: carbon.ts
## ARIMA(0,2,1)
##
## Coefficients:
## ma1
## -0.8881
## s.e. 0.0705
##
## sigma^2 estimated as 35922: log likelihood=-360.1
## AIC=724.2 AICc=724.44 BIC=728.18checkresiduals(carbon.ts)## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
pacf(carbon.ts)
ndiffs(carbon.ts)## [1] 2The graphs and ndiffs() function support the (0,2,1) model that auto.arima selected.
- Do residual diagnostic checking of your ARIMA model. Are the residuals white noise?
carbon.arima<- auto.arima(carbon.ts)
checkresiduals(carbon.arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,1)
## Q* = 7.496, df = 9, p-value = 0.5856
##
## Model df: 1. Total lags used: 10The residuals appear to be white noise per the Ljung-Box test and the ACF plot.
- Use your chosen ARIMA model to forecast the next four years.
carbon.arima.fc<- forecast(carbon.arima, h=4)
autoplot(carbon.arima.fc)
e. Now try to identify an appropriate ETS model.
carbon.ets<-ets(carbon.ts)
carbon.ets## ETS(A,N,N)
##
## Call:
## ets(y = carbon.ts)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 3427.7308
##
## sigma: 190.3948
##
## AIC AICc BIC
## 817.2823 817.7438 823.3583- Do residual diagnostic checking of your ETS model. Are the residuals white noise?
checkresiduals(carbon.ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 5.1865, df = 8, p-value = 0.7375
##
## Model df: 2. Total lags used: 10The Ljung-Box test and the ACF plot both indicate that the residuals are white noise in this model as well.
e.Use your chosen ETS model to forecast the next four years.
carbon.ets.fc<-forecast(carbon.ets, h=4)
autoplot(carbon.ets.fc)
h. Which of the two models do you prefer?
autoplot(carbon.arima.fc) + autolayer(carbon.ets.fc, series = "ETS", alpha=0.5)
I prefer the ARIMA model for this dataset as the ETS model appears to be similar to a naive model in this example. The yearly data means that there could be no seasonal component, and thus the ETS is somewhat subpar. The downward forecast for the ARIMA model also seems more appropriate as more European countries strive to meet their goals set at the Paris Climate Agreement.