Chapter 7- Exponential Smoothing
##Packages used
library(readr)
library(ggplot2)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(fpp2)
## Loading required package: fma
##
## Attaching package: 'fma'
## The following object is masked from 'package:GGally':
##
## pigs
## Loading required package: expsmooth
1. Consider the pigs series — the number of pigs slaughtered in Victoria each month.
## 1a. Use the ses() function in R to find the optimal values of α and ℓ0, and generate forecasts for the next four months.
str(pigs)
## Time-Series [1:188] from 1980 to 1996: 76378 71947 33873 96428 105084 ...
pigs.ses <- ses(pigs,h=4)
summary(pigs.ses)
##
## 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
So, α=0.2971 and l0= 77260.0561
## 1b. 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.
pt <- 98816.41
s <-sd(residuals(pigs.ses))
paste('Point forecast:', pt)
## [1] "Point forecast: 98816.41"
paste('Lo 95:',pt - 1.96 * s)
## [1] "Lo 95: 78679.9711418255"
paste('Hi 95:',pt + 1.96 * s)
## [1] "Hi 95: 118952.848858174"
It looks like the intervals are slightly larger. The lower 95% is higher than the one calculated by R, and the upper 95% is lower.
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()?
pigs.f <- ses(pigs,h=4, alpha = 0.1, initial = "simple")
pigs.f
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98389.03 84025.49 112752.6 76421.90 120356.2
## Oct 1995 98389.03 83953.85 112824.2 76312.34 120465.7
## Nov 1995 98389.03 83882.57 112895.5 76203.31 120574.7
## Dec 1995 98389.03 83811.63 112966.4 76094.83 120683.2
forecast(pigs.f,h=4)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98389.03 84025.49 112752.6 76421.90 120356.2
## Oct 1995 98389.03 83953.85 112824.2 76312.34 120465.7
## Nov 1995 98389.03 83882.57 112895.5 76203.31 120574.7
## Dec 1995 98389.03 83811.63 112966.4 76094.83 120683.2
I did not get the same forecast as ses(), the forecasts are slightly different.
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?
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.
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.
## 5a. Plot the series and discuss the main features of the data.
autoplot(books, ylab= "Book Sales", xlab="Days", main="Book Sales over Time")
We observe seasonal pattern in the sales of hardcover and paperback book sales. This simply means that as the sales for hardcover books increase, there is a decrease paperback books (vice versa).
## 5b. Use the ses() function to forecast each series, and plot the forecasts.
books <- data.frame(books)
Hardcover <- ts(books$Hardcover)
Paperback <- ts(books$Paperback)
##Predicting the next four days' sales of hardcover and paperback books
Hardcover.fcast <- ses(Hardcover, h=4)
Hardcover.fcast
## 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
Paperback.fcast <- ses(Paperback, h=4)
Paperback.fcast
## 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
##Plotting the forecasts
autoplot(Hardcover.fcast)
autoplot(Paperback.fcast)
## 5c.Compute the RMSE values for the training data in each case.
round(accuracy(Hardcover.fcast),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.17 31.93 26.77 2.64 13.39 0.8 -0.14
round(accuracy(Paperback.fcast),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.18 33.64 27.84 0.47 15.58 0.7 -0.21
Hardcover RSME = 31.93
Paperback RSME = 33.64
6.We will continue with the daily sales of paperback and hardcover books in data set books.
##6a. Apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
#Hardcover
Hardcover.holt <- holt(Hardcover,h=4)
Hardcover.holt
## 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.3334
#Paperback
Paperback.holt <- holt(Paperback, h=4)
Paperback.holt
## 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.7735
##6b. 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.
round(accuracy(Hardcover.holt),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.14 27.19 23.16 -2.11 12.16 0.69 -0.03
round(accuracy(Paperback.holt),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -3.72 31.14 26.18 -5.51 15.58 0.66 -0.18
In this case the RSME for Hardcover is 27.19 which compared to the previous RSME, it looks like using holt's method improved by 4.74. RSME for Paperback is 31.14, which improved by 2.5.
From what we've observed, Holt's method is much better for this dataset than the SES for this specific dataset as it shows an upward trend making the Holt's method the most appropriate method to be used here.
##6c. Compare the forecasts for the two series using both methods. Which do you think is best?
# Holt's method is much better because it takes into account the trend elements of all time series datasets while SES does not account for the tres. We also observed a much lower RSME in the holt's method comapred to the SES method.
##6d. 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.
#Upper - Hardcover
Hardcover.holt$mean[1]+round(accuracy(Hardcover.holt),2)*1.96
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 249.8995 303.4663 295.5675 246.0383 274.0075 251.5263 250.1151
#Lower - Hardcover
Hardcover.holt$mean[1]-round(accuracy(Hardcover.holt),2)*1.96
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 250.4483 196.8815 204.7803 254.3095 226.3403 248.8215 250.2327
#Upper - Paperback
Paperback.holt$mean[1]+round(accuracy(Paperback.holt),2)*1.96
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 202.1756 270.5012 260.7796 198.6672 240.0036 210.7604 209.114
#Lower - Paperback
Paperback.holt$mean[1]-round(accuracy(Paperback.holt),2)*1.96
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 216.758 148.4324 158.154 220.2664 178.93 208.1732 209.8196
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.
[Hint: use h=100 when calling holt() so you can clearly see the differences between the various options when plotting the forecasts.]
Which model gives the best RMSE?
default <- holt(eggs, h=100)
damped <- holt(eggs, h=100, damped = T)
exponential <- holt(eggs, h=100, exponential = T)
lambda <- holt(eggs, h=100, lambda = 'auto', biasadj = T)
da_ex <- holt(eggs, h=100, exponential = T, damped = T)
da_la <- holt(eggs, h=100, damped = T, lambda = 'auto', biasadj = T)
autoplot(eggs) +
autolayer(default, series='Default', PI=F) +
autolayer(damped, series='Damped', PI=F) +
autolayer(exponential, series='Exponential', PI=F) +
autolayer(lambda, series='Box-Cox Transformed', PI=F) +
autolayer(da_ex, series='Damped & Exponential', PI=F) +
autolayer(da_la, series='Damped & Box-Cox', PI=F) +
ggtitle('Forecast of US Eggs Prices') +
xlab('Year') +
ylab('Price of Dozen Eggs')
##Accuracy of the forecasts
accuracy(default)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## ACF1
## Training set 0.01348202
accuracy(damped)
## ME RMSE MAE MPE MAPE MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
## ACF1
## Training set -0.003195358
accuracy(exponential)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.4918791 26.49795 19.29399 -1.263235 9.766049 0.9517436 0.0103908
accuracy(lambda)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.2015298 26.38689 18.99362 -1.63043 9.713172 0.9369265 0.0383996
accuracy(da_ex)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.9089678 26.59113 19.54973 -2.125756 10.02328 0.964359
## ACF1
## Training set 0.01376115
accuracy(da_la)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.806213 26.58589 19.55896 -2.58425 10.1176 0.9648141 0.005322063
## From above, the box-cox forcast has the lowest RSME.
8.Recall your retail time series data (from Exercise 3 in Section 2.10).
##Import retail data
library(readxl)
retail <- read_excel("/Users/dorothymensah/Desktop/retail.xlsx")
## New names:
## * `` -> ...1
#8a. Why is multiplicative seasonality necessary for this series
retail.ts <- ts(retail[,20], frequency = 12, start = c(1982,4))
## Warning in data.matrix(data): NAs introduced by coercion
autoplot(retail.ts)
Multiplicative Seasonality is necessary for this series because as we can see from above, the dataset has an increasing trend.
##8b.Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
#damped=F
holt.fit1 <- hw(retail.ts,seasonal = "multiplicative",damped = FALSE)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Missing values encountered. Using longest contiguous portion of time series
summary(holt.fit1)
##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = retail.ts, seasonal = "multiplicative", damped = FALSE)
##
## Smoothing parameters:
## alpha = 0.6742
## beta = 5e-04
## gamma = 0.2035
##
## Initial states:
## l = 83.2951
## b = 1.19
## s = 0.994 1.0082 1.0681 1.143 0.9848 0.9569
## 0.9179 0.9751 0.9961 0.9696 1.0118 0.9745
##
## sigma: 0.0448
##
## AIC AICc BIC
## 3907.523 3909.209 3974.550
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4692346 9.398919 6.907787 -0.4923675 3.416299 0.3403441
## ACF1
## Training set 0.06806468
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Feb 2014 392.1741 369.6635 414.6846 357.7472 426.6009
## Mar 2014 343.7213 319.9340 367.5087 307.3417 380.1010
## Apr 2014 382.7157 352.3878 413.0436 336.3331 429.0983
## May 2014 389.6814 355.3421 424.0207 337.1639 442.1988
## Jun 2014 382.8142 346.0026 419.6257 326.5158 439.1125
## Jul 2014 379.8759 340.5374 419.2144 319.7128 440.0389
## Aug 2014 397.7296 353.8049 441.6544 330.5526 464.9067
## Sep 2014 397.5302 351.0595 444.0010 326.4593 468.6011
## Oct 2014 392.0839 343.8564 440.3114 318.3263 465.8414
## Nov 2014 404.8268 352.6841 456.9694 325.0814 484.5721
## Dec 2014 399.2696 345.6330 452.9061 317.2396 481.2996
## Jan 2015 447.4498 384.9691 509.9306 351.8938 543.0059
## Feb 2015 405.3665 345.2714 465.4617 313.4589 497.2742
## Mar 2015 355.2521 300.8821 409.6221 272.1003 438.4039
## Apr 2015 395.5194 333.1517 457.8871 300.1362 490.9025
## May 2015 402.6825 337.3753 467.9897 302.8037 502.5612
## Jun 2015 395.5513 329.6751 461.4275 294.8024 496.3003
## Jul 2015 392.4809 325.4519 459.5099 289.9689 494.9929
## Aug 2015 410.8913 339.0225 482.7601 300.9774 520.8052
## Sep 2015 410.6497 337.1717 484.1278 298.2747 523.0247
## Oct 2015 404.9887 330.9347 479.0428 291.7328 518.2446
## Nov 2015 418.1152 340.0589 496.1715 298.7384 537.4920
## Dec 2015 412.3404 333.8179 490.8630 292.2505 532.4303
## Jan 2016 462.0588 372.3752 551.7424 324.8996 599.2180
autoplot(retail.ts)+autolayer(holt.fit1)+ylab('Retail Sales')+xlab("Year")
#damped=T
holt.fit2 <- hw(retail.ts,seasonal = "multiplicative",damped = TRUE)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Missing values encountered. Using longest contiguous portion of time series
summary(holt.fit2)
##
## Forecast method: Damped Holt-Winters' multiplicative method
##
## Model Information:
## Damped Holt-Winters' multiplicative method
##
## Call:
## hw(y = retail.ts, seasonal = "multiplicative", damped = TRUE)
##
## Smoothing parameters:
## alpha = 0.7222
## beta = 0.0065
## gamma = 0.1546
## phi = 0.98
##
## Initial states:
## l = 83.4463
## b = 0.4216
## s = 1.0229 0.968 1.064 1.0739 0.975 0.9727
## 0.9694 0.9971 0.9986 0.9687 1.0018 0.9879
##
## sigma: 0.0441
##
## AIC AICc BIC
## 3891.866 3893.756 3962.837
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.7419341 9.371794 6.802753 0.2471587 3.270906 0.3351692
## ACF1
## Training set 0.02072814
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Feb 2014 394.1059 371.8093 416.4025 360.0062 428.2056
## Mar 2014 346.3371 322.0905 370.5836 309.2552 383.4190
## Apr 2014 384.5351 353.2224 415.8477 336.6465 432.4237
## May 2014 389.7116 354.0165 425.4068 335.1206 444.3027
## Jun 2014 380.5691 342.1775 418.9606 321.8543 439.2839
## Jul 2014 374.9845 333.9202 416.0487 312.1821 437.7868
## Aug 2014 393.0053 346.7774 439.2333 322.3058 463.7049
## Sep 2014 390.0108 341.1307 438.8910 315.2551 464.7666
## Oct 2014 384.7794 333.7226 435.8362 306.6947 462.8640
## Nov 2014 398.3317 342.6599 454.0034 313.1890 483.4743
## Dec 2014 390.9617 333.6526 448.2709 303.3150 478.6085
## Jan 2015 436.3933 369.5415 503.2451 334.1523 538.6343
## Feb 2015 397.5615 333.0400 462.0830 298.8844 496.2386
## Mar 2015 349.3125 290.4728 408.1521 259.3250 439.2999
## Apr 2015 387.7719 320.1239 455.4199 284.3133 491.2306
## May 2015 392.9260 322.0683 463.7836 284.5586 501.2933
## Jun 2015 383.6448 312.2499 455.0397 274.4557 492.8339
## Jul 2015 377.9541 305.4806 450.4277 267.1154 488.7928
## Aug 2015 396.0552 317.9098 474.2006 276.5421 515.5683
## Sep 2015 392.9768 313.2908 472.6628 271.1075 514.8461
## Oct 2015 387.6470 306.9556 468.3383 264.2402 511.0538
## Nov 2015 401.2409 315.5922 486.8896 270.2525 532.2293
## Dec 2015 393.7601 307.6492 479.8710 262.0649 525.4554
## Jan 2016 439.4545 341.0831 537.8259 289.0085 589.9006
autoplot(retail.ts)+autolayer(holt.fit2)+ylab('Retail Sales')+xlab("Year")
##8c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
accuracy(holt.fit1)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4692346 9.398919 6.907787 -0.4923675 3.416299 0.3403441
## ACF1
## Training set 0.06806468
accuracy(holt.fit2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.7419341 9.371794 6.802753 0.2471587 3.270906 0.3351692
## ACF1
## Training set 0.02072814
In comparing the accuracy tables above, the model wit the trend damped is slightly better as it shows to have a lower RSME of 9.371
#8d. Check that the residuals from the best method look like white noise.
plot(holt.fit2$residuals,ylab = "Residuals",xlab = "Year", main = "Damped Model")
## This appears to look like white noise.
checkresiduals(holt.fit2)
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 23.257, df = 7, p-value = 0.001537
##
## Model df: 17. Total lags used: 24
#8e. 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?
##Therefore fitting the training set with 3 differetn methods: Seasonal naive, holt winter's multiplicativ and holt winters additve with box cox?
train<- window(retail.ts, end = c(2010,12))
test<- window(retail.ts, start = c(2011,1))
autoplot(retail.ts)+autolayer(train, series="Training")+
autolayer(test, series="Test")+ ylab("Turnover")
## Warning: Removed 1 row(s) containing missing values (geom_path).
fit_snaive <- snaive(train, h=36)
fit1_hw <- hw(train, h=36, seasonal='multiplicative', damped=F)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Missing values encountered. Using longest contiguous portion of time series
fit2_hw <- hw(train, h=36, seasonal='additive', damped=F, lambda='auto')
## Warning in ets(x, "AAA", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Missing values encountered. Using longest contiguous portion of time series
autoplot(test, series='Ground Truth') +
autolayer(fit_snaive, series='Seasonal Naive Forecast', PI=F) +
autolayer(fit1_hw, series="Holt-Winter's Forecast 1", PI=F) +
autolayer(fit2_hw, series="Holt-Winter's Forecast 2", PI=F) +
guides(colour=guide_legend(title="Legend")) +
ggtitle('Test Set Forecast') +
ylab('Turnover')
round(accuracy(fit_snaive),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 11.42 26.08 19.17 4.8 9.6 1 0.89
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?
10.
For this exercise use data set ukcars, the quarterly UK passenger vehicle production data from 1977Q1–2005Q1.
#10a. Plot the data and describe the main features of the series.
autoplot(ukcars)+ylab("Production")+xlab("Year")
There is an overall postive trend in the production of cars
#10b. Decompose the series using STL and obtain the seasonally adjusted data.
ukcars.stl <- stl(ukcars,s.window = "periodic")
summary(ukcars.stl)
## Call:
## stl(x = ukcars, s.window = "periodic")
##
## Time.series components:
## seasonal trend remainder
## Min. :-45.12476 Min. :215.7849 Min. :-46.73093
## 1st Qu.: -1.55632 1st Qu.:263.2765 1st Qu.: -9.08993
## Median : 21.01432 Median :329.0151 Median : -0.00269
## Mean : 0.22714 Mean :333.1863 Mean : 0.06435
## 3rd Qu.: 25.66677 3rd Qu.:403.8199 3rd Qu.: 8.88187
## Max. : 25.66677 Max. :455.5570 Max. : 35.75430
## IQR:
## STL.seasonal STL.trend STL.remainder data
## 27.22 140.54 17.97 124.83
## % 21.8 112.6 14.4 100.0
##
## Weights: all == 1
##
## Other components: List of 5
## $ win : Named num [1:3] 1131 7 5
## $ deg : Named int [1:3] 0 1 1
## $ jump : Named num [1:3] 114 1 1
## $ inner: int 2
## $ outer: int 0
#10c. 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.)
stlf(ukcars, etsmodel="AAN", damped=TRUE, h=8)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 415.0289 384.4576 445.6001 368.2742 461.7836
## 2005 Q3 364.2543 326.8700 401.6386 307.0799 421.4287
## 2005 Q4 400.8059 357.6702 443.9416 334.8355 466.7762
## 2006 Q1 432.5663 384.3596 480.7731 358.8404 506.2922
## 2006 Q2 415.0250 362.2312 467.8189 334.2838 495.7663
## 2006 Q3 364.2507 307.2369 421.2646 277.0556 451.4459
## 2006 Q4 400.8026 339.8596 461.7456 307.5983 494.0069
## 2007 Q1 432.5633 367.9289 497.1976 333.7136 531.4130
#10d.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).
stlf(ukcars, etsmodel="AAN", damped=FALSE, h=8)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 416.1915 385.7950 446.5881 369.7040 462.6791
## 2005 Q3 366.3343 328.7907 403.8780 308.9163 423.7524
## 2005 Q4 403.8032 360.2690 447.3375 337.2233 470.3831
## 2006 Q1 436.4809 387.6847 485.2771 361.8535 511.1083
## 2006 Q2 419.8568 366.3120 473.4016 337.9671 501.7465
## 2006 Q3 369.9996 312.0932 427.9060 281.4394 458.5598
## 2006 Q4 407.4685 345.5056 469.4313 312.7045 502.2325
## 2007 Q1 440.1461 374.3755 505.9167 339.5587 540.7336
#10e. Now use ets() to choose a seasonal model for the data
ets(ukcars)
## ETS(A,N,A)
##
## Call:
## ets(y = ukcars)
##
## Smoothing parameters:
## alpha = 0.6199
## gamma = 1e-04
##
## Initial states:
## l = 314.2568
## s = -1.7579 -44.9601 21.1956 25.5223
##
## sigma: 25.9302
##
## AIC AICc BIC
## 1277.752 1278.819 1296.844
## A,N,A was chosen
#10f. 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?
#ETS model
model1 <- ets(ukcars)
round(accuracy(model1),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.31 25.23 20.18 -0.16 6.63 0.66 0.03
# STL decomposition
ukcars.stl1 <- stl(ukcars, s.window = "periodic")
ukcarsSTL.ts <- ts(ukcars.stl1)
round(accuracy(ukcarsSTL.ts),2)
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## ME RMSE MAE MPE MAPE
## Training set NaN NaN NaN NaN NaN
#10g. Compare the forecasts from the three approaches? Which seems most reasonable?
forecast(model1,h=8)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 427.4885 394.2576 460.7195 376.6662 478.3109
## 2005 Q3 361.3329 322.2353 400.4305 301.5383 421.1275
## 2005 Q4 404.5358 360.3437 448.7280 336.9497 472.1219
## 2006 Q1 431.8154 383.0568 480.5741 357.2455 506.3854
## 2006 Q2 427.4885 374.5571 480.4200 346.5369 508.4401
## 2006 Q3 361.3329 304.5345 418.1313 274.4672 448.1986
## 2006 Q4 404.5358 344.1174 464.9542 312.1338 496.9378
## 2007 Q1 431.8154 367.9809 495.6500 334.1890 529.4419
#ETS model2
model2<-stlf(ukcars, etsmodel="AAN",damped=FALSE, h=8)
forecast(model2, h=8)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 416.1915 385.7950 446.5881 369.7040 462.6791
## 2005 Q3 366.3343 328.7907 403.8780 308.9163 423.7524
## 2005 Q4 403.8032 360.2690 447.3375 337.2233 470.3831
## 2006 Q1 436.4809 387.6847 485.2771 361.8535 511.1083
## 2006 Q2 419.8568 366.3120 473.4016 337.9671 501.7465
## 2006 Q3 369.9996 312.0932 427.9060 281.4394 458.5598
## 2006 Q4 407.4685 345.5056 469.4313 312.7045 502.2325
## 2007 Q1 440.1461 374.3755 505.9167 339.5587 540.7336
#STL decompostition
forecast(ukcars.stl1, h=8)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 Q2 427.2046 394.5838 459.8254 377.3154 477.0938
## 2005 Q3 361.0655 322.7354 399.3955 302.4447 419.6862
## 2005 Q4 404.6339 361.3411 447.9267 338.4232 470.8446
## 2006 Q1 431.8570 384.1145 479.5995 358.8412 504.8729
## 2006 Q2 427.2046 375.3931 479.0160 347.9658 506.4433
## 2006 Q3 361.0655 305.4822 416.6487 276.0582 446.0727
## 2006 Q4 404.6339 345.5190 463.7489 314.2254 495.0424
## 2007 Q1 431.8570 369.4098 494.3042 336.3523 527.3618
## THE RESULTS SEEM FAIRLY SIMILAR.
#10h. Check the residuals of your preferred model.
m.fcast <- forecast(model2,h=8)
m.fcast$residuals
## Qtr1 Qtr2 Qtr3 Qtr4
## 1977 -28.835239869 34.792950811 -17.843831364 31.143100613
## 1978 -18.175720851 1.789747863 -29.218141165 -65.825405579
## 1979 25.998433199 1.346804631 -74.484984962 29.002956755
## 1980 8.317191758 -40.520338332 -12.229254637 -28.393509990
## 1981 5.050525600 6.027685989 48.253212322 -9.425344894
## 1982 -22.321958202 -27.783679532 -0.012727342 16.652487401
## 1983 8.042045233 33.012673112 2.153599029 6.699556156
## 1984 -21.962669518 -38.468462137 6.355632624 -25.140545039
## 1985 53.635699541 10.703052762 -5.015190761 -10.836159374
## 1986 -21.242991619 7.760589156 13.877718528 17.882577013
## 1987 -7.407826721 14.778064951 23.680829001 -3.593278240
## 1988 -7.537614439 10.031461191 2.939348431 34.807075618
## 1989 9.840156810 2.642922463 -21.056725817 -19.263914257
## 1990 -3.993613576 -6.622635114 17.983390170 43.258441001
## 1991 -26.141859081 -11.319129411 -25.370095016 -20.330656797
## 1992 26.319778794 2.231717462 9.759212652 0.260887790
## 1993 14.823014726 13.807557287 4.426125704 -25.630712216
## 1994 3.574426241 24.304018286 17.680371091 16.389172378
## 1995 14.033379103 -3.362987339 -30.777150586 5.482016102
## 1996 21.166106881 11.769191549 24.587674772 31.871314218
## 1997 -35.179172026 0.422917653 -2.609546302 23.437070977
## 1998 -3.659459264 13.668629606 2.380732555 -41.199023082
## 1999 10.246771423 4.210664071 17.338356788 19.168463908
## 2000 3.052622471 -46.312369844 -55.885831158 -15.707943051
## 2001 -24.933899266 -6.688247260 14.202562977 24.358003624
## 2002 28.743000922 -24.122914045 34.619407687 -36.074702184
## 2003 -2.915573212 24.231270146 15.264612974 -16.733418686
## 2004 0.006014205 -0.769632650 0.489235653 -22.342510884
## 2005 -0.071449002
summary(m.fcast$residuals)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -74.4850 -17.8438 2.2317 -0.3413 14.8230 53.6357
11. For this exercise use data set visitors, the monthly Australian short-term overseas visitors data, May 1985–April 2005.
##11a. Make a time plot of your data and describe the main features of the series.
str(visitors)
## Time-Series [1:240] from 1985 to 2005: 75.7 75.4 83.1 82.9 77.3 ...
autoplot(visitors)+ylab("Short-term Oversea Visitors")+xlab("Year")
There is an overall positive trend with seasonal trends throughout the dataset.
#11b. 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.
#test set
visitor.test<- window(visitors, start = 2004)
#Forecasting 1 year out
visit.fcast <- hw(visitor.test,seasonal = "multiplicative", damped = FALSE, h=12)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Seasonal component could not be estimated
##Training set
#forecasting 1 year out
visit.fcast1 <- hw(visitors, seasonal = "multiplicative", damped = FALSE, h=12)
plot(visit.fcast1)
11c. Why is multiplicative seasonality necessary here?
# Multiplicative seasonality is necessary here because of the seasonal trend we observed in the original plot of the dataset.
#11d. Forecast the two-year test set using each of the following methods:
#i.an ETS model:
v.ets <- ets(visitor.test)
v.ets_fcast <-forecast(v.ets,h=24)
plot(v.ets_fcast)
checkresiduals(v.ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 1.3521, df = 3, p-value = 0.7168
##
## Model df: 2. Total lags used: 5
#ii.an additive ETS model applied to a Box-Cox transformed series:
v.bc <-ets(visitor.test, model = "ZZZ", lambda = "auto")
v.bc_fcast <- forecast(v.bc)
plot(v.bc_fcast)
checkresiduals(v.bc)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 1.3525, df = 3, p-value = 0.7167
##
## Model df: 2. Total lags used: 5
#iii. a seasonal naïve method:
v.snaive <- snaive(visitor.test, h=24)
v.snaive_fcast<- forecast(v.snaive, h=24)
plot(v.snaive_fcast)
checkresiduals(v.snaive)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 1.4962, df = 3, p-value = 0.6832
##
## Model df: 0. Total lags used: 3
#iv.an STL decomposition applied to the Box-Cox transformed data followed by an ETS model applied to the seasonally adjusted (transformed) data:
v.stl <- stl(visitors, s.window = "periodic")
v.stl_fcast <- forecast(v.stl,h=24)
plot(v.stl_fcast)
11e. Which method gives the best forecasts? Does it pass the residual tests?
## The seasonal naive method appears to give the best forecasts. It however,does not pass the residual test, which i think its due to how small the sample size happens to be.
##11f.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?
12. The fets() function below returns ETS forecasts.
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:
#a. Ausbeer
aus.test <- window(ausbeer,start = 2008)
aus.ets <-ets(aus.test)
round(accuracy(aus.ets),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.01 35.9 26.57 -0.68 6.12 3.39 -0.07
aus.snaive <- snaive(aus.test)
round(accuracy(aus.snaive))
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2 11 8 -1 2 1 0
aus.stl <- stlf(aus.test)
round(accuracy(aus.stl),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.86 5.88 4.36 -0.46 1.07 0.56 0.02
### The stlf() model gives the best forecast results for this data set
#b Bricksq
b.test <- window(bricksq,start = 1992)
b.ets <-ets(b.test)
round(accuracy(b.ets),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 12.93 33 28.85 2.5 6.45 1.27 -0.02
b.snaive <- snaive(b.test)
round(accuracy(b.snaive))
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 23 26 23 5 5 1 0
b.stl <- stlf(b.test)
round(accuracy(b.stl),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4.9 10.96 8.44 1.07 1.91 0.37 -0.3
##The stlf() model gives the best forecasting results for this dataset
#c. Dole
dole.test <- window(dole,start = 1990)
dole.ets <-ets(dole.test)
round(accuracy(dole.ets),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -583.69 12272.63 9763.9 -0.22 1.7 0.04 0.47
dole.snaive <- snaive(dole.test)
round(accuracy(dole.snaive))
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 221412 223754 221412 31 31 1 1
dole.stl <- stlf(dole.test)
round(accuracy(dole.stl),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -533.24 7483.45 6008.59 -0.05 1.03 0.03 0.1
## The stlf() and ets() models for this data set produce interestign outputs, however, i belive that the seasonal naive model gives the best forecasting results
#d. a10
a.test <- window(a10,start = 2006)
a.ets <-ets(a.test)
round(accuracy(a.ets),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0 3.74 2.86 -2.62 14.96 0.81 -0.08
a.snaive <- snaive(a.test)
round(accuracy(a.snaive))
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3 4 4 14 16 1 0
a.stl <- stlf(a.test)
round(accuracy(a.stl),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.03 1.15 0.94 -0.55 4.81 0.27 -0.28
## The ets()model gives the best forecasting results for this dataset
#e. h02
h.test <- window(h02,start = 2006)
h.ets <-ets(h.test)
round(accuracy(h.ets),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01 0.05 0.04 0.49 5.04 0.6 -0.33
h.snaive <- snaive(h.test)
round(accuracy(h.snaive))
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0 0 0 3 8 1 0
h.stl <- stlf(h.test)
round(accuracy(h.stl),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01 0.05 0.04 0.62 4.97 0.59 -0.32
## The ets() model gives the best forecasting results for this dataset
#f. usmelec
u.test <- window(usmelec,start = 2006)
u.ets <-ets(u.test)
round(accuracy(u.ets),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.17 8.49 6.62 -0.1 1.93 0.63 0.17
u.snaive <- snaive(u.test)
round(accuracy(u.snaive))
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0 14 11 0 3 1 1
u.stl <- stlf(u.test)
round(accuracy(u.stl),2)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.2 8.11 6.4 0.03 1.87 0.6 0.11
###The stlf() model gives the best forecasting results for this dataset
14.Use ets() on the following series:
#a. bicoal
ets(bicoal)
## ETS(M,N,N)
##
## Call:
## ets(y = bicoal)
##
## Smoothing parameters:
## alpha = 0.8205
##
## Initial states:
## l = 542.665
##
## sigma: 0.1262
##
## AIC AICc BIC
## 595.2499 595.7832 600.9253
#b. chicken
ets(chicken)
## ETS(M,N,N)
##
## Call:
## ets(y = chicken)
##
## Smoothing parameters:
## alpha = 0.98
##
## Initial states:
## l = 159.8322
##
## sigma: 0.1691
##
## AIC AICc BIC
## 635.2382 635.6018 641.9836
#c. dole
ets(dole)
## ETS(M,Ad,M)
##
## Call:
## ets(y = dole)
##
## Smoothing parameters:
## alpha = 0.697
## beta = 0.124
## gamma = 0.303
## phi = 0.902
##
## Initial states:
## l = 2708.6621
## b = 836.017
## s = 1.0404 0.8893 0.9103 1.0301 1.0576 1.0584
## 0.9801 0.9632 1.021 0.9838 1.0145 1.0514
##
## sigma: 0.0935
##
## AIC AICc BIC
## 10602.67 10604.30 10676.19
#d. usdeaths
ets(usdeaths)
## ETS(A,N,A)
##
## Call:
## ets(y = usdeaths)
##
## Smoothing parameters:
## alpha = 0.5972
## gamma = 0.0019
##
## Initial states:
## l = 9195.6403
## s = -62.6129 -270.0351 263.3823 -89.4907 1005.529 1662.647
## 795.2585 333.326 -551.161 -737.5102 -1552.872 -796.4611
##
## sigma: 294.4663
##
## AIC AICc BIC
## 1141.016 1149.587 1175.166
#e. lynx
ets(lynx)
## ETS(M,N,N)
##
## Call:
## ets(y = lynx)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 2372.8047
##
## sigma: 0.9594
##
## AIC AICc BIC
## 2058.138 2058.356 2066.346
#f. ibmclose
ets(ibmclose)
## ETS(A,N,N)
##
## Call:
## ets(y = ibmclose)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 459.9339
##
## sigma: 7.2637
##
## AIC AICc BIC
## 3648.450 3648.515 3660.182
#g. eggs
ets(eggs)
## ETS(M,N,N)
##
## Call:
## ets(y = eggs)
##
## Smoothing parameters:
## alpha = 0.8198
##
## Initial states:
## l = 278.8889
##
## sigma: 0.1355
##
## AIC AICc BIC
## 1043.286 1043.553 1050.916
#14b Find an example where it does not work well. Can you figure out why?
# It seems the ets() model did not work well on the usdeaths datasets. I'm not really sure why, but perhaps it can be explained by the patterns of the dataset?
15. Show that the point forecasts from an ETS(M,A,M) model are the same as those obtained using Holt-Winters’ multiplicative method.
#The three letter codes used in the ETS model matches the process used with the Holt-Winters multiplicative method.
16.Show that the forecast variance for an ETS(A,N,N) model is given by σ2[1+α2(h−1)].
## With the ETS(A,N,N) model, the traditional parameter region is 0<α<1 however, the admissible region is 0<α<2 . In this case, σ^2[1+α^2(h-1) shows that the same forecast variance can be found in the ETS(A,N,N) model.
17.Write down 95% prediction intervals for an ETS(A,N,N) model as a function of ℓT,α,h and σ, assuming normally distributed errors.
CHAPTER 8-ARIMA MODELS
1. Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers
#1a. Explain the differences among these figures. Do they all indicate that the data are white noise?
# The figures show a correlation between the lags of. The y-axis, which is the correlation has the same sale for all the plots, but the x-axis shows an increasing number of lags the longer the series gets. They do indicate white noise and this can be seen at each lag point between the dashed blue lines.
#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?
# We know that the critical value region is dependent on the sample size used. In this case, as the sample size increases, the critical values decrease. This explains why the critical value region gets smaller as the sample size increases.
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) + ylab("Close Price")+xlab("Day")
par(mfrow=c(1,2))
plot(acf(ibmclose))
plot(pacf(ibmclose))
There is a trend element throughout the plots. The ACF plots portray significant autocorrelations in the time series. This is observed in the pattern shown for the correlations when the lag increases- the correlation value does decrease as the lag increases.
The PACF plots show a strong outlier and the point after this outlier point are still in the critical value region?
3.For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
#3a. Usnetelec
plot(BoxCox(usnetelec,lambda = 0), main = "US Net Electricity Generation") #<- log transformation
#3b. Usgdp
plot(BoxCox(usgdp,lambda = 1/2),main = "US GDP" ) #<- Square root transformation
#3c. mcopper
plot(BoxCox(mcopper,lambda = 1/3), main = "Monthly Copper Prices") #<- cube root transformation
#3d. enplanements
plot(BoxCox(enplanements,lambda = 0), main = "Monthly US Domestic Enplanements") #<- Log transformation
#3e. visitors
plot(BoxCox(visitors,lambda = 1/3),main = "Monthly Australian Overseas Visitors") #<- cube root transformation
4. For the enplanements data, write down the differences you chose above using backshift operator notation.
acf(enplanements)
pacf(enplanements)
acf(enplanements - lag(enplanements,1))
pacf(enplanements - lag(enplanements,1))
There is a clear sign of seasonality that is present.This can be seen in the ACF plot. After the lag1 of the series is backshift by 1 period is taken, most of it is found to be covered by differencing of order 1.
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.
##Retail data has already been imported.
retailts <- ts(retail)
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
retail.transformation <-BoxCox(retailts,lambda = 0) #<- log transformation
6.Use R to simulate and plot some data from simple ARIMA models.
#6a. 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]
#6b. Produce a time plot for the series. How does the plot change as you change ϕ1?
plot(y)
## we find that as we decreasee ϕ1, the data becomes more condensed
y2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y2[i] <- 0.1*y2[i-1] + e[i]
plot(y2)
#6c. Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1.
input1 <- 0.6
yO <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
yO[i] <- input1*yO[i-1] + e[i]
plot(yO)
#6d. Produce a time plot for the series. How does the plot change as you change θ1?
input1 <- 0.6
input2 <- 0.1
input3 <- 1.5
## Test 1
y1 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y1[i] <- 0.6*y1[i-1] + e[i]
plot(y1)
## Test 2
y2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y2[i] <- 0.6*y2[i-1] + e[i]
plot(y2)
## Test 3
y3 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y3[i] <- 1.5*y3[i-1] + e[i]
plot(y3)
## The plot becomes less condensed as the theta vlaue increases.
#6e.Generate data from an ARMA(1,1) model with ϕ1=0.6, θ1=0.6 and σ^2=1.
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
arima(y)
##
## Call:
## arima(x = y)
##
## Coefficients:
## intercept
## -0.1458
## s.e. 0.1275
##
## sigma^2 estimated as 1.626: log likelihood = -166.21, aic = 336.43
yAr <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
yAr[i] <- 0.6*yAr[i-1] + 0.6*e[i-1] + e[i]
plot(yAr)
#6f. 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*e[i-1] + e[i]
plot(yAr2)
#6g.Graph the latter two series and compare them.
plot(yAr)
plot(yAr2)
## The second ARIMA model shows a much higher volatility in the data as compared to the first ARIMA model.
7.Consider wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.
#7a. By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data.
autoplot(wmurders) + ylab("Female Murder Rate")+ xlab("Years")
wmurders1 <- ets(wmurders, model="ZZZ", damped=NULL, alpha=NULL, beta=NULL,
gamma=NULL, phi=NULL, lambda=NULL, biasadj=FALSE,
additive.only=FALSE, restrict=TRUE,
allow.multiplicative.trend=FALSE)
wmurders1
## ETS(M,N,N)
##
## Call:
## ets(y = wmurders, model = "ZZZ", damped = NULL, alpha = NULL,
##
## Call:
## beta = NULL, gamma = NULL, phi = NULL, additive.only = FALSE,
##
## Call:
## lambda = NULL, biasadj = FALSE, restrict = TRUE, allow.multiplicative.trend = FALSE)
##
## Smoothing parameters:
## alpha = 0.9599
##
## Initial states:
## l = 2.4177
##
## sigma: 0.0617
##
## AIC AICc BIC
## 49.36113 49.83172 55.38313
## ETS autoselection chooses M,N,N
wmurders2 <- wmurders %>%
Arima(order=c(0,1,1), seasonal=c(0,1,1)) %>%
residuals() %>% ggtsdisplay()
wmurders2
wmurders3 <- arima(wmurders, order=c(3,1,1))
wmurders3
##
## Call:
## arima(x = wmurders, order = c(3, 1, 1))
##
## Coefficients:
## ar1 ar2 ar3 ma1
## 0.0683 0.3011 -0.0608 -0.1179
## s.e. 1.7595 0.1635 0.5328 1.7569
##
## sigma^2 estimated as 0.04104: log likelihood = 9.5, aic = -8.99
##Auto ARIMA suggested (1,2,1)
auto.arima(wmurders)
## 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.97
#7b. Should you include a constant in the model? Explain.
## I dont think so because including a constant may cause a drift/changes to the ARima model suggestions made
#7c. Write this model in terms of the backshift operator
wmAuto <- auto.arima(wmurders)
wmAuto
## 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.97
acf(wmurders)
pacf(wmurders)
wmurders.ts <- ts(wmurders)
## pacf(wmurders - lag(wmurders,1))
## pacf(wmurders.ts - lag(wmurders.ts,1))
## 7d. Fit the model using R and examine the residuals. Is the model satisfactory?
##
#7e. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.
forecast(wmurders,h=3)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 2.592549 2.387657 2.797442 2.279194 2.905905
## 2006 2.592549 2.308273 2.876826 2.157786 3.027313
## 2007 2.592549 2.246455 2.938644 2.063243 3.121855
#7f. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
plot(forecast(wmurders, h=3))
##7g. Does auto.arima() give the same model you have chosen? If not, which model do you think is better?
auto.arima(wmurders)
## 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.97
autoplot(forecast(auto.arima(wmurders),h=3))
## These two plots showcase similar forecasts, howeverm some differences can be seen in the forecasting range.
8.Consider austa, the total international visitors to Australia (in millions) for the period 1980-2015.
#8a. Use auto.arima() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise.
aus <- auto.arima(austa)
aus
## 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
checkresiduals(aus)
##
## 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: 7
## Auto Arima selected (0,1,1) as the proper model for this dataset
##Plot forecasts for the next 10 periods.
plot(forecast(aus, h=10))
##8b.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.
aus2 <- forecast(arima(austa, order = c(0,1,1)),h=10)
plot(aus2)
#8c. Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.
##aus3 <- forecast(arima(austa, order = c(2,1,3)),h=10)
##plot(aus3)
## Got an error message foe this commad above says non-stationary AR part from CSS
#8d. Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.
aus4 <- forecast(arima(austa, order=c(0,0,1)),h=10)
plot(aus4)
#removing MA term
plot(forecast(arima(austa),h=10))
#8e. Plot forecasts from an ARIMA(0,2,1) model with no constant.
aus5 <- forecast(arima(austa, order=c(0,2,1)),h=10)
plot(aus5)
9. For the usgdp series:
#9a. if necessary, find a suitable Box-Cox transformation for the data
plot(BoxCox(usgdp,lambda = 0), main = "US GDP")
usgdp1 <- BoxCox(usgdp,lambda = 0)
#9b.fit a suitable ARIMA model to the transformed data using auto.arima()
auto.arima(usgdp1)
## Series: usgdp1
## ARIMA(2,1,2) with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 drift
## 1.3594 -0.7732 -1.1002 0.6102 0.0084
## s.e. 0.1489 0.1774 0.2197 0.2285 0.0007
##
## sigma^2 estimated as 8.503e-05: log likelihood=773.76
## AIC=-1535.52 AICc=-1535.15 BIC=-1514.73
#9c.try some other plausible models by experimenting with the orders chosen
arima(usgdp1,order = c(0,0,1))
##
## Call:
## arima(x = usgdp1, order = c(0, 0, 1))
##
## Coefficients:
## ma1 intercept
## 1.0000 8.3967
## s.e. 0.0186 0.0373
##
## sigma^2 estimated as 0.08289: log likelihood = -43.93, aic = 93.87
arima(usgdp1,order = c(1,1,1))
##
## Call:
## arima(x = usgdp1, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.8512 -0.4140
## s.e. 0.0579 0.1117
##
## sigma^2 estimated as 9.845e-05: log likelihood = 753.49, aic = -1500.97
#9d. choose what you think is the best model and check the residual diagnostics
usgdp.best <- auto.arima(usgdp1)
checkresiduals(usgdp.best)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,2) with drift
## Q* = 6.3724, df = 3, p-value = 0.09483
##
## Model df: 5. Total lags used: 8
#9e. produce forecasts of your fitted model. Do the forecasts look reasonable?
forecast(usgdp.best,h=25)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2006 Q2 9.350972 9.339154 9.362789 9.332898 9.369045
## 2006 Q3 9.360099 9.341097 9.379101 9.331038 9.389161
## 2006 Q4 9.368788 9.343212 9.394364 9.329673 9.407904
## 2007 Q1 9.377005 9.345848 9.408162 9.329355 9.424655
## 2007 Q2 9.384918 9.349435 9.420401 9.330651 9.439185
## 2007 Q3 9.392785 9.354123 9.431446 9.333656 9.451913
## 2007 Q4 9.400821 9.359781 9.441862 9.338055 9.463587
## 2008 Q1 9.409127 9.366114 9.452139 9.343345 9.474909
## 2008 Q2 9.417665 9.372764 9.462565 9.348995 9.486334
## 2008 Q3 9.426311 9.379400 9.473223 9.354567 9.498056
## 2008 Q4 9.434927 9.385815 9.484039 9.359816 9.510037
## 2009 Q1 9.443414 9.391979 9.494849 9.364751 9.522077
## 2009 Q2 9.451753 9.398015 9.505491 9.369568 9.533938
## 2009 Q3 9.459988 9.404096 9.515881 9.374508 9.545469
## 2009 Q4 9.468198 9.410357 9.526039 9.379738 9.556658
## 2010 Q1 9.476452 9.416851 9.536054 9.385300 9.567605
## 2010 Q2 9.484788 9.423549 9.546027 9.391131 9.578445
## 2010 Q3 9.493199 9.430373 9.556025 9.397115 9.589283
## 2010 Q4 9.501650 9.437234 9.566066 9.403134 9.600166
## 2011 Q1 9.510097 9.444065 9.576129 9.409110 9.611084
## 2011 Q2 9.518507 9.450845 9.586169 9.415027 9.621987
## 2011 Q3 9.526871 9.457594 9.596148 9.420921 9.632821
## 2011 Q4 9.535200 9.464354 9.606047 9.426850 9.643551
## 2012 Q1 9.543518 9.471165 9.615871 9.432863 9.654173
## 2012 Q2 9.551847 9.478049 9.625644 9.438983 9.664710
autoplot(forecast(usgdp.best,h=25))
#9f. compare the results with what you would obtain using ets() (with no transformation).
ets(usgdp1)
## ETS(A,A,N)
##
## Call:
## ets(y = usgdp1)
##
## Smoothing parameters:
## alpha = 0.9959
## beta = 1e-04
##
## Initial states:
## l = 7.3524
## b = 0.0084
##
## sigma: 0.0099
##
## AIC AICc BIC
## -883.5955 -883.3357 -866.2552
usgdp.ets <- forecast(ets(usgdp1),h=25)
autoplot(usgdp.ets)
##The results look fairly similar between the two models
10. Consider austourists, the quarterly visitor nights (in millions) spent by international tourists to Australia for the period 1999–2015.
#10a. Describe the time plot.
autoplot(austourists)
##The plot above portrays a seasonal data with a positive trend
#10b.What can you learn from the ACF graph?
acf(austourists)
## Many of the values in this ACF plot lie outside of the critical value region indication that this may not have any white noise.
#10c.What can you learn from the PACF graph?
pacf(austourists)
##We see similar results to the acf in terms of values passing the critical region, signaling that in most cases there is not white noise present in the series.
#10d.Produce plots of the seasonally differenced data(1−B4)Yt. What model do these graphs suggest?
plot(diff(austourists, lag=4, differences=1), xlab="Year", ylab="Total Night Visitors", main="International Tourists to Australia", col="blue")
##We observe seasonality in the data but no direction in the trend of the data.
#10e.Does auto.arima() give the same model that you chose? If not, which model do you think is better?
auto.arima(austourists)
## 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.6
## Auto arima does select the same seaonal model.
#10f.Write the model in terms of the backshift operator, then without using the backshift operator.
### (1−B^4)Yt = BY(t-1) + e[i]
11.onsider 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.
## a. Examine the 12-month moving average of this series to see what kind of trend is involved.
mel1 <- ma(usmelec, order=12)
plot(mel1, col="red", main="Electricity Monthly Total Net Generation", xlab="Year", ylab="Electricity Generation (billions kWh")
## b. Do the data need transforming? If so, find a suitable transformation.
## The data does not seem to need transformation however, a log transformation will be used to examine the differences.
melLog <- BoxCox(usmelec, lambda=0)
plot(BoxCox(mel1, lambda=0), main="Electricity Monthly Total Net Generation", xlab="Year", ylab="Electricity Generation (billions kWh")
## c. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
## The data shows some small seasonal patterns throughout.
## Original set
##library(aTSA)
tsdisplay(diff(usmelec, 1))
## Log Transform
tsdisplay(diff(melLog, 1))
## d. 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?
## Autoselection chooses ARIMA(1,0,2)(0,1,1)[12] with drift
auto.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.22
## Alternate model
arima(usmelec, order=c(1,1,1))
##
## Call:
## arima(x = usmelec, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## -0.1626 0.3242
## s.e. 0.2454 0.2344
##
## sigma^2 estimated as 555.7: log likelihood = -2220.85, aic = 4447.69
## The autoselection model produces better results in terms of AIC.
## e. 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.
autoMel <- auto.arima(usmelec)
checkresiduals(autoMel)
##
## 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: 24
tsdisplay(autoMel$residuals)
melARIMAauto <- arima(usmelec)
autoMelF <- forecast(melARIMAauto)
autoMelF
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2013 259.5576 171.3822 347.733 124.705 394.4102
## Aug 2013 259.5576 171.3822 347.733 124.705 394.4102
## Sep 2013 259.5576 171.3822 347.733 124.705 394.4102
## Oct 2013 259.5576 171.3822 347.733 124.705 394.4102
## Nov 2013 259.5576 171.3822 347.733 124.705 394.4102
## Dec 2013 259.5576 171.3822 347.733 124.705 394.4102
## Jan 2014 259.5576 171.3822 347.733 124.705 394.4102
## Feb 2014 259.5576 171.3822 347.733 124.705 394.4102
## Mar 2014 259.5576 171.3822 347.733 124.705 394.4102
## Apr 2014 259.5576 171.3822 347.733 124.705 394.4102
## May 2014 259.5576 171.3822 347.733 124.705 394.4102
## Jun 2014 259.5576 171.3822 347.733 124.705 394.4102
## Jul 2014 259.5576 171.3822 347.733 124.705 394.4102
## Aug 2014 259.5576 171.3822 347.733 124.705 394.4102
## Sep 2014 259.5576 171.3822 347.733 124.705 394.4102
## Oct 2014 259.5576 171.3822 347.733 124.705 394.4102
## Nov 2014 259.5576 171.3822 347.733 124.705 394.4102
## Dec 2014 259.5576 171.3822 347.733 124.705 394.4102
## Jan 2015 259.5576 171.3822 347.733 124.705 394.4102
## Feb 2015 259.5576 171.3822 347.733 124.705 394.4102
## Mar 2015 259.5576 171.3822 347.733 124.705 394.4102
## Apr 2015 259.5576 171.3822 347.733 124.705 394.4102
## May 2015 259.5576 171.3822 347.733 124.705 394.4102
## Jun 2015 259.5576 171.3822 347.733 124.705 394.4102
## adf.test(autoMel$residuals)
## The autoselection model does not appear to be stationary based on the results
melA <- arima(usmelec, order=c(1,1,1))
tsdisplay(melA$residuals)
12.For the mcopper data:
#12a. if necessary, find a suitable Box-Cox transformation for the data;
autoplot(mcopper, col="dodgerblue") + ylab("Monthly Copper Prices") + xlab("Year")
mc <-BoxCox(mcopper, lambda=0)
autoplot(mc, col="dodgerblue") + ylab("Monthly Copper Prices") + xlab("Year")
#12b.fit a suitable ARIMA model to the transformed data using auto.arima();
auto.arima(mc)
## Series: mc
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## 0.3756
## s.e. 0.0385
##
## sigma^2 estimated as 0.003634: log likelihood=782.84
## AIC=-1561.69 AICc=-1561.66 BIC=-1553.02
#Arima (0,1,1) was selected
#12c.try some other plausible models by experimenting with the orders chosen;
#Model2
arima(mc,order = c(0,0,1))
##
## Call:
## arima(x = mc, order = c(0, 0, 1))
##
## Coefficients:
## ma1 intercept
## 0.9451 6.7305
## s.e. 0.0091 0.0267
##
## sigma^2 estimated as 0.1062: log likelihood = -169.03, aic = 344.06
#model3
arima(mc,order = c(1,1,1))
##
## Call:
## arima(x = mc, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## -0.0083 0.3825
## s.e. 0.1007 0.0912
##
## sigma^2 estimated as 0.003628: log likelihood = 782.85, aic = -1559.69
#12d. choose what you think is the best model and check the residual diagnostics;
##I think model 1 is the best model based on its AIC and BIC values.
#12e.produce forecasts of your fitted model. Do the forecasts look reasonable?
mc1 <- arima(mc)
mc.fcast <- forecast(mc)
mc.fcast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2007 8.133809 8.052148 8.215470 8.008919 8.258699
## Feb 2007 8.126612 8.004955 8.248269 7.940553 8.312670
## Mar 2007 8.120296 7.964292 8.276299 7.881709 8.358883
## Apr 2007 8.114754 7.927141 8.302367 7.827825 8.401683
## May 2007 8.109891 7.892420 8.327362 7.777298 8.442485
## Jun 2007 8.105624 7.859604 8.351644 7.729369 8.481879
## Jul 2007 8.101880 7.828388 8.375371 7.683611 8.520149
## Aug 2007 8.098594 7.798571 8.398618 7.639748 8.557441
## Sep 2007 8.095711 7.770005 8.421418 7.597586 8.593837
## Oct 2007 8.093182 7.742577 8.443787 7.556978 8.629386
## Nov 2007 8.090962 7.716193 8.465731 7.517802 8.664122
## Dec 2007 8.089014 7.690772 8.487256 7.479956 8.698072
## Jan 2008 8.087305 7.666247 8.508363 7.443352 8.731258
## Feb 2008 8.085805 7.642553 8.529057 7.407909 8.763701
## Mar 2008 8.084489 7.619635 8.549343 7.373556 8.795422
## Apr 2008 8.083334 7.597443 8.569226 7.340227 8.826441
## May 2008 8.082321 7.575929 8.588713 7.307861 8.856781
## Jun 2008 8.081432 7.555051 8.607813 7.276402 8.886462
## Jul 2008 8.080652 7.534770 8.626534 7.245797 8.915506
## Aug 2008 8.079967 7.515049 8.644885 7.216000 8.943935
## Sep 2008 8.079366 7.495855 8.662877 7.186963 8.971770
## Oct 2008 8.078839 7.477158 8.680521 7.158647 8.999032
## Nov 2008 8.078377 7.458928 8.697826 7.131011 9.025743
## Dec 2008 8.077971 7.441138 8.714803 7.104020 9.051922
autoplot(mc.fcast)+ ylab("Monthly Copper Prices") + xlab("Year")
#12f.ompare the results with what you would obtain using ets() (with no transformation).
mc2<- ets(mc)
mc.fcast1<- forecast(mc2)
mc.fcast1
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2007 8.133809 8.052148 8.215470 8.008919 8.258699
## Feb 2007 8.126612 8.004955 8.248269 7.940553 8.312670
## Mar 2007 8.120296 7.964292 8.276299 7.881709 8.358883
## Apr 2007 8.114754 7.927141 8.302367 7.827825 8.401683
## May 2007 8.109891 7.892420 8.327362 7.777298 8.442485
## Jun 2007 8.105624 7.859604 8.351644 7.729369 8.481879
## Jul 2007 8.101880 7.828388 8.375371 7.683611 8.520149
## Aug 2007 8.098594 7.798571 8.398618 7.639748 8.557441
## Sep 2007 8.095711 7.770005 8.421418 7.597586 8.593837
## Oct 2007 8.093182 7.742577 8.443787 7.556978 8.629386
## Nov 2007 8.090962 7.716193 8.465731 7.517802 8.664122
## Dec 2007 8.089014 7.690772 8.487256 7.479956 8.698072
## Jan 2008 8.087305 7.666247 8.508363 7.443352 8.731258
## Feb 2008 8.085805 7.642553 8.529057 7.407909 8.763701
## Mar 2008 8.084489 7.619635 8.549343 7.373556 8.795422
## Apr 2008 8.083334 7.597443 8.569226 7.340227 8.826441
## May 2008 8.082321 7.575929 8.588713 7.307861 8.856781
## Jun 2008 8.081432 7.555051 8.607813 7.276402 8.886462
## Jul 2008 8.080652 7.534770 8.626534 7.245797 8.915506
## Aug 2008 8.079967 7.515049 8.644885 7.216000 8.943935
## Sep 2008 8.079366 7.495855 8.662877 7.186963 8.971770
## Oct 2008 8.078839 7.477158 8.680521 7.158647 8.999032
## Nov 2008 8.078377 7.458928 8.697826 7.131011 9.025743
## Dec 2008 8.077971 7.441138 8.714803 7.104020 9.051922
autoplot(mc.fcast1)+ ylab("Monthly Copper Prices") + xlab("Year")
##The results appear to be similar however, slight differences can be observed. The ets model portrays a more negative trend compared to the ARIMA model
13. Choose one of the following seasonal time series: hsales, auscafe, qauselec, qcement, qgas.
#13a.Do the data need transforming? If so, find a suitable transformation.
str(hsales)
## Time-Series [1:275] from 1973 to 1996: 55 60 68 63 65 61 54 52 46 42 ...
autoplot(hsales)+ ylab("Monthly One-Family Home Sales") + xlab("Year")
hsales.bc <- BoxCox(hsales, lambda=0) ## Log Transform
autoplot(hsales.bc) + ylab("Monthly One-Family Home Sales") + xlab("Year")
#13b.Are the data stationary? If not, find an appropriate differencing which yields stationary data.
home <- BoxCox(hsales,lambda = 3)
autoplot(home) + ylab("Monthly One-Family Home Sales") + xlab("Year")
home2 <-BoxCox(hsales, lambda=.5)
autoplot(home2) + ylab("Monthly One-Family Home Sales") + xlab("Year")
home3 <- BoxCox(hsales, lambda=-2)
autoplot(home3) + ylab("Monthly One-Family Home Sales") + xlab("Year")
## None of the lambda values tested provided stationary transdformations
#13c.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 selection
auto.arima(hsales)
## Series: hsales
## ARIMA(1,0,0)(1,1,0)[12] with drift
##
## Coefficients:
## ar1 sar1 drift
## 0.8867 -0.4320 -0.0228
## s.e. 0.0294 0.0569 0.1642
##
## sigma^2 estimated as 27.92: log likelihood=-811.38
## AIC=1630.76 AICc=1630.92 BIC=1645.05
##(0,0,1)
arima(hsales, order=c(0,0,1))
##
## Call:
## arima(x = hsales, order = c(0, 0, 1))
##
## Coefficients:
## ma1 intercept
## 0.7465 52.2712
## s.e. 0.0309 0.8445
##
## sigma^2 estimated as 64.5: log likelihood = -963.54, aic = 1933.08
## (1,1,1)
arima(hsales, order=c(1,1,1))
##
## Call:
## arima(x = hsales, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.1713 -0.0142
## s.e. 0.1945 0.1902
##
## sigma^2 estimated as 40.03: log likelihood = -894.29, aic = 1794.58
## the autoselectiom model is the best according to its AIC value
#13d.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.
hsales.est <- auto.arima(hsales)
checkresiduals(hsales.est)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0)(1,1,0)[12] with drift
## Q* = 39.66, df = 21, p-value = 0.008177
##
## Model df: 3. Total lags used: 24
acf(hsales)
pacf(hsales)
##Many of the values seem to lie outside of the critical region as well as having noe white noise preseent
#13e.Forecast the next 24 months of data using your preferred model.
hsales.fcast <- forecast(hsales.est,h=24)
autoplot(hsales.fcast)
#13f.Compare the forecasts obtained using ets().
hsales.ets <- ets(hsales)
hets.fcast <- forecast(hsales.est,h=24)
## the ets model's forecasting region seems to be much wider than the ARIMA model.
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?
hsales.stlf <- stlf(hsales, method="arima")
hsales.stlf_fcast <- forecast(hsales.stlf,h=24)
plot(hsales.stlf_fcast)
##This methods seems to reveal a more focused forescast, which seems to be much more preferred.
15.or your retail time series (Exercise 5 above):
#15a develop an appropriate seasonal ARIMA model;
##auto.arima(retailts)
#15b. compare the forecasts with those you obtained in earlier chapters;
#15c. 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?
16. Consider sheep, the sheep population of England and Wales from 1867–1939
#16a. Produce a time plot of the time series.
autoplot(sheep) + ylab("Sheep Population") + xlab("Year")
#16b. 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)?
auto.arima(sheep)
## Series: sheep
## ARIMA(3,1,0)
##
## Coefficients:
## ar1 ar2 ar3
## 0.4210 -0.2018 -0.3044
## s.e. 0.1193 0.1363 0.1243
##
## sigma^2 estimated as 4991: log likelihood=-407.56
## AIC=823.12 AICc=823.71 BIC=832.22
## Arima (3,1,0 )was chosen
#16c. By examining the ACF and PACF of the differenced data, explain why this model is appropriate.
sheep1 <- auto.arima(sheep)
forecast(arima(sheep))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1940 1856.671 1573.106 2140.237 1422.995 2290.347
## 1941 1856.671 1573.106 2140.237 1422.995 2290.347
## 1942 1856.671 1573.106 2140.237 1422.995 2290.347
## 1943 1856.671 1573.106 2140.237 1422.995 2290.347
## 1944 1856.671 1573.106 2140.237 1422.995 2290.347
## 1945 1856.671 1573.106 2140.237 1422.995 2290.347
## 1946 1856.671 1573.106 2140.237 1422.995 2290.347
## 1947 1856.671 1573.106 2140.237 1422.995 2290.347
## 1948 1856.671 1573.106 2140.237 1422.995 2290.347
## 1949 1856.671 1573.106 2140.237 1422.995 2290.347
#ACF
acf(sheep)
#PACF
pacf(sheep)
##This model might be appropriate as it captures the seasonal trends that is eveident in the data.
## 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).
## 1940
## ϕ1=0.42, ϕ2=−0.20, ϕ3=−0.30
##σ^2[1+α^2(h-1)]
## 1941
## ϕ1=0.42, ϕ2=−0.20, ϕ3=−0.30
## 1942
## ϕ1=0.42, ϕ2=−0.20, ϕ3=−0.30
## e. Now fit the model in R and obtain the forecasts using forecast. How are they different from yours? Why?
sheep.ts <- ts(sheep)
forecast(arima(sheep.ts))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 74 1856.671 1573.106 2140.237 1422.995 2290.347
## 75 1856.671 1573.106 2140.237 1422.995 2290.347
## 76 1856.671 1573.106 2140.237 1422.995 2290.347
## 77 1856.671 1573.106 2140.237 1422.995 2290.347
## 78 1856.671 1573.106 2140.237 1422.995 2290.347
## 79 1856.671 1573.106 2140.237 1422.995 2290.347
## 80 1856.671 1573.106 2140.237 1422.995 2290.347
## 81 1856.671 1573.106 2140.237 1422.995 2290.347
## 82 1856.671 1573.106 2140.237 1422.995 2290.347
## 83 1856.671 1573.106 2140.237 1422.995 2290.347
#The results are somewhat similar
17. The annual bituminous coal production in the United States from 1920 to 1968 is in data set bicoal.
#17a.Produce a time plot of the data.
autoplot(bicoal) + ylab("Bituminous Coal Production") + xlab("Year")
#17b. 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)?
auto.arima(bicoal)
## Series: bicoal
## ARIMA(4,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ar4 mean
## 0.8334 -0.3443 0.5525 -0.3780 481.5221
## s.e. 0.1366 0.1752 0.1733 0.1414 21.0591
##
## sigma^2 estimated as 2795: log likelihood=-262.05
## AIC=536.1 AICc=538.1 BIC=547.45
##The autoselection model chose ARIMA (4,0,0)
#17c. Explain why this model was chosen using the ACF and PACF.
coal1 <- auto.arima(bicoal)
coal.fcast <- forecast(arima(bicoal))
#ACF
acf(bicoal)
#PACF
pacf(bicoal)
##As stated before, this model would be appropriate as it captures the seasonality in the dataset.
#17d. 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).
## 1969
## ϕ1=0.83, ϕ2=−0.34, ϕ3=0.55, ϕ4=−0.38
##σ^2[1+α^2(h-1)]
## 1970
## ϕ1=0.83, ϕ2=−0.34, ϕ3=0.55, ϕ4=−0.38
## 1971
## ϕ1=0.83, ϕ2=−0.34, ϕ3=0.55, ϕ4=−0.38
#17e. Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?
coal.ts <- ts(bicoal)
forecast(arima(coal.ts))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 50 478.1837 378.8997 577.4677 326.3419 630.0254
## 51 478.1837 378.8997 577.4677 326.3419 630.0254
## 52 478.1837 378.8997 577.4677 326.3419 630.0254
## 53 478.1837 378.8997 577.4677 326.3419 630.0254
## 54 478.1837 378.8997 577.4677 326.3419 630.0254
## 55 478.1837 378.8997 577.4677 326.3419 630.0254
## 56 478.1837 378.8997 577.4677 326.3419 630.0254
## 57 478.1837 378.8997 577.4677 326.3419 630.0254
## 58 478.1837 378.8997 577.4677 326.3419 630.0254
## 59 478.1837 378.8997 577.4677 326.3419 630.0254
#The results are similar from the ARIMA model but may be subject to some constraints within R.
18. Before doing this exercise, you will need to install the Quandl package in R using
library(Quandl)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
#18a. Select a time series from Quandl. Then copy its short URL and import the data using y <- Quandl("?????", api_key="?????", type="ts") (Replace each ????? with the appropriate values.)