ADEC7460 Homework #2
Jordan Harrop
Chapter Seven Excersises
7.1 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 (data set books).
head(books)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## Paperback Hardcover
## 1 199 139
## 2 172 128
## 3 111 172
## 4 209 139
## 5 161 191
## 6 119 1687.1a Plot the series and discuss the main features of the data.
plot(books, main = "Daily Sales or Paperback & Hardcover Books")
From the plots, both Hardcover and Paperback sales seem to be stationary with no clear trend or seasonality - This is required for simple exportential smoothing.
7.1b Use simple exponential smoothing with the ses function (setting initial=“simple”) and explore different values of alpha for the paperback series. Record the within-sample SSE for the one-step forecasts. Plot SSE against alpha and find which value of alpha works best. What is the effect of alpha on the forecasts?
ses.auto <- function(x, hz){
require(fma)
sq <- seq(0.1, 0.99, length.out = 18)
SSE <- c()
Alpha <- c()
for (i in sq) {
fit <- ses(x, alpha = i, intial="simple", h = hz)
error <- sum(residuals(fit)^2)
Alpha <- append(Alpha, i)
SSE <- append(SSE, error) }
df <- data.frame(Alpha = Alpha, SSE = SSE)
return(df)
}
y <- books[,1]
alpha.test <- ses.auto(y, 3)
plot(y = alpha.test$SSE, x = alpha.test$Alpha, ylab="SSE", xlab="Alpha", type = "o")
In this situation, as alpha increases the SSE also increases. However, this is only true once alpha reaches 0.16. From the graph, the optimal alpha should be between 0.13 and 0.16. (Assuming SSE was calculated correctly)
7.1c Now let ses select the optimal value of alpha. Use this value to generate forecasts for the next four days. Compare your results with 2.
fit1 <- ses(books[,1], h=4)
summary(fit1)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, 1], h = 4)
##
## Smoothing parameters:
## alpha = 0.1685
##
## Initial states:
## l = 170.8257
##
## sigma: 33.6377
##
## AIC AICc BIC
## 318.9747 319.8978 323.1783
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 7.176212 33.63769 27.8431 0.4737524 15.57782 0.7021303
## ACF1
## Training set -0.2117579
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 207.1098 164.0013 250.2182 141.1811 273.0384
## 32 207.1098 163.3934 250.8261 140.2513 273.9682
## 33 207.1098 162.7937 251.4258 139.3342 274.8853
## 34 207.1098 162.2021 252.0174 138.4294 275.7901 The ses function selected roughly the same alpha that was determined in step two. The only differece is the interval between tested alpha’s was smaller than what was used in the function used above.
7.1d Repeat but with initial=“optimal”. How much difference does an optimal initial level make?
fit2 <- ses(books[,1], h=4, intial="optimal")
summary(fit2)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, 1], h = 4, intial = "optimal")
##
## Smoothing parameters:
## alpha = 0.1685
##
## Initial states:
## l = 170.8257
##
## sigma: 33.6377
##
## AIC AICc BIC
## 318.9747 319.8978 323.1783
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 7.176212 33.63769 27.8431 0.4737524 15.57782 0.7021303
## ACF1
## Training set -0.2117579
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 207.1098 164.0013 250.2182 141.1811 273.0384
## 32 207.1098 163.3934 250.8261 140.2513 273.9682
## 33 207.1098 162.7937 251.4258 139.3342 274.8853
## 34 207.1098 162.2021 252.0174 138.4294 275.7901The `initial=“optimal” had no effect on the outcome of the model.
7.1e Repeat steps (b)-(d) with the hardcover series.
y <- books[,2]
alpha.test <- ses.auto(y, 3)
plot(y = alpha.test$SSE, x = alpha.test$Alpha, ylab="SSE", xlab="Alpha", type = "o")
fit3 <- ses(books[,2], h=4)
summary(fit3)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, 2], h = 4)
##
## Smoothing parameters:
## alpha = 0.3283
##
## Initial states:
## l = 149.2836
##
## sigma: 31.931
##
## AIC AICc BIC
## 315.8506 316.7737 320.0542
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 9.166918 31.93101 26.7731 2.636328 13.39479 0.7987858
## ACF1
## Training set -0.1417817
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 239.5602 198.6390 280.4815 176.9766 302.1439
## 32 239.5602 196.4905 282.6299 173.6908 305.4297
## 33 239.5602 194.4443 284.6762 170.5613 308.5591
## 34 239.5602 192.4869 286.6336 167.5677 311.5527fit4 <- ses(books[,2], h=4, intial="optimal")
summary(fit4)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, 2], h = 4, intial = "optimal")
##
## Smoothing parameters:
## alpha = 0.3283
##
## Initial states:
## l = 149.2836
##
## sigma: 31.931
##
## AIC AICc BIC
## 315.8506 316.7737 320.0542
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 9.166918 31.93101 26.7731 2.636328 13.39479 0.7987858
## ACF1
## Training set -0.1417817
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 239.5602 198.6390 280.4815 176.9766 302.1439
## 32 239.5602 196.4905 282.6299 173.6908 305.4297
## 33 239.5602 194.4443 284.6762 170.5613 308.5591
## 34 239.5602 192.4869 286.6336 167.5677 311.55277.2 Apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
7.2a Compare the SSE measures of Holt’s method for the two series to those of simple exponential smoothing in the previous question. Discuss the merits of the two forecasting methods for these data sets.
# SES Paperbacks - Initial = Simple
sum(residuals(fit1)^2)## [1] 33944.82# SES Hardcovers - Initial = Simple
sum(residuals(fit3)^2)## [1] 30587.69# Holt Paperbacks - Initial = Simple
sum(residuals(fit5)^2)## [1] 46917.39# Holt Hardcovers - Initial = Simple
sum(residuals(fit6)^2)## [1] 36842.1 Both the Hardcover and the Paperback time series do not show a strong trend - subjectivly, the Hardcover time series shows the strongest growth trend of the the two series. Both ses models preformed better than the Holt trend method. It is interesting to see that difference between the simple exponential smoothing SSE and the Holt trend method SSE for the Hardcover series was not as great as the Paperback - maybe because Hardcover showed more of a trend.
7.2b Compare the forecasts for the two series using both methods. Which do you think is best?
fcst.p.ses <- forecast(fit1, h=4); accuracy(fcst.p.ses)## ME RMSE MAE MPE MAPE MASE
## Training set 7.176212 33.63769 27.8431 0.4737524 15.57782 0.7021303
## ACF1
## Training set -0.2117579fcst.h.ses <- forecast(fit3, h=4); accuracy(fcst.h.ses)## ME RMSE MAE MPE MAPE MASE
## Training set 9.166918 31.93101 26.7731 2.636328 13.39479 0.7987858
## ACF1
## Training set -0.1417817fcst.p.holt <- forecast(fit5, h=4); accuracy(fcst.p.holt)## ME RMSE MAE MPE MAPE MASE
## Training set 7.769844 39.54634 33.5377 1.633306 18.19621 0.8457332
## ACF1
## Training set -0.1088681fcst.h.holt <- forecast(fit6, h=4); accuracy(fcst.h.holt)## ME RMSE MAE MPE MAPE MASE
## Training set 7.193267 35.04383 27.99174 2.423793 14.18241 0.8351445
## ACF1
## Training set -0.07743714Personally, I think the ses method seems to work best for this data - particuraly because of the lack of strong trend present in the data. The simple expontential models have better RMSE and MAE.
7.2c Calculate a 95% prediction interval for the first forecast for each series using both methods, assuming normal errors. Compare your forecasts with those produced by R.
# The formula for PI = point forecast +/- 1.96*Sigma
sigma <- sqrt(fit1$model$sigma2)
adjst <-1.96*sigma
fcst.p.ses## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 207.1098 164.0013 250.2182 141.1811 273.0384
## 32 207.1098 163.3934 250.8261 140.2513 273.9682
## 33 207.1098 162.7937 251.4258 139.3342 274.8853
## 34 207.1098 162.2021 252.0174 138.4294 275.7901a <- 207.1098 + adjst
b <- 207.1098 - adjst
results <- c(b,a); results## [1] 141.1799 273.0397sigma <- sqrt(fit3$model$sigma2)
adjst <-1.96*sigma
fcst.h.ses## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 239.5602 198.6390 280.4815 176.9766 302.1439
## 32 239.5602 196.4905 282.6299 173.6908 305.4297
## 33 239.5602 194.4443 284.6762 170.5613 308.5591
## 34 239.5602 192.4869 286.6336 167.5677 311.5527a <- 239.5602 + adjst
b <- 239.5602 - adjst
results <- c(b,a); results## [1] 176.9754 302.1450sigma <- sqrt(fit5$model$sigma2)
adjst <-1.96*sigma
fcst.p.holt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 222.0201 171.3394 272.7007 144.51068 299.5295
## 32 229.6904 164.8872 294.4935 130.58245 328.7983
## 33 237.3606 145.1175 329.6038 96.28696 378.4343
## 34 245.0309 115.5211 374.5407 46.96280 443.0991a <- 222.0201 + adjst
b <- 222.0201 - adjst
results <- c(b,a); results## [1] 144.5093 299.5309sigma <- sqrt(fit6$model$sigma2)
adjst <-1.96*sigma
fcst.h.holt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 250.7889 205.8784 295.6993 182.1042 319.4735
## 32 254.7003 202.4087 306.9918 174.7273 334.6733
## 33 258.6117 196.3181 320.9052 163.3419 353.8815
## 34 262.5231 187.9903 337.0558 148.5350 376.5111a <- 250.7889 + adjst
b <- 250.7889 - adjst
results <- c(b,a); results## [1] 182.1030 319.4748The calculated prediction interval and the R output prediction interval are the same. It is important to note that the prediction interval grows as the the point forecasts moove further into the future - this is due to the alpha2 being mutiplied by each step h.
7.3 For this exercise, use the quarterly UK passenger vehicle production data from 1977:1–2005:1 (data set ukcars).
7.3a Plot the data and describe the main features of the series.
data(ukcars)
plot(ukcars)
For this time series, there looks to be seasonality, cyclical behavior and trend.
7.3b Decompose the series using STL and obtain the seasonally adjusted data.
decomp.fit <- stl(ukcars, s.window = "periodic", robust = TRUE)
plot(decomp.fit)
7.3c Forecast the next two years of the series using an additive damped trend method applied to the seasonally adjusted data. Then reseasonalize the forecasts. Record the parameters of the method and report the RMSE of the one-step forecasts from your method.
ukcars.sadjst <- seasadj(decomp.fit)
ukcars.stl.fit <- stlf(ukcars.sadjst, etsmodel="AAN", damped=TRUE)
ukcars.stl.fcst <- forecast(ukcars.stl.fit, h=2); ukcars.stl.fcst## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 391.9193 364.4475 419.3912 349.9047 433.9339
## 2005 Q3 418.4930 383.2300 453.7559 364.5629 472.4230accuracy(ukcars.stl.fit)## ME RMSE MAE MPE MAPE MASE
## Training set 1.368317 21.43641 16.67588 0.09849256 5.394512 0.5434588
## ACF1
## Training set 0.0090413477.3d Forecast the next two years of the series using Holt’s linear method applied to the seasonally adjusted data. Then reseasonalize the forecasts. Record the parameters of the method and report the RMSE of of the one-step forecasts from your method.
ukcars.hw.fit <- hw(ukcars.sadjst, damped=FALSE)
ukcars.hw.fcst <- forecast(ukcars.hw.fit, h=2); ukcars.hw.fcst## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 404.3676 371.5179 437.2173 354.1283 454.6068
## 2005 Q3 408.1908 369.6463 446.7352 349.2422 467.1394accuracy(ukcars.hw.fit)## ME RMSE MAE MPE MAPE MASE
## Training set 2.03986 25.63276 20.80823 0.3334681 6.710032 0.67813
## ACF1
## Training set 0.032014687.3e Now use ets() to choose a seasonal model for the data.
ukcars.ets.fit <- ets(ukcars)
summary(ukcars.ets.fit)## ETS(A,N,A)
##
## Call:
## ets(y = ukcars)
##
## Smoothing parameters:
## alpha = 0.6267
## gamma = 0.0001
##
## Initial states:
## l = 313.0916
## s=-1.1271 -44.606 21.5553 24.1778
##
## sigma: 25.2579
##
## AIC AICc BIC
## 1277.980 1279.047 1297.072
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1.324962 25.25792 20.12508 -0.1634983 6.609629 0.6558666
## ACF1
## Training set 0.01909295accuracy(ukcars.ets.fit)## ME RMSE MAE MPE MAPE MASE
## Training set 1.324962 25.25792 20.12508 -0.1634983 6.609629 0.6558666
## ACF1
## Training set 0.019092957.3f Compare the RMSE of the fitted model with the RMSE of the model you obtained using an STL decomposition with Holt’s method. Which gives the better in-sample fits?
After comparing both method’s accuracy metrics, the STL decomposition method fits the data best.
7.3g Compare the forecasts from the two approaches? Which seems most reasonable?
plot(ukcars.stl.fcst)
plot(ukcars.hw.fcst)
From looking at the plots, both forecasts look reasonable. If I had to pick one, the STL + ETS method would be the one I would go with.
7.4 For this exercise, use the monthly Australian short-term overseas visitors data, May 1985–April 2005. (Data set: visitors.)
7.4a Make a time plot of your data and describe the main features of the series.
data("visitors")
plot(visitors)
The series shows multiplicative seasonality with an upward trend.
7.4b Forecast the next two years using Holt-Winters’ multiplicative method.
visitors.hw.fit <- hw(visitors, seasonal = "multiplicative")
visitors.hw.fcst <- forecast(visitors.hw.fit, h=2); visitors.hw.fcst## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 369.3175 343.3002 395.3348 329.5275 409.1076
## Jun 2005 395.5080 365.2767 425.7393 349.2733 441.74277.4c Why is multiplicative seasonality necessary here?
The sesonal peaks and valleys are growing as the series growth trends upward.
7.4d Experiment with making the trend exponential and/or damped.
visitors.hw.fit.exp <- hw(visitors, seasonal = "multiplicative", exponential = TRUE)
visitors.hw.fcst.exp <- forecast(visitors.hw.fit.exp, h=2)
visitors.hw.fcst.exp## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 363.2966 338.4841 390.1402 324.3823 402.6815
## Jun 2005 391.2661 359.9892 424.2877 343.7333 441.0104visitors.hw.fit.exp.d <- hw(visitors, seasonal = "multiplicative", exponential = TRUE, damped = TRUE)
visitors.hw.fcst.exp.d <- forecast(visitors.hw.fit.exp.d, h=2)
visitors.hw.fcst.exp.d## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 355.2463 329.9847 380.5688 315.5618 393.091
## Jun 2005 382.6635 351.8544 415.0455 335.0439 431.0917.4e Compare the RMSE of the one-step forecasts from the various methods. Which do you prefer?
accuracy(visitors.hw.fit)## ME RMSE MAE MPE MAPE MASE
## Training set -0.9498442 14.8295 10.96716 -0.8150922 4.271167 0.4050069
## ACF1
## Training set 0.2223887accuracy(visitors.hw.fit.exp)## ME RMSE MAE MPE MAPE MASE
## Training set 0.6442177 14.49416 10.62951 0.2554469 4.0328 0.3925378
## ACF1
## Training set 0.07595792accuracy(visitors.hw.fit.exp.d)## ME RMSE MAE MPE MAPE MASE
## Training set 0.7230142 14.45533 10.72791 0.03798703 4.090931 0.3961716
## ACF1
## Training set 0.01218167From comparing the metrics, the exponential and the exponential damped models are better than just the multiplicative hw model. However, both the exponential and damped models are close to each other. Picking one would depend on how long the forecasts need to be - for shot-term forecasts the exponential only model may be best.
7.4f Now fit each of the following models to the same data: a multiplicative Holt-Winters’ method; an ETS model; an additive ETS model applied to a Box-Cox transformed series; a seasonal naive method applied to the Box-Cox transformed series; an STL decomposition applied to the Box-Cox transformed data followed by an ETS model applied to the seasonally adjusted (transformed) data.
visitors.ets.fit <- ets(visitors)
lam <- BoxCox.lambda(visitors)
visitors.ets2.fit <- ets(visitors, additive.only = TRUE, lambda = lam)
visitors.naive.fit <- snaive(visitors, lambda = lam)
visitors.boxcox <- BoxCox(visitors, lambda = lam)
decomp.fit <- stl(visitors.boxcox, s.window = "periodic", robust = TRUE)
visitors.sadjst <- seasadj(decomp.fit)
visitors.stl.fit <- stlf(visitors.sadjst, etsmodel="AAN", damped=TRUE)7.4g For each model, look at the residual diagnostics and compare the forecasts for the next two years. Which do you prefer?
forecast(visitors.hw.fit, h=2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 369.3175 343.3002 395.3348 329.5275 409.1076
## Jun 2005 395.5080 365.2767 425.7393 349.2733 441.7427forecast(visitors.ets.fit, h=2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 360.8599 337.0861 384.6337 324.501 397.2187
## Jun 2005 395.8905 364.7028 427.0782 348.193 443.5880forecast(visitors.ets2.fit, h=2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 364.4376 343.0577 386.7638 332.1144 398.9744
## Jun 2005 399.2435 372.0686 427.8242 358.2372 443.5378forecast(visitors.naive.fit, h=2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 341.3 292.9977 395.1076 269.5329 425.9462
## Jun 2005 367.3 316.3105 423.9776 291.4957 456.4110forecast(visitors.stl.fit, h=2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 15.74975 15.50668 15.99282 15.37800 16.1215
## Jun 2005 15.99054 15.69292 16.28815 15.53537 16.4457plot(visitors.hw.fit$residuals)
plot(visitors.ets.fit$residuals)
plot(visitors.ets2.fit$residuals)
plot(visitors.naive.fit$residuals)
plot(visitors.stl.fit$residuals)
From looking at the residual plots, the BOX-COX transformed ETS model seems to be best representing white noise.
Chapter Eight Excersises
8.1 Figure 8.24 shows the ACFs for 36 random numbers, 360 random numbers and for 1,000 random numbers.
8.1a Explain the differences among these figures. Do they all indicate the data are white noise?
Yes, they do all indicate that the data is white noise. The chart to the left indicates that there is some autocorrelations starting at 5 to 7. The middle chart shows more autocorrelation between 5 and 10, but none of the observations break the critical barrier - there isn’t any significance. The right chart shows mostly all white noise.
8.1b 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 critical values are determined by the data and it’s a reasonable assumption to think that some observatins will break the critical values while still being white noise.
8.2 A classic example of a non-stationary series is the daily closing IBM stock prices (data set ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows the series is non-stationary and should be differenced.
data("ibmclose")
tsdisplay(ibmclose)
The ACF plot shows strong autocorrelation out to a 25 obseration lag, however, with differencing the PACF plot shows that the autocorrelation can be corrected.
8.3 For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
# Custom Function
transform.plot <- function(x){
lam <- BoxCox.lambda(x)
transformed <- BoxCox(x, lambda = lam)
tsdisplay(transformed)
return(paste("Box-Cox Lambda = ", round(lam, 4)))
}8.3a usnetelec
data("usnetelec")
transform.plot(usnetelec)
## [1] "Box-Cox Lambda = 0.5168"The series needs one differencing in order to be stationary
8.3b usgdp
data("usgdp")
transform.plot(usgdp)
## [1] "Box-Cox Lambda = 0.3664"The series needs one differencing in order to be stationary
8.3c mcopper
data("mcopper")
transform.plot(mcopper)
## [1] "Box-Cox Lambda = 0.1919"The series needs two steps of differencing in order to be stationary
8.3d enplanements
data("enplanements")
transform.plot(enplanements)
## [1] "Box-Cox Lambda = -0.2269"The series needs one differencing and seasonal differencing in order to be stationary
8.3e visitors
data("visitors")
transform.plot(visitors)
## [1] "Box-Cox Lambda = 0.2775"The series needs one differencing and seasonal differencing in order to be stationary
8.4 For the enplanements data, write down the differences you chose above using backshift operator notation.
Since we are taking seasonal differencing, after transformation the backshift notation should be y(t) = (1-B^12) * y(t)
8.5 Use R to simulate and plot some data from simple ARIMA models.
8.5a Use the following R code to generate data from an AR(1) model with beta(1) = 0.6 and sigma2 = 1. The process starts with y0 = 0.
set.seed(1234)
y1 <- ts(numeric(100))
y2 <- ts(numeric(100))
y3 <- ts(numeric(100))
y4 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100) {
y1[i] <- 0.1*y1[i-1] + e[i]
y2[i] <- 0.6*y2[i-1] + e[i]
y3[i] <- 1*y3[i-1] + e[i]
y4[i] <- 1.6*y4[i-1] + e[i]
}8.5b Produce a time plot for the series. How does the plot change as you change beta(1)?
tsdisplay(y1)
tsdisplay(y2)
tsdisplay(y3)
tsdisplay(y4)
As beta(1) increase the more autocorrelation is present in the series. The first plot show white noise and as beta(1) is increased more autocorrelation is present.
8.5c Write your own code to generate data from an MA(1) model with beta(1) = 0.6 and sigma2 = 1.
y1 <- ts(numeric(100))
y2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y1[i] <- 3 + e[i] + 0.6*e[i-1]
}
for(i in 2:100){
y2[i] <- 3 + e[i] + 1.1*e[i-1]
}8.5d Produce a time plot for the series. How does the plot change as you change beta(1)
tsdisplay(y1)
tsdisplay(y2)
In theory, the series should remain similiar with beta(1) changing.
8.5e Generate data from an ARMA(1,1) model.
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <- 0.6*y[i-1] + e[i] + 0.6*e[i-1]
}
tsdisplay(y)
8.5f Generate data from an AR(2) model.
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 3:100){
y[i] <- -0.8*y[i-1]+ 0.3*y[i-2] + e[i]
}
tsdisplay(y)
8.5g Graph the latter two series and compare them.
The ARMA(1) model looks like it may have seasonality to it, but overall, stationary. The AR(2) model is unique with a multiplicative (osculating) pattern happening at the tail.
8.6 Consider the number of women murdered each year (per 100,000 standard population) in the United States (data set wmurders).
data("wmurders")8.6a By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data.
transform.plot(wmurders)
## [1] "Box-Cox Lambda = -0.0953"wmurders.diff <- diff(wmurders)
tsdisplay(wmurders.diff)
Maybe a (2,1,0) ARIMA will work best.
8.6b Should you include a constant in the model? Explain.
No, the plot above does not show the average drastically changing.
8.6c Write this model in terms of the backshift operator.
(1-B-B2)(1-B)yt
8.6d Fit the model using R and examine the residuals. Is the model satisfactory?
fit <- Arima(wmurders, order = c(2,1,0))
plot(fit$residuals)
The models is satisfactory becasue there is no clear pattern in the resdiuals.
8.6e Forecast three times ahead. Check your forecasts by hand to make sure you know how they have been calculated.
fcst <- forecast(fit, h=3)
fcst## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 2.553355 2.288694 2.818015 2.148591 2.958118
## 2006 2.533802 2.170056 2.897548 1.977501 3.090104
## 2007 2.524231 2.033825 3.014637 1.774220 3.2742428.6f Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
plot(fcst)
8.6g Does auto.arima give the same model you have chosen? If not, which model do you think is better?
fit2 <- auto.arima(wmurders)
fit2## Series: wmurders
## ARIMA(1,2,1)
##
## Coefficients:
## ar1 ma1
## -0.2434 -0.8261
## s.e. 0.1553 0.1143
##
## sigma^2 estimated as 0.04632: log likelihood=6.44
## AIC=-6.88 AICc=-6.39 BIC=-0.97The models are different - Given my confidence in my ability to chose a correct model, I would assume the auto.arima is a better model.
8.7 Consider the quarterly number of international tourists to Australia for the period 1999-2010. (Data set austourists.)
data("austourists")8.7a Describe the time plot.
tsdisplay(austourists)
The time series looks to have multplicative seasonality with a growth trend.
8.7b What can you learn from the ACF graph?
There is some autocorellation in the lagged observations.
8.7c What can you learn from the PACF graph?
Seasonal differencing is required for this data.
8.7d Produce plots of the seasonally differenced data (1-B4)Yt. What model do these graphs suggest?
tourists.decom <- decompose(austourists, "multiplicative")
austourists.adjst <- austourists / tourists.decom$seasonal
austourists.adjst.diff <- diff(austourists.adjst)
tsdisplay(austourists.adjst.diff)
Maybe (1,0,0)?
8.7e Does auto.arima give the same model that you chose? If not, which model do you think is better?
fit <- auto.arima(austourists)
fit## Series: austourists
## ARIMA(1,0,0)(1,1,0)[4] with drift
##
## Coefficients:
## ar1 sar1 drift
## 0.4493 -0.5012 0.4665
## s.e. 0.1368 0.1293 0.1055
##
## sigma^2 estimated as 5.606: log likelihood=-99.47
## AIC=206.95 AICc=207.97 BIC=214.09Yes. But, I would still go with the auto.arima.
8.8 Consider the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period 1985-1996). (Data set usmelec.) In general there are two peaks per year: in mid-summer and mid-winter.
data("usmelec")
tsdisplay(usmelec)
8.8a Examine the 12-month moving average of this series to see what kind of trend is involved.
usmelec.ma <- ma(usmelec, order = 12)
plot(usmelec.ma)
8.8b Do the data need transforming? If so, find a suitable transformation.
usmelec.lam <- BoxCox.lambda(usmelec)
usmelec.boxcox <- BoxCox(usmelec, usmelec.lam)
tsdisplay(usmelec.boxcox)
8.8c Are the data stationary? If not, find an appropriate differencing which yields stationary data.
usmelec.bcx.diff <- diff(usmelec.boxcox)
usmelec.bcx.decom <- decompose(usmelec.bcx.diff, "multiplicative")
usmelec.adjst <- usmelec.bcx.diff / usmelec.bcx.decom$seasonal
tsdisplay(usmelec.adjst)
8.8d 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?
fit <- auto.arima(usmelec, trace = TRUE, ic ="aic", lambda = usmelec.lam)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0) with drift : -3324.434
## ARIMA(1,1,0)(1,0,0)[12] with drift : Inf
## ARIMA(0,1,1)(0,0,1)[12] with drift : -3626.366
## ARIMA(0,1,0) : -3326.287
## ARIMA(0,1,1)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,1) with drift : -3324.355
## ARIMA(0,1,1)(0,0,2)[12] with drift : -3783.239
## ARIMA(1,1,1)(0,0,2)[12] with drift : -3836.101
## ARIMA(1,1,0)(0,0,2)[12] with drift : -3776.41
## ARIMA(1,1,2)(0,0,2)[12] with drift : -3873.503
## ARIMA(2,1,3)(0,0,2)[12] with drift : -3906.568
## ARIMA(2,1,3)(0,0,2)[12] : -3897.851
## ARIMA(2,1,3)(1,0,2)[12] with drift : Inf
## ARIMA(2,1,3)(0,0,1)[12] with drift : -3809.203
## ARIMA(1,1,3)(0,0,2)[12] with drift : -3886.669
## ARIMA(3,1,3)(0,0,2)[12] with drift : Inf
## ARIMA(2,1,2)(0,0,2)[12] with drift : -3898.42
## ARIMA(2,1,4)(0,0,2)[12] with drift : -3903.185
## ARIMA(3,1,4)(0,0,2)[12] with drift : Inf
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(2,1,3)(0,0,2)[12] with drift : Inf
## ARIMA(2,1,4)(0,0,2)[12] with drift : Inf
## ARIMA(2,1,2)(0,0,2)[12] with drift : -3916.965
##
## Best model: ARIMA(2,1,2)(0,0,2)[12] with driftfit## Series: usmelec
## ARIMA(2,1,2)(0,0,2)[12] with drift
## Box Cox transformation: lambda= -0.4772402
##
## Coefficients:
## ar1 ar2 ma1 ma2 sma1 sma2 drift
## 0.849 -0.3689 -1.2181 0.2639 0.7653 0.4396 0.0001
## s.e. 0.114 0.0728 0.1135 0.1109 0.0529 0.0360 0.0001
##
## sigma^2 estimated as 0.000009869: log likelihood=1966.48
## AIC=-3916.96 AICc=-3916.64 BIC=-3884.04The best model is ARIMA(2,1,2)
8.8e 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.
plot(fit$residuals)
8.8f Forecast the next 15 years of generation of electricity by the U.S. electric industry. Get the latest figures from http://data.is/zgRWCO to check on the accuracy of your forecasts.
hh <- 12*15
fcst <- forecast(fit, h=hh)
plot(fcst)
8.8g How many years of forecasts do you think are sufficiently accurate to be usable?
Probably only a year to two years ahead will be useful.
8.9 For the mcopper data:
data(mcopper)
plot(mcopper)
8.9a If necessary, find a suitable Box-Cox transformation for the data;
lam <- BoxCox.lambda(mcopper)
mcopper.bcx <- BoxCox(mcopper, lambda = lam)
tsdisplay(mcopper.bcx)
8.9b Fit a suitable ARIMA model to the transformed data using auto.arima();
fit <- auto.arima(mcopper, trace = TRUE, ic ="aic", lambda = lam)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2)(1,0,1)[12] with drift : -66.32021
## ARIMA(0,1,0) with drift : -13.82248
## ARIMA(1,1,0)(1,0,0)[12] with drift : -62.06849
## ARIMA(0,1,1)(0,0,1)[12] with drift : -82.25982
## ARIMA(0,1,0) : -12.82906
## ARIMA(0,1,1)(1,0,1)[12] with drift : -71.92506
## ARIMA(0,1,1) with drift : -83.13952
## ARIMA(1,1,1) with drift : -80.30584
## ARIMA(0,1,2) with drift : -81.16609
## ARIMA(1,1,2) with drift : -78.27629
## ARIMA(0,1,1) : -83.33279
## ARIMA(0,1,1)(1,0,0)[12] : -74.05872
## ARIMA(0,1,1)(0,0,1)[12] : -82.603
## ARIMA(0,1,1)(1,0,1)[12] : -72.07316
## ARIMA(1,1,1) : -80.55931
## ARIMA(0,1,2) : -81.33929
## ARIMA(1,1,2) : -80.62448
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(0,1,1) : -86.0969
##
## Best model: ARIMA(0,1,1)fit## 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.438.9c Try some other plausible models by experimenting with the orders chosen;
fit2 <- Arima(mcopper, order = c(0,1,2), lambda = lam)
fit3 <- Arima(mcopper, order = c(1,1,2), lambda = lam)
accuracy(fit)## ME RMSE MAE MPE MAPE MASE
## Training set 3.480533 77.27254 44.92858 0.166202 4.303677 0.2021433
## ACF1
## Training set -0.08442198accuracy(fit2)## ME RMSE MAE MPE MAPE MASE
## Training set 3.505681 77.25525 44.94118 0.167361 4.303859 0.2022
## ACF1
## Training set -0.08273006accuracy(fit3)## ME RMSE MAE MPE MAPE MASE
## Training set 3.392046 77.82026 44.9212 0.1627576 4.298004 0.2021101
## ACF1
## Training set -0.10688678.9d Choose what you think is the best model and check the residual diagnostics;
#Fit - From auto.arima
plot(fit$residuals)
8.9e Produce forecasts of your fitted model. Do the forecasts look reasonable?
fcst <- forecast(fit, h=10)
plot(fcst)
Yes, the forecasts do look reasonable given the nature of the series.
8.9f Compare the results with what you would obtain using ets() (with no transformation).
fit4 <- ets(mcopper); fit4## ETS(M,Ad,N)
##
## Call:
## ets(y = mcopper)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.2174
## phi = 0.8
##
## Initial states:
## l = 265.142
## b = -3.4724
##
## sigma: 0.0629
##
## AIC AICc BIC
## 8052.055 8052.206 8078.065fcst2 <- forecast(fit4, h=10)
plot(fcst2)
8.10 Choose one of the following seasonal time series: condmilk, hsales, uselec
data("condmilk")
data("hsales")
data("uselec")
tsdisplay(condmilk)
tsdisplay(hsales)
tsdisplay(uselec)
8.10a Do the data need transforming? If so, find a suitable transformation.
lam <- BoxCox.lambda(condmilk)
lam## [1] -0.4026742condmilk.bcx <- BoxCox(condmilk, lambda = lam)
tsdisplay(condmilk.bcx)
8.10b Are the data stationary? If not, find an appropriate differencing which yields stationary data.
tsdisplay(diff(diff(condmilk.bcx,4)))
8.10c 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?
fit <- Arima(condmilk.bcx, order=c(0,1,1), seasonal=c(0,1,1))
fit2 <- Arima(condmilk.bcx, order=c(0,1,2), seasonal=c(0,1,2))
tsdisplay(residuals(fit))
tsdisplay(residuals(fit2))
8.10d 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.
fit3 <- auto.arima(condmilk, trace = TRUE, ic ="aic", lambda = lam)##
## ARIMA(2,0,2)(1,0,1)[12] with non-zero mean : Inf
## ARIMA(0,0,0) with non-zero mean : -331.2422
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean : -502.0451
## ARIMA(0,0,1)(0,0,1)[12] with non-zero mean : -454.303
## ARIMA(0,0,0) with zero mean : 517.6466
## ARIMA(1,0,0) with non-zero mean : -451.2644
## ARIMA(1,0,0)(2,0,0)[12] with non-zero mean : -516.0346
## ARIMA(1,0,0)(2,0,1)[12] with non-zero mean : Inf
## ARIMA(0,0,0)(2,0,0)[12] with non-zero mean : Inf
## ARIMA(2,0,0)(2,0,0)[12] with non-zero mean : Inf
## ARIMA(1,0,1)(2,0,0)[12] with non-zero mean : -514.1255
## ARIMA(2,0,1)(2,0,0)[12] with non-zero mean : Inf
## ARIMA(1,0,0)(2,0,0)[12] with zero mean : Inf
##
## Best model: ARIMA(1,0,0)(2,0,0)[12] with non-zero meantsdisplay(residuals(fit3))
8.10e Forecast the next 24 months of data using your preferred model.
fcst <- forecast(fit3, h=10)
plot(fcst)
8.10f Compare the forecasts obtained using ets().
fit4 <- ets(condmilk)
fcst2 <- forecast(fit4, h=10)
plot(fcst2)
Both forecasts are very similiar. However, the ARIMA model seems to be more conservitive.
8.11 For the same time series you used in exercise Q10, 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 exercise Q10. Which do you think is the best approach?
condmilk.decomp <- stl(condmilk.bcx, t.window=13, s.window="periodic")
condmilk.sadjst <- seasadj(condmilk.decomp)
condmilk.stl <- stlf(condmilk.sadjst, method="arima")
fcst3 <- forecast(condmilk.stl, h=10)
plot(fcst3)
The best approach depends a lot on the situation, but I do like the stlf() approach.