Week 5 - Chapter 7 and 8 Exercises
Bryce Anderson
11/12/2019
## R Packages Used
library(readr)
library(ggplot2)
library(forecast)
library(GGally)
library(fpp2)
library(seasonal)
library(apricom)
library(psych)
library(dplyr)Chapter 7 - Exponential Smoothing
- Consider the pigs series — the number of pigs slaughtered in Victoria each month.
## a. 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 ...pigsSES <- ses(pigs)
pigsSES## 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
## Jan 1996 98816.41 83448.52 114184.3 75313.25 122319.6
## Feb 1996 98816.41 82955.06 114677.8 74558.56 123074.2
## Mar 1996 98816.41 82476.49 115156.3 73826.66 123806.2
## Apr 1996 98816.41 82011.54 115621.3 73115.58 124517.2
## May 1996 98816.41 81559.12 116073.7 72423.66 125209.2
## Jun 1996 98816.41 81118.26 116514.6 71749.42 125883.4pigsSES2 <- ses(pigs, h=4) ## 4 quarters per year
## b. 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.
pigsSES$upper[1,"95%"]## 95%
## 119020.8pigsSES$lower[1,"95%"]## 95%
## 78611.97res <- sd(pigsSES$residuals)
res## [1] 10273.69pigsSES$mean[1] + 1.96*res## [1] 118952.8pigsSES$mean[1] - 1.96*res## [1] 78679.97- 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()?
pigsf1 <- ses(pigs, h=4, alpha=0.1, initial="simple")
pigsf1## 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.2forecast(pigsf1, 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## The forecasts appear to give slightly different forecasts when incorporating the forecast() functionModify 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?
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.
Data set books contains the daily sales of paperback and hardcover books at the same store. The task is to forecast the next four days’ sales for paperback and hardcover books.
## a. Plot the series and discuss the main features of the data.
autoplot(books, ylab="Book sales", xlab="Number of Days", main="Book Sales over Time")
## The paperback and hardcover book sales show seasonal patterns and typically act as the inverse of each other. This means that typically, when paperback sales increase, hardcover sales decrease and vice versa.
## b. Use the ses() function to forecast each series, and plot the forecasts.
books <- data.frame(books)
hardTS <- ts(books$Hardcover)
paperTS <- ts(books$Paperback)
## Predicting 4 days of Hardcover sales
hF <- ses(hardTS, h=4)
hF## 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## Predicting 4 days of Paperback sales
pF <- ses(paperTS, h=4)
pF## 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(hF)
autoplot(pF)
## c. Compute the RMSE values for the training data in each case.
round(forecast::accuracy(pF),2)[2]## [1] 33.64round(forecast::accuracy(hF),2)[2]## [1] 31.93## Hardcover RMSE = 33.64
## Paperback RMSE = 31.93- Data set books contains the daily sales of paperback and hardcover books at the same store. The task is to forecast the next four days’ sales for paperback and hardcover books.
## a. Now apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
## 4 day predictions for each book type
h1 <- holt(hardTS, h=4)
h1## 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.3334h2 <- holt(paperTS, h=4)
h1## 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## b. Compare the RMSE measures of Holt’s method for the two series to those of simple exponential smoothing in the previous question. (Remember that Holt’s method is using one more parameter than SES.) Discuss the merits of the two forecasting methods for these data sets.
## Hardcover
round(forecast::accuracy(h1),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.14 27.19 23.16 -2.11 12.16 0.69 -0.03## Paperback
round(forecast::accuracy(h2),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -3.72 31.14 26.18 -5.51 15.58 0.66 -0.18## c. Compare the forecasts for the two series using both methods. Which do you think is best?
## It appears that the Holt's method performed best overall and had the lowest RMSE value on Hardcover.
## d. 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
h1$mean[1] + round(forecast::accuracy(h1),2)*1.96## ME RMSE MAE MPE MAPE MASE
## Training set 249.8995 303.4663 295.5675 246.0383 274.0075 251.5263
## ACF1
## Training set 250.1151## Lower Hardcover
h1$mean[1] - round(forecast::accuracy(h1),2)*1.96## ME RMSE MAE MPE MAPE MASE
## Training set 250.4483 196.8815 204.7803 254.3095 226.3403 248.8215
## ACF1
## Training set 250.2327## Upper Paperback
h2$mean[1] + round(forecast::accuracy(h2),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
h2$mean[1] - round(forecast::accuracy(h2),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- 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?
## Holt Method
holt(eggs, h=100)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1994 59.78553313 24.970286 94.60078 6.540207 113.0309
## 1995 57.06372643 12.206005 101.92145 -11.540238 125.6677
## 1996 54.34191973 1.308668 107.37517 -26.765440 135.4493
## 1997 51.62011302 -8.488401 111.72863 -40.307926 143.5482
## 1998 48.89830632 -17.537664 115.33428 -52.706742 150.5034
## 1999 46.17649962 -26.035964 118.38896 -64.262933 156.6159
## 2000 43.45469292 -34.106500 121.01589 -75.164916 162.0743
## 2001 40.73288622 -41.832449 123.29822 -85.539898 167.0057
## 2002 38.01107951 -49.273098 125.29526 -95.478551 171.5007
## 2003 35.28927281 -56.472472 127.05102 -105.048206 175.6268
## 2004 32.56746611 -63.464327 128.59926 -114.300487 179.4354
## 2005 29.84565941 -70.275216 129.96654 -123.276007 182.9673
## 2006 27.12385271 -76.926479 131.17418 -132.007398 186.2551
## 2007 24.40204600 -83.435566 132.23966 -140.521350 189.3254
## 2008 21.68023930 -89.816966 133.17744 -148.840022 192.2005
## 2009 18.95843260 -96.082866 133.99973 -156.982051 194.8989
## 2010 16.23662590 -102.243631 134.71688 -164.963291 197.4365
## 2011 13.51481920 -108.308165 135.33780 -172.797358 199.8270
## 2012 10.79301249 -114.284184 135.87021 -180.496053 202.0821
## 2013 8.07120579 -120.178426 136.32084 -188.069681 204.2121
## 2014 5.34939909 -125.996817 136.69562 -195.527304 206.2261
## 2015 2.62759239 -131.744601 136.99979 -202.876944 208.1321
## 2016 -0.09421431 -137.426446 137.23802 -210.125737 209.9373
## 2017 -2.81602102 -143.046526 137.41448 -217.280071 211.6480
## 2018 -5.53782772 -148.608596 137.53294 -224.345686 213.2700
## 2019 -8.25963442 -154.116046 137.59678 -231.327766 214.8085
## 2020 -10.98144112 -159.571946 137.60906 -238.231008 216.2681
## 2021 -13.70324782 -164.979092 137.57260 -245.059687 217.6532
## 2022 -16.42505453 -170.340036 137.48993 -251.817706 218.9676
## 2023 -19.14686123 -175.657116 137.36339 -258.508640 220.2149
## 2024 -21.86866793 -180.932478 137.19514 -265.135773 221.3984
## 2025 -24.59047463 -186.168101 136.98715 -271.702130 222.5212
## 2026 -27.31228133 -191.365811 136.74125 -278.210504 223.5859
## 2027 -30.03408804 -196.527300 136.45912 -284.663482 224.5953
## 2028 -32.75589474 -201.654137 136.14235 -291.063466 225.5517
## 2029 -35.47770144 -206.747783 135.79238 -297.412688 226.4573
## 2030 -38.19950814 -211.809598 135.41058 -303.713228 227.3142
## 2031 -40.92131484 -216.840851 134.99822 -309.967029 228.1244
## 2032 -43.64312155 -221.842733 134.55649 -316.175908 228.8897
## 2033 -46.36492825 -226.816355 134.08650 -322.341569 229.6117
## 2034 -49.08673495 -231.762762 133.58929 -328.465610 230.2921
## 2035 -51.80854165 -236.682939 133.06586 -334.549533 230.9324
## 2036 -54.53034835 -241.577809 132.51711 -340.594753 231.5341
## 2037 -57.25215506 -246.448244 131.94393 -346.602603 232.0983
## 2038 -59.97396176 -251.295068 131.34714 -352.574344 232.6264
## 2039 -62.69576846 -256.119059 130.72752 -358.511164 233.1196
## 2040 -65.41757516 -260.920954 130.08580 -364.414190 233.5790
## 2041 -68.13938186 -265.701450 129.42269 -370.284491 234.0057
## 2042 -70.86118857 -270.461210 128.73883 -376.123079 234.4007
## 2043 -73.58299527 -275.200863 128.03487 -381.930915 234.7649
## 2044 -76.30480197 -279.921006 127.31140 -387.708914 235.0993
## 2045 -79.02660867 -284.622209 126.56899 -393.457946 235.4047
## 2046 -81.74841537 -289.305013 125.80818 -399.178839 235.6820
## 2047 -84.47022208 -293.969936 125.02949 -404.872385 235.9319
## 2048 -87.19202878 -298.617469 124.23341 -410.539337 236.1553
## 2049 -89.91383548 -303.248085 123.42041 -416.180415 236.3527
## 2050 -92.63564218 -307.862233 122.59095 -421.796308 236.5250
## 2051 -95.35744888 -312.460345 121.74545 -427.387675 236.6728
## 2052 -98.07925559 -317.042831 120.88432 -432.955146 236.7966
## 2053 -100.80106229 -321.610088 120.00796 -438.499326 236.8972
## 2054 -103.52286899 -326.162494 119.11676 -444.020793 236.9751
## 2055 -106.24467569 -330.700413 118.21106 -449.520103 237.0308
## 2056 -108.96648239 -335.224193 117.29123 -454.997790 237.0648
## 2057 -111.68828910 -339.734170 116.35759 -460.454367 237.0778
## 2058 -114.41009580 -344.230666 115.41047 -465.890327 237.0701
## 2059 -117.13190250 -348.713991 114.45019 -471.306143 237.0423
## 2060 -119.85370920 -353.184443 113.47702 -476.702272 236.9949
## 2061 -122.57551590 -357.642310 112.49128 -482.079153 236.9281
## 2062 -125.29732261 -362.087868 111.49322 -487.437211 236.8426
## 2063 -128.01912931 -366.521384 110.48313 -492.776852 236.7386
## 2064 -130.74093601 -370.943117 109.46124 -498.098471 236.6166
## 2065 -133.46274271 -375.353314 108.42783 -503.402447 236.4770
## 2066 -136.18454941 -379.752215 107.38312 -508.689149 236.3200
## 2067 -138.90635612 -384.140052 106.32734 -513.958929 236.1462
## 2068 -141.62816282 -388.517050 105.26072 -519.212132 235.9558
## 2069 -144.34996952 -392.883424 104.18349 -524.449088 235.7491
## 2070 -147.07177622 -397.239385 103.09583 -529.670117 235.5266
## 2071 -149.79358292 -401.585134 101.99797 -534.875530 235.2884
## 2072 -152.51538963 -405.920870 100.89009 -540.065628 235.0348
## 2073 -155.23719633 -410.246781 99.77239 -545.240700 234.7663
## 2074 -157.95900303 -414.563051 98.64505 -550.401029 234.4830
## 2075 -160.68080973 -418.869861 97.50824 -555.546889 234.1853
## 2076 -163.40261643 -423.167382 96.36215 -560.678543 233.8733
## 2077 -166.12442314 -427.455784 95.20694 -565.796249 233.5474
## 2078 -168.84622984 -431.735228 94.04277 -570.900257 233.2078
## 2079 -171.56803654 -436.005874 92.86980 -575.990809 232.8547
## 2080 -174.28984324 -440.267874 91.68819 -581.068139 232.4885
## 2081 -177.01164994 -444.521379 90.49808 -586.132476 232.1092
## 2082 -179.73345665 -448.766534 89.29962 -591.184043 231.7171
## 2083 -182.45526335 -453.003480 88.09295 -596.223054 231.3125
## 2084 -185.17707005 -457.232353 86.87821 -601.249720 230.8956
## 2085 -187.89887675 -461.453288 85.65553 -606.264245 230.4665
## 2086 -190.62068345 -465.666414 84.42505 -611.266828 230.0255
## 2087 -193.34249016 -469.871857 83.18688 -616.257661 229.5727
## 2088 -196.06429686 -474.069741 81.94115 -621.236934 229.1083
## 2089 -198.78610356 -478.260186 80.68798 -626.204828 228.6326
## 2090 -201.50791026 -482.443308 79.42749 -631.161524 228.1457
## 2091 -204.22971696 -486.619220 78.15979 -636.107193 227.6478
## 2092 -206.95152367 -490.788035 76.88499 -641.042007 227.1390
## 2093 -209.67333037 -494.949859 75.60320 -645.966131 226.6195lambda <- BoxCox.lambda(eggs)
holt1 <- holt(eggs, h=100)
holt2 <- holt(eggs, h=100, lambda=lambda)
holt3 <- holt(eggs, h=100, lambda=lambda, damped=TRUE)
autoplot(holt(eggs,h=100)) + autolayer(holt1, series="holts") + autolayer(holt2, series="boxcox") + autolayer(holt3, series="boxcox_damped") + ylab("Paperback") + xlab("Year")
- Recall your retail time series data (from Exercise 3 in Section 2.10).
## Importing the Australian Retail data set
library(readr)
retailR <- read_csv("retailCSV.csv")
library(readxl)
retailXL <- read_excel("retail.xlsx")
## a. Why is multiplicative seasonality necessary for this series?
retailTS <- ts(retailXL[,20], frequency=12, start=c(1982,4))
autoplot(retailTS)
## Multiplicative seasonality is necessary for this series becuase of the increasing trend seen in the autoplot graph,
## b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
## Damped = FALSE
holtFit1 <- hw(retailTS, seasonal="multiplicative", damped=FALSE, h=30)
autoplot(retailTS) + autolayer(holtFit1) + ylab("Retail Sales") + xlab("Year")
## Damped = TRUE
holtFit2 <- hw(retailTS, seasonal="multiplicative", damped=TRUE, h=30)
autoplot(retailTS) + autolayer(holtFit2) + ylab("Retail Sales") + xlab("Year")
## c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
holtFit1$mean[1]## [1] 392.1741holtFit2$mean[1]## [1] 394.1059round(forecast::accuracy(holtFit1),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.47 9.4 6.91 -0.49 3.42 0.34 0.07round(forecast::accuracy(holtFit2),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.74 9.37 6.8 0.25 3.27 0.34 0.02## In comparing the above mean values, the damped fit shows slightly better results. The damped model also shows a slightly lower RMSE value.
## d. Check that the residuals from the best method look like white noise.
plot(holtFit2$residuals, ylab="Residuals", xlab="Year", main="Damped Model Resdiuals")
checkresiduals(holtFit2)
##
## 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## e. 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?
testRetail <- window(retailTS, start=2011)
## Seasonal Naïve (12 months * 3 years = 36 periods (Jan 2011-Dec 2013))
retailSN <- snaive(testRetail, h=36)
round(forecast::accuracy(retailSN),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 5.25 27.31 21.93 1.29 5.91 1 0.63## RSME value of 27.31- 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?
lambdaTS <- BoxCox.lambda(retailTS)
bc1 <- BoxCox(retailTS, lambdaTS)
bc1## Jan Feb Mar Apr May Jun Jul
## 1982 NA 4.171794 4.165599 4.122013
## 1983 4.285306 4.246614 4.168702 4.154121 4.111036 4.126364 4.111036
## 1984 4.213951 4.288901 4.220778 4.261562 4.207068 4.244726 4.191123
## 1985 4.349467 4.308405 4.188100 4.252250 4.258779 4.277161 4.239040
## 1986 4.480012 4.402869 4.290692 4.355293 4.387039 4.415323 4.361903
## 1987 4.575353 4.458994 4.394990 4.402869 4.410674 4.410674 4.396572
## 1988 4.574073 4.553985 4.477144 4.535421 4.499822 4.515066 4.494209
## 1989 4.578543 4.565707 4.461193 4.504702 4.474266 4.494209 4.562467
## 1990 4.728110 4.767675 4.655198 4.704320 4.703223 4.704320 4.691054
## 1991 4.742381 4.681528 4.577906 4.638757 4.643487 4.628016 4.617743
## 1992 4.760013 4.754865 4.682092 4.757443 4.726509 4.729708 4.663875
## 1993 4.698266 4.632209 4.528014 4.481443 4.546079 4.492800 4.538098
## 1994 4.737651 4.805454 4.646430 4.774259 4.647018 4.655198 4.654616
## 1995 4.743429 4.776275 4.669039 4.756928 4.756928 4.751762 4.717380
## 1996 5.006971 5.036951 4.952566 4.998147 4.975867 4.957441 4.921857
## 1997 5.000074 5.011158 4.877582 4.970701 4.939840 4.954601 4.903046
## 1998 5.012296 4.972692 4.871767 4.901748 4.782293 4.755897 4.766657
## 1999 4.819909 4.787274 4.660415 4.743952 4.707604 4.722762 4.703772
## 2000 4.741857 4.697713 4.618958 4.692168 4.723835 4.703772 4.736068
## 2001 4.986887 4.968306 4.883358 4.984930 4.844387 4.808849 4.835526
## 2002 4.987668 5.027698 4.829881 4.941080 4.933194 4.901748 4.888214
## 2003 5.029186 5.014946 4.828463 4.918469 5.020972 5.032153 4.989619
## 2004 5.174021 5.153496 5.031041 5.075427 5.157005 5.147720 5.118223
## 2005 5.198904 5.136374 5.043550 5.097615 5.109839 5.097615 5.104768
## 2006 5.319718 5.234730 5.109839 5.190692 5.177444 5.170270 5.156687
## 2007 5.243377 5.195872 5.085208 5.183941 5.174021 5.174332 5.185788
## 2008 5.309463 5.237911 5.209428 5.256462 5.248807 5.264605 5.232118
## 2009 5.385715 5.353505 5.210622 5.305751 5.316576 5.331136 5.285571
## 2010 5.552045 5.524641 5.376501 5.427775 5.487176 5.476025 5.484616
## 2011 5.622297 5.526064 5.381606 5.432571 5.410180 5.369153 5.364960
## 2012 5.502996 5.426629 5.292083 5.366936 5.429148 5.402185 5.409946
## 2013 5.550273 5.415778 5.306813 5.435072 5.430519 5.416242 5.411116
## 2014 5.562793
## Aug Sep Oct Nov Dec
## 1982 4.140355 4.122013 4.137147 4.191123 4.192129
## 1983 4.164562 4.169734 4.160401 4.115444 4.145676
## 1984 4.205091 4.220778 4.210024 4.280791 4.290692
## 1985 4.280791 4.311904 4.311904 4.366008 4.416095
## 1986 4.363547 4.400513 4.379817 4.426072 4.457525
## 1987 4.483584 4.471378 4.489267 4.498422 4.504702
## 1988 4.508863 4.503310 4.481443 4.510246 4.513004
## 1989 4.563116 4.584889 4.581087 4.624404 4.641717
## 1990 4.679271 4.696053 4.570865 4.646430 4.647018
## 1991 4.634596 4.659258 4.624404 4.674171 4.678706
## 1992 4.636382 4.634000 4.647605 4.677008 4.667894
## 1993 4.576630 4.605496 4.642898 4.694390 4.698818
## 1994 4.689939 4.693835 4.701574 4.654616 4.650533
## 1995 4.803507 4.817038 4.815599 4.921012 4.934444
## 1996 4.917619 4.937354 4.900448 4.967906 4.944378
## 1997 4.939840 4.948073 4.957036 4.975867 4.952566
## 1998 4.793214 4.722226 4.715759 4.781793 4.751762
## 1999 4.673033 4.716300 4.738704 4.750725 4.744475
## 2000 4.875797 4.873560 4.927755 4.993506 4.894790
## 2001 4.860012 4.855447 4.857733 4.929432 4.869069
## 2002 4.949710 4.922703 4.912932 4.986104 4.927335
## 2003 5.046102 5.058746 5.041722 5.116219 5.070847
## 2004 5.180235 5.131791 5.152856 5.178686 5.115549
## 2005 5.204934 5.142877 5.158913 5.255332 5.236756
## 2006 5.243950 5.225992 5.211815 5.178065 5.179306
## 2007 5.266558 5.266558 5.235310 5.268229 5.253635
## 2008 5.273222 5.263486 5.212708 5.295588 5.249091
## 2009 5.378206 5.383059 5.397920 5.458125 5.430291
## 2010 5.561627 5.522196 5.542150 5.549287 5.504874
## 2011 5.351497 5.344935 5.374793 5.373080 5.336529
## 2012 5.465595 5.436433 5.442079 5.433254 5.444999
## 2013 5.413683 5.455030 5.418098 5.465595 5.486324
## 2014lambdaTS <- BoxCox.lambda(testRetail)
bc1 <- BoxCox(testRetail, lambdaTS)
bc1## Jan Feb Mar Apr May Jun Jul
## 2011 5061.528 4285.012 3340.600 3646.792 3508.894 3269.857 3246.380
## 2012 4117.653 3609.668 2864.726 3257.422 3625.360 3460.951 3507.481
## 2013 4468.125 3542.860 2937.995 3662.534 3633.928 3545.695 3514.548
## 2014 4565.935
## Aug Sep Oct Nov Dec
## 2011 3172.157 3136.611 3301.703 3292.000 3091.674
## 2012 3860.340 3671.129 3707.010 3651.083 3725.711
## 2013 3530.111 3790.669 3557.042 3860.340 4000.857
## 2014retailSTL <- stl(bc1[,1], s.window="periodic")
## h=36 (3 years * 12 months)
forecast1 <- forecast(retailSTL, h=36)
round(forecast::accuracy(forecast1),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -14 159.53 137.79 -0.45 3.83 0.44 -0.01## This method does not perform as well as the previous test performed on the created test set.- For this exercise use data set ukcars, the quarterly UK passenger vehicle production data from 1977Q1–2005Q1.
## a. Plot the data and describe the main features of the series.
autoplot(ukcars) + ylab("Vehicle Production") + xlab("Year")
## We see an overall positive trend in vehicle production with a seasonal components.
## b. Decompose the series using STL and obtain the seasonally adjusted data.
ukcarsSTL <- stl(ukcars, s.window="periodic")
ukcarsSTL## Call:
## stl(x = ukcars, s.window = "periodic")
##
## Components
## seasonal trend remainder
## 1977 Q1 25.66677 318.5295 -13.825249626
## 1977 Q2 21.01432 324.6948 25.341911916
## 1977 Q3 -45.12476 330.9708 -15.176049732
## 1977 Q4 -1.55632 334.7004 10.735966254
## 1978 Q1 25.66677 334.0139 -1.189641048
## 1978 Q2 21.01432 320.7959 21.011758342
## 1978 Q3 -45.12476 300.9207 5.485038711
## 1978 Q4 -1.55632 288.6423 -46.730931160
## 1979 Q1 25.66677 276.1832 23.532046833
## 1979 Q2 21.01432 267.2862 28.399513279
## 1979 Q3 -45.12476 261.5016 -45.223827211
## 1979 Q4 -1.55632 252.2466 6.526756134
## 1980 Q1 25.66677 247.9784 24.481824797
## 1980 Q2 21.01432 239.4564 -9.006725034
## 1980 Q3 -45.12476 223.2180 3.461793086
## 1980 Q4 -1.55632 215.7849 -21.630623541
## 1981 Q1 25.66677 220.2357 -0.250430089
## 1981 Q2 21.01432 233.6016 -9.089932014
## 1981 Q3 -45.12476 241.8626 28.523188447
## 1981 Q4 -1.55632 240.7936 -1.026285436
## 1982 Q1 25.66677 230.3050 1.413217582
## 1982 Q2 21.01432 222.3072 -14.860538991
## 1982 Q3 -45.12476 222.2015 -1.705761122
## 1982 Q4 -1.55632 231.0511 -3.032794821
## 1983 Q1 25.66677 244.9445 -4.461265566
## 1983 Q2 21.01432 257.0965 9.140190429
## 1983 Q3 -45.12476 263.2765 7.731214593
## 1983 Q4 -1.55632 257.4187 9.450666648
## 1984 Q1 25.66677 245.2947 1.797509364
## 1984 Q2 21.01432 233.9046 -20.784913580
## 1984 Q3 -45.12476 228.4033 13.183489842
## 1984 Q4 -1.55632 237.1828 -30.075524618
## 1985 Q1 25.66677 247.7529 17.863310026
## 1985 Q2 21.01432 258.0678 5.339865349
## 1985 Q3 -45.12476 258.1793 8.516496453
## 1985 Q4 -1.55632 250.2516 2.001681939
## 1986 Q1 25.66677 246.4844 -18.394218794
## 1986 Q2 21.01432 250.3747 -4.372980363
## 1986 Q3 -45.12476 259.5206 5.992168272
## 1986 Q4 -1.55632 267.0536 12.303689570
## 1987 Q1 25.66677 274.7015 -17.135282062
## 1987 Q2 21.01432 282.7800 -1.722305893
## 1987 Q3 -45.12476 290.1326 14.712201563
## 1987 Q4 -1.55632 294.2213 4.992997392
## 1988 Q1 25.66677 295.1384 -14.676141546
## 1988 Q2 21.01432 300.0450 1.046637055
## 1988 Q3 -45.12476 313.0266 -11.178866207
## 1988 Q4 -1.55632 324.8864 18.546947724
## 1989 Q1 25.66677 331.9778 -1.640598597
## 1989 Q2 21.01432 329.0151 11.510606843
## 1989 Q3 -45.12476 321.2409 -5.683114815
## 1989 Q4 -1.55632 312.5455 0.115820802
## 1990 Q1 25.66677 308.5850 -7.563781787
## 1990 Q2 21.01432 315.1534 -9.108686479
## 1990 Q3 -45.12476 326.9998 -7.618043743
## 1990 Q4 -1.55632 333.4080 35.754296494
## 1991 Q1 25.66677 332.0422 -11.545969032
## 1991 Q2 21.01432 318.1429 9.053818020
## 1991 Q3 -45.12476 307.8729 -12.740133598
## 1991 Q4 -1.55632 306.5401 -12.465738791
## 1992 Q1 25.66677 310.8608 6.790443375
## 1992 Q2 21.01432 318.7590 3.655635693
## 1992 Q3 -45.12476 325.3152 -4.804421022
## 1992 Q4 -1.55632 332.1973 -0.893989077
## 1993 Q1 25.66677 340.5537 -1.699468477
## 1993 Q2 21.01432 344.6912 12.742515376
## 1993 Q3 -45.12476 344.0280 1.894724342
## 1993 Q4 -1.55632 343.3536 -10.040311206
## 1994 Q1 25.66677 348.0348 -11.165536675
## 1994 Q2 21.01432 359.2368 8.881872441
## 1994 Q3 -45.12476 374.3526 -5.905812714
## 1994 Q4 -1.55632 385.3718 8.016551486
## 1995 Q1 25.66677 388.6002 7.378976268
## 1995 Q2 21.01432 384.3830 11.425715058
## 1995 Q3 -45.12476 381.6207 -24.782895640
## 1995 Q4 -1.55632 383.7478 -0.289469197
## 1996 Q1 25.66677 394.5880 2.727254896
## 1996 Q2 21.01432 411.2334 -4.525738213
## 1996 Q3 -45.12476 424.2350 -2.260256000
## 1996 Q4 -1.55632 428.3269 31.809443631
## 1997 Q1 25.66677 427.3670 -16.808791145
## 1997 Q2 21.01432 423.5993 -3.126575508
## 1997 Q3 -45.12476 427.6065 -12.915722425
## 1997 Q4 -1.55632 434.9783 17.301020647
## 1998 Q1 25.66677 442.4171 -5.641897664
## 1998 Q2 21.01432 442.2694 4.948265854
## 1998 Q3 -45.12476 437.2456 11.515205896
## 1998 Q4 -1.55632 433.4005 -17.896137200
## 1999 Q1 25.66677 431.3057 3.523496436
## 1999 Q2 21.01432 439.1628 -11.245141015
## 1999 Q3 -45.12476 450.8969 2.014845258
## 1999 Q4 -1.55632 455.5570 15.407293417
## 2000 Q1 25.66677 444.3540 24.290258982
## 2000 Q2 21.01432 421.1754 -8.949762594
## 2000 Q3 -45.12476 395.4301 -15.199289613
## 2000 Q4 -1.55632 375.7606 4.590727142
## 2001 Q1 25.66677 367.5362 -6.102990303
## 2001 Q2 21.01432 369.1569 -17.776203042
## 2001 Q3 -45.12476 380.2567 0.658094924
## 2001 Q4 -1.55632 394.1493 4.486990545
## 2002 Q1 25.66677 403.5459 20.542371093
## 2002 Q2 21.01432 408.8754 -27.637700761
## 2002 Q3 -45.12476 404.8623 32.109499840
## 2002 Q4 -1.55632 404.8806 -17.434230990
## 2003 Q1 25.66677 406.8758 -8.217613506
## 2003 Q2 21.01432 412.2684 -0.002686131
## 2003 Q3 -45.12476 418.0638 18.273925233
## 2003 Q4 -1.55632 419.2933 -8.996989753
## 2004 Q1 25.66677 416.2195 3.571765783
## 2004 Q2 21.01432 413.9253 -6.737593135
## 2004 Q3 -45.12476 410.1435 14.029258400
## 2004 Q4 -1.55632 407.2154 -11.617080187
## 2005 Q1 25.66677 403.8199 3.309325330summary(ukcarsSTL)## 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## c. Forecast the next two years of the series using an additive damped trend method applied to the seasonally adjusted data. (This can be done in one step using stlf() with arguments etsmodel="AAN", damped=TRUE.)
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## d. Forecast the next two years of the series using Holt’s linear method applied to the seasonally adjusted data (as before but with damped=FALSE).
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## e. 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## Chose ETS(A,N,A)
## f. 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
m1 <- ets(ukcars)
round(forecast::accuracy(m1),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 Model
ukcarsSTL <- stl(ukcars, s.window="periodic")
ukcarsSTLTS <- ts(ukcarsSTL)
round(forecast::accuracy(ukcarsSTLTS),2)## ME RMSE MAE MPE MAPE
## Training set NaN NaN NaN NaN NaN## g. Compare the forecasts from the three approaches? Which seems most reasonable?
## ETS Model Forecast
forecast(m1, 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 Model 2 Forecast
m2 <- stlf(ukcars, etsmodel="AAN", damped=FALSE, h=8)
forecast(m2, 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## SSTL Decomposition Model Forecast
forecast(ukcarsSTL, 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 are all fairly similar. I decided to choose model 1 because the ets() function decided the best model formulation for this dataset.
## h. Check the residuals of your preferred model.
mForecast1 <- forecast(m2, h=8)
mForecast1$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.071449002summary(mForecast1$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -74.4850 -17.8438 2.2317 -0.3413 14.8230 53.6357- For this exercise use data set visitors, the monthly Australian short-term overseas visitors data, May 1985–April 2005.
## a. Make a time plot of your data and describe the main features of the series.
autoplot(visitors) + ylab("Short-term Overseas Visitors") + xlab("Year")
## We see an overall positive trend with some seasonal trends at different points within each year.
## b. 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.
## Testing set
testVisitor <- window(visitors, start=2004)
## Forecasting one year (12 months)
hw(testVisitor, seasonal="multiplicative", damped=FALSE, h=12)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 441.3111 354.24112 528.3810 308.149069 574.4731
## Jun 2005 439.6645 325.75496 553.5740 265.454900 613.8740
## Jul 2005 438.0179 294.86439 581.1714 219.083511 656.9522
## Aug 2005 436.3713 261.77064 610.9719 169.342651 703.3999
## Sep 2005 434.7247 226.62773 642.8216 116.467867 752.9815
## Oct 2005 433.0781 189.55935 676.5968 60.648333 805.5078
## Nov 2005 431.4315 150.66814 712.1948 2.041018 860.8219
## Dec 2005 429.7849 110.04129 749.5284 -59.220727 918.7904
## Jan 2006 428.1383 67.75421 788.5223 -123.021579 979.2981
## Feb 2006 426.4917 23.87303 829.1103 -189.260392 1042.2437
## Mar 2006 424.8450 -21.54358 871.2337 -257.847447 1107.5375
## Apr 2006 423.1984 -68.44309 914.8400 -328.702399 1175.0993## Training Set
## Forecasting one year (12 months)
hw(visitors, seasonal="multiplicative", damped=FALSE, h=12)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 363.6434 337.2901 389.9967 323.3395 403.9474
## Jun 2005 389.5974 356.8764 422.3184 339.5550 439.6399
## Jul 2005 482.7709 437.0247 528.5172 412.8081 552.7338
## Aug 2005 428.5001 383.4951 473.5050 359.6709 497.3293
## Sep 2005 420.4548 372.1236 468.7860 346.5386 494.3710
## Oct 2005 475.5963 416.3290 534.8636 384.9547 566.2379
## Nov 2005 503.4598 435.9471 570.9725 400.2080 606.7116
## Dec 2005 608.3779 521.1122 695.6436 474.9165 741.8393
## Jan 2006 455.1525 385.6597 524.6452 348.8725 561.4324
## Feb 2006 509.0590 426.6702 591.4478 383.0562 635.0618
## Mar 2006 496.8101 411.8771 581.7431 366.9162 626.7039
## Apr 2006 441.9586 362.3923 521.5249 320.2724 563.6448## c. Why is multiplicative seasonality necessary here?
## Multiplicative seasonality is necessary here becuase of the increasing trend seen in data. The multiplicative method captures trends in the data that the additive wouldn't.
## d. Forecast the two-year test set using each of the following methods:
## i. an ETS model
vis1 <- ets(testVisitor)
forecast(vis1, h=24)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 465.1445 385.2825 545.0065 343.0061 587.2829
## Jun 2005 465.1445 380.5627 549.7263 335.7878 594.5012
## Jul 2005 465.1445 376.0927 554.1963 328.9515 601.3375
## Aug 2005 465.1445 371.8366 558.4524 322.4424 607.8466
## Sep 2005 465.1445 367.7664 562.5226 316.2175 614.0715
## Oct 2005 465.1445 363.8596 566.4294 310.2426 620.0464
## Nov 2005 465.1445 360.0980 570.1910 304.4897 625.7993
## Dec 2005 465.1445 356.4665 573.8225 298.9358 631.3532
## Jan 2006 465.1445 352.9525 577.3365 293.5616 636.7274
## Feb 2006 465.1445 349.5452 580.7438 288.3507 641.9383
## Mar 2006 465.1445 346.2356 584.0534 283.2890 647.0000
## Apr 2006 465.1445 343.0156 587.2734 278.3645 651.9245
## May 2006 465.1445 339.8784 590.4106 273.5665 656.7225
## Jun 2006 465.1445 336.8178 593.4712 268.8857 661.4033
## Jul 2006 465.1445 333.8285 596.4605 264.3141 665.9749
## Aug 2006 465.1445 330.9058 599.3832 259.8442 670.4448
## Sep 2006 465.1445 328.0454 602.2436 255.4696 674.8195
## Oct 2006 465.1445 325.2435 605.0455 251.1844 679.1047
## Nov 2006 465.1445 322.4966 607.7924 246.9833 683.3057
## Dec 2006 465.1445 319.8016 610.4875 242.8616 687.4274
## Jan 2007 465.1445 317.1556 613.1334 238.8150 691.4740
## Feb 2007 465.1445 314.5562 615.7328 234.8395 695.4495
## Mar 2007 465.1445 312.0008 618.2882 230.9315 699.3576
## Apr 2007 465.1445 309.4874 620.8016 227.0875 703.2015 ## ii. an additive ETS model applied to a Box-Cox transformed series
vis2 <- ets(testVisitor, model="AZZ", lambda="auto")
forecast(vis2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 465.1511 385.2891 545.0131 343.0127 587.2895
## Jun 2005 465.1511 380.5737 549.7286 335.8011 594.5011
## Jul 2005 465.1511 376.1076 554.1946 328.9708 601.3314
## Aug 2005 465.1511 371.8551 558.4471 322.4672 607.8350
## Sep 2005 465.1511 367.7882 562.5141 316.2473 614.0549
## Oct 2005 465.1511 363.8844 566.4178 310.2771 620.0252
## Nov 2005 465.1511 360.1257 570.1765 304.5286 625.7737
## Dec 2005 465.1511 356.4969 573.8053 298.9788 631.3234
## Jan 2006 465.1511 352.9854 577.3168 293.6085 636.6937
## Feb 2006 465.1511 349.5806 580.7216 288.4013 641.9009
## Mar 2006 465.1511 346.2733 584.0290 283.3431 646.9591
## Apr 2006 465.1511 343.0555 587.2467 278.4220 651.8803
## May 2006 465.1511 339.9204 590.3819 273.6272 656.6750
## Jun 2006 465.1511 336.8618 593.4404 268.9496 661.3526
## Jul 2006 465.1511 333.8745 596.4277 264.3809 665.9213
## Aug 2006 465.1511 330.9537 599.3485 259.9139 670.3883
## Sep 2006 465.1511 328.0951 602.2071 255.5421 674.7601
## Oct 2006 465.1511 325.2950 605.0072 251.2596 679.0426
## Nov 2006 465.1511 322.5498 607.7524 247.0612 683.2410
## Dec 2006 465.1511 319.8564 610.4458 242.9421 687.3601
## Jan 2007 465.1511 317.2121 613.0901 238.8980 691.4042
## Feb 2007 465.1511 314.6143 615.6879 234.9249 695.3773
## Mar 2007 465.1511 312.0605 618.2417 231.0192 699.2830
## Apr 2007 465.1511 309.5486 620.7536 227.1776 703.1246 ## iii. a seasonal naïve method
vis3 <- snaive(testVisitor, h=24)
forecast(vis3, h=24)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 341.3 279.6464 402.9536 247.0090 435.5910
## Jun 2005 367.3 305.6464 428.9536 273.0090 461.5910
## Jul 2005 472.0 410.3464 533.6536 377.7090 566.2910
## Aug 2005 405.8 344.1464 467.4536 311.5090 500.0910
## Sep 2005 395.6 333.9464 457.2536 301.3090 489.8910
## Oct 2005 449.9 388.2464 511.5536 355.6090 544.1910
## Nov 2005 479.9 418.2464 541.5536 385.6090 574.1910
## Dec 2005 593.1 531.4464 654.7536 498.8090 687.3910
## Jan 2006 462.4 400.7464 524.0536 368.1090 556.6910
## Feb 2006 501.6 439.9464 563.2536 407.3090 595.8910
## Mar 2006 504.7 443.0464 566.3536 410.4090 598.9910
## Apr 2006 409.5 347.8464 471.1536 315.2090 503.7910
## May 2006 341.3 254.1087 428.4913 207.9524 474.6476
## Jun 2006 367.3 280.1087 454.4913 233.9524 500.6476
## Jul 2006 472.0 384.8087 559.1913 338.6524 605.3476
## Aug 2006 405.8 318.6087 492.9913 272.4524 539.1476
## Sep 2006 395.6 308.4087 482.7913 262.2524 528.9476
## Oct 2006 449.9 362.7087 537.0913 316.5524 583.2476
## Nov 2006 479.9 392.7087 567.0913 346.5524 613.2476
## Dec 2006 593.1 505.9087 680.2913 459.7524 726.4476
## Jan 2007 462.4 375.2087 549.5913 329.0524 595.7476
## Feb 2007 501.6 414.4087 588.7913 368.2524 634.9476
## Mar 2007 504.7 417.5087 591.8913 371.3524 638.0476
## Apr 2007 409.5 322.3087 496.6913 276.1524 542.8476summary(vis3)##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = testVisitor, h = 24)
##
## Residual sd: 27.039
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 42.025 48.10855 42.025 8.609435 8.609435 1 -0.4227005
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 341.3 279.6464 402.9536 247.0090 435.5910
## Jun 2005 367.3 305.6464 428.9536 273.0090 461.5910
## Jul 2005 472.0 410.3464 533.6536 377.7090 566.2910
## Aug 2005 405.8 344.1464 467.4536 311.5090 500.0910
## Sep 2005 395.6 333.9464 457.2536 301.3090 489.8910
## Oct 2005 449.9 388.2464 511.5536 355.6090 544.1910
## Nov 2005 479.9 418.2464 541.5536 385.6090 574.1910
## Dec 2005 593.1 531.4464 654.7536 498.8090 687.3910
## Jan 2006 462.4 400.7464 524.0536 368.1090 556.6910
## Feb 2006 501.6 439.9464 563.2536 407.3090 595.8910
## Mar 2006 504.7 443.0464 566.3536 410.4090 598.9910
## Apr 2006 409.5 347.8464 471.1536 315.2090 503.7910
## May 2006 341.3 254.1087 428.4913 207.9524 474.6476
## Jun 2006 367.3 280.1087 454.4913 233.9524 500.6476
## Jul 2006 472.0 384.8087 559.1913 338.6524 605.3476
## Aug 2006 405.8 318.6087 492.9913 272.4524 539.1476
## Sep 2006 395.6 308.4087 482.7913 262.2524 528.9476
## Oct 2006 449.9 362.7087 537.0913 316.5524 583.2476
## Nov 2006 479.9 392.7087 567.0913 346.5524 613.2476
## Dec 2006 593.1 505.9087 680.2913 459.7524 726.4476
## Jan 2007 462.4 375.2087 549.5913 329.0524 595.7476
## Feb 2007 501.6 414.4087 588.7913 368.2524 634.9476
## Mar 2007 504.7 417.5087 591.8913 371.3524 638.0476
## Apr 2007 409.5 322.3087 496.6913 276.1524 542.8476 ## iv. an STL decomposition applied to the Box-Cox transformed data followed by an ETS model applied to the seasonally adjusted (transformed) data.
vis4 <- stl(visitors, s.window="periodic")
forecast(vis4, h=24)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2005 401.2570 375.1438 427.3701 361.3204 441.1935
## Jun 2005 416.4592 387.3886 445.5298 371.9996 460.9189
## Jul 2005 471.8865 440.1316 503.6413 423.3216 520.4514
## Aug 2005 443.2371 409.0069 477.4674 390.8865 495.5878
## Sep 2005 435.2928 398.7535 471.8320 379.4108 491.1747
## Oct 2005 472.7590 434.0474 511.4706 413.5547 531.9633
## Nov 2005 490.5152 449.7460 531.2844 428.1641 552.8663
## Dec 2005 560.3010 517.5724 603.0296 494.9533 625.6488
## Jan 2006 464.0268 419.4240 508.6296 395.8128 532.2408
## Feb 2006 496.6757 450.2737 543.0777 425.7099 567.6414
## Mar 2006 488.9345 440.7998 537.0692 415.3188 562.5502
## Apr 2006 454.2768 404.4689 504.0847 378.1022 530.4514
## May 2006 418.2968 366.8695 469.7241 339.6456 496.9481
## Jun 2006 433.4991 380.5012 486.4969 352.4459 514.5523
## Jul 2006 488.9263 434.4025 543.4502 405.5394 572.3133
## Aug 2006 460.2770 404.2682 516.2858 374.6189 545.9351
## Sep 2006 452.3326 394.8766 509.7887 364.4612 540.2040
## Oct 2006 489.7989 430.9306 548.6671 399.7677 579.8301
## Nov 2006 507.5551 447.3072 567.8030 415.4139 599.6963
## Dec 2006 577.3409 515.7436 638.9381 483.1360 671.5457
## Jan 2007 481.0667 418.1485 543.9848 384.8416 577.2917
## Feb 2007 513.7155 449.5031 577.9280 415.5111 611.9200
## Mar 2007 505.9744 440.4927 571.4560 405.8288 606.1199
## Apr 2007 471.3167 404.5894 538.0439 369.2662 573.3672## e. Which method gives the best forecasts? Does it pass the residual tests?
## The Seasoanl naive method appears to give the best results in this case. The sample size is small with just two additional years forecasted. It does not pass the residual tests, but would likely be improved with a larger set.
checkresiduals(vis3)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 1.4962, df = 3, p-value = 0.6832
##
## Model df: 0. Total lags used: 3summary(vis3$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 6.70 33.33 44.50 42.02 53.20 72.40 12## f. 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?
## Producing errors when trying different time series with the tsCV() function
## tsCV(tsCV(visitors, forecastfunction, h = 1, window = NULL, xreg = NULL,initial = 0))- The fets() function below returns ETS forecasts.
fets <- function(y, h) {
forecast(ets(y), h = h)
}
fets(y=qcement, h=4)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2014 Q2 2.521670 2.368825 2.674515 2.287913 2.755427
## 2014 Q3 2.610968 2.412978 2.808958 2.308169 2.913767
## 2014 Q4 2.571906 2.344302 2.799510 2.223816 2.919996
## 2015 Q1 2.268522 2.042601 2.494443 1.923006 2.614039## a. Apply tsCV() for a forecast horizon of h=4, for both ETS and seasonal naïve methods to the qcement data, (Hint: use the newly created fets() and the existing snaive() functions as your forecast function arguments.)
f1 <- fets(y=qcement, h=4)
f1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2014 Q2 2.521670 2.368825 2.674515 2.287913 2.755427
## 2014 Q3 2.610968 2.412978 2.808958 2.308169 2.913767
## 2014 Q4 2.571906 2.344302 2.799510 2.223816 2.919996
## 2015 Q1 2.268522 2.042601 2.494443 1.923006 2.614039f2 <- snaive(qcement, h=4)
f2## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2014 Q2 2.528 2.356456 2.699544 2.265646 2.790354
## 2014 Q3 2.637 2.465456 2.808544 2.374646 2.899354
## 2014 Q4 2.565 2.393456 2.736544 2.302646 2.827354
## 2015 Q1 2.229 2.057456 2.400544 1.966646 2.491354## ETS Model
qcementTS <- ts(qcement)
tsCV1 <- tsCV(qcementTS, fets)
## Seasonal Naive Model
fets2 <- function(y, h) {
forecast(snaive(y), h = h)
}
tsCV2 <- tsCV(qcementTS, fets2)
## b. Compute the MSE of the resulting 4-step-ahead errors. (Hint: make sure you remove missing values.) Why are there missing values? Comment on which forecasts are more accurate. Is this what you expected?
f1na <- na.omit(f1)
f2na <- na.omit(f2)
round(forecast::accuracy(f1na),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0 0.08 0.06 -0.09 3.61 0.55 -0.04round(forecast::accuracy(f2na),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.03 0.13 0.1 2.34 6.77 1 0.61summary(f1na)##
## Forecast method: ETS(M,A,M)
##
## Model Information:
## ETS(M,A,M)
##
## Call:
## ets(y = y)
##
## Smoothing parameters:
## alpha = 0.7505
## beta = 0.003
## gamma = 1e-04
##
## Initial states:
## l = 0.5043
## b = 0.0081
## s = 1.0297 1.0488 1.0164 0.9052
##
## sigma: 0.0473
##
## AIC AICc BIC
## 6.783846 7.591021 37.843192
##
## Error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.0006237224 0.07829922 0.05603263 -0.08782872 3.606505
## MASE ACF1
## Training set 0.5485176 -0.03632985
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2014 Q2 2.521670 2.368825 2.674515 2.287913 2.755427
## 2014 Q3 2.610968 2.412978 2.808958 2.308169 2.913767
## 2014 Q4 2.571906 2.344302 2.799510 2.223816 2.919996
## 2015 Q1 2.268522 2.042601 2.494443 1.923006 2.614039summary(f2na)##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = qcement, h = 4)
##
## Residual sd: 0.1297
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03419651 0.1338564 0.1021528 2.338509 6.768085 1
## ACF1
## Training set 0.613902
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2014 Q2 2.528 2.356456 2.699544 2.265646 2.790354
## 2014 Q3 2.637 2.465456 2.808544 2.374646 2.899354
## 2014 Q4 2.565 2.393456 2.736544 2.302646 2.827354
## 2015 Q1 2.229 2.057456 2.400544 1.966646 2.491354## The missing values are the result of values not populated in the forecast or data. Missing values can cause issues in forecasting and comparing descriptiive statistics of models.
## The f1 model (ETS model) seems to perform better based on the error measures produced above.- Compare ets(), snaive() and stlf() on the following six time series. For stlf(), you might need to use a Box-Cox transformation. Use a test set of three years to decide what gives the best forecasts. ausbeer, bricksq, dole, a10, h02, usmelec.
## a. ausbeer
testAus <- window(ausbeer, start=2008)
ets1 <- ets(testAus)
round(forecast::accuracy(ets1),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.01 35.9 26.57 -0.68 6.12 3.39 -0.07sn1 <- snaive(testAus)
round(forecast::accuracy(sn1),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.17 11.16 7.83 -0.62 2 1 0.12st1 <- stlf(testAus)
round(forecast::accuracy(st1),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() function provides the best forecasting model of the three for ausbeer.
## b. bricksq
testBri <- window(bricksq, start=1992)
ets2 <- ets(testBri)
round(forecast::accuracy(ets2),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 12.93 33 28.85 2.5 6.45 1.27 -0.02sn2 <- snaive(testBri)
round(forecast::accuracy(sn2),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 22.71 26.07 22.71 5 5 1 -0.3st2 <- stlf(testBri)
round(forecast::accuracy(st2),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() function provides the best forecasting model of the three for brickq.
## c. dole
testDol <- window(dole, start=1990)
ets3 <- ets(testDol)
round(forecast::accuracy(ets3),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -583.69 12272.63 9763.9 -0.22 1.7 0.04 0.47sn3 <- snaive(testDol)
round(forecast::accuracy(sn3),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 221411.8 223754 221411.8 31.02 31.02 1 0.74st3 <- stlf(testDol)
round(forecast::accuracy(st3),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 dole test set produces interesting results in both the stlf() and ets() models. If the negative values for ME are interpreted as helpful, the stlf() appears to produce the best results. Outside of that, the seasonal naive shows the best results.
## d. a10
testa10 <- window(a10, start=2006)
ets4 <- ets(testa10)
round(forecast::accuracy(ets4),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0 3.74 2.86 -2.62 14.96 0.81 -0.08sn4 <- snaive(testa10)
round(forecast::accuracy(sn4),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3.21 3.9 3.51 14.43 16.08 1 -0.14st4 <- stlf(testa10)
round(forecast::accuracy(st4),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 shows the best results for the a10 dataset.
## e. h02
testh02 <- window(h02, start=2006)
ets5 <- ets(testh02)
round(forecast::accuracy(ets5),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01 0.05 0.04 0.49 5.04 0.6 -0.33sn5 <- snaive(testh02)
round(forecast::accuracy(sn5),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.03 0.09 0.07 2.87 8.11 1 -0.19st5 <- stlf(testh02)
round(forecast::accuracy(st5),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 shows slightly better results for forecasting the h02 dataset.
## f. usmelec
testUs <- window(usmelec, start=2010)
ets6 <- ets(testUs)
round(forecast::accuracy(ets6),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.07 6.38 5.18 -0.35 1.54 0.7 0.15sn6 <- snaive(testUs)
round(forecast::accuracy(sn6),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.18 9.73 7.39 -0.63 2.2 1 0.2st6 <- stlf(testUs)
round(forecast::accuracy(st6),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.92 6.23 5.05 -0.31 1.5 0.68 0.14## The stlf() model appears to show the best results for forecasting the usmelec dataset.- Use ets() on the following series:
## a. bicoal, chicken, dole, usdeaths, lynx, ibmclose, eggs. Does it always give good forecasts?
## 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## b. Find an example where it does not work well. Can you figure out why?
## The ets() function did not perform overly well on the usdeaths dataset. This can be explained by a number of factors in the dataset. One reason could be that there are not clear patterns that exist in the data set as seen below. The time series range is also fairly short.
usdeaths## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## 1973 9007 8106 8928 9137 10017 10826 11317 10744 9713 9938 9161
## 1974 7750 6981 8038 8422 8714 9512 10120 9823 8743 9129 8710
## 1975 8162 7306 8124 7870 9387 9556 10093 9620 8285 8433 8160
## 1976 7717 7461 7776 7925 8634 8945 10078 9179 8037 8488 7874
## 1977 7792 6957 7726 8106 8890 9299 10625 9302 8314 8850 8265
## 1978 7836 6892 7791 8129 9115 9434 10484 9827 9110 9070 8633
## Dec
## 1973 8927
## 1974 8680
## 1975 8034
## 1976 8647
## 1977 8796
## 1978 9240- 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 code in an ETS model indiciates the ETS classification and notification. The possible inputs are “N” for none, “A” for additive, “M” for multiplicative, or “Z” for automatic selection. In this case, ETS(M,A,M) utilizes the multiplicative and additive processes. This matches the process used within the Holt-Winters' multiplicative method.- Show that the forecast variance for an ETS(A,N,N) model is given by: σ^2[1+α^2(h-1)]
## for the ETS(A,N,N) model, the traditional parameter region is 0<α<1 but the admissible region is 0<α<2 .
## In this case, σ^2[1+α^2(h-1) shows the same forecast variance that can be found in the ETS(A,N,N) model.- Write down 95% prediction intervals for an ETS(A,N,N) model as a function of ℓT, α, h and σ, assuming normally distributed errors.
## Assuming that the forecast errors are normally distributed, a 95% prediction interval for the h-step forecast is ^yT+h|T ± 1.96^σh, where ^σh is an estimate of the standard deviation of the h-step forecast distribution. More generally, a prediction interval can be written as ^yT+h|T ± c^σhChapter 8 - ARIMA Models
- Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.
## a. Explain the differences among these figures. Do they all indicate that the data are white noise?
## The figures show the correlation between a series of lags in the time series. It is a correlation coefficient so the scale is from -1 to 1 as seen on the y-axis. All of the correlations at each lag point lies within the dashed blue line. This indicates white noise within each of the figures.
## b. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?
## The critical value region is dependent on the sample size of numbers used. In this case, we can see the series of numbers increased dramatically from series x1 to x3. As the series number increases, the critical region decreases, indicating that less is required for the correlation to leave the critical value range. - 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")
plot(acf(ibmclose))
plot(pacf(ibmclose))
## The first ACF plot shows a strong trend of autocorrelation in the time series. This can be seen in the pattern shown for the correlations as the lag increases. The correlation value slowly decreases as lag increases.
## The Patrial ACF plot shows a strong outlier at the Lag 1 position with a correlation of almost 1. The values past this point are within the critical region but are still affected by the strong correlation at lag 1.- For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
## a. usnetelec
plot(BoxCox(usnetelec, lambda=0), main="US Net Electricity Generation")
## b. usgdp
plot(BoxCox(usgdp, lambda=1/2), main="US GDP")
## c. mcopper
plot(BoxCox(usgdp, lambda=1/3), main="Monthly Copper Prices")
## d. enplanements
plot(BoxCox(enplanements, lambda=0), main="Monthly US Domestic Enplanements")
## e. visitors
plot(BoxCox(visitors, lambda=0), main="Monthly Australian Overseas Visitors")
- For the enplanements data, write down the differences you chose above using backshift operator notation.
## I chose a Log transformation for enplanements
enp <- plot(BoxCox(usgdp, lambda=0), main="Monthly US Domestic Enplanements")
enp## NULLacf(enplanements)
acf(BoxCox(enplanements, lambda=0))
pacf(enplanements)
## We do not see a major difference in the acf plots here. If we adjust the lambda value further we see similar results
acf(enplanements)
acf(BoxCox(enplanements, lambda=1/3))
pacf(BoxCox(enplanements, lambda=1/3))
## If we incorporate a backshift operater of 2 we see the following:
## pacf(enplanements - lag(enplanements,2)) - 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.
## Importing the retail dataset
retail <- read_excel("retail.xlsx")
retailTS <- ts(retail)
retailTransform <- BoxCox(retailTS, lambda = 0) ## Log transformation- Use R to simulate and plot some data from simple ARIMA models.
## a. 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]
## b. Produce a time plot for the series. How does the plot change as you change ϕ1?
autoplot(y)
## As we decrease ϕ1 we see the data become more condensed and volatile.
y2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y2[i] <- 0.1*y2[i-1] + e[i]
plot(y2)
## c. Write your own code to generate data from an MA(1) model with θ1=0.6 and σ^2=1.
## Select θ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)
## d. 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 and volatile as the theta value increases. As we hit 1.5, the graph becomes uninterpretable to its original form.
## e. 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.3455
## s.e. 0.1139
##
## sigma^2 estimated as 1.298: log likelihood = -154.94, aic = 313.89yAr <- 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)
arima(yAr)##
## Call:
## arima(x = yAr)
##
## Coefficients:
## intercept
## 0.4602
## s.e. 0.1925
##
## sigma^2 estimated as 3.704: log likelihood = -207.37, aic = 418.74## f. Generate data from an AR(2) model with ϕ1=−0.8, ϕ2=0.3 and σ^2=1. (Note that these parameters will give a non-stationary series.)
yAr2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
yAr2[i] <- -0.8*yAr2[i-1] + 0.3*e[i-1] + e[i]
plot(yAr2)
## g. Graph the latter two series and compare them.
plot(yAr)
plot(yAr2)
## The second ARIMA model shows much higher volatility in the data compared to the first Arima model- Consider wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.
## a. 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("Year")
describe(wmurders)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 55 3.39 0.75 3.46 3.4 0.97 2.23 4.49 2.26 -0.14 -1.57
## se
## X1 0.1summary(wmurders)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.233 2.609 3.463 3.393 4.055 4.492wm1 <- 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)
wm1## 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## The ETS autoselection picks ETS(M,N,N)
wm2 <- wmurders %>%
Arima(order=c(0,1,1), seasonal=c(0,1,1)) %>%
residuals() %>% ggtsdisplay()
wm2
wm3 <- arima(wmurders, order=c(3,1,1))
wm3##
## 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 suggests (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## b. Should you include a constant in the model? Explain.
## A constant value could cause unintentional changes to the suggested ARIMA model for the dataset
## c. 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.97acf(wmurders)
pacf(wmurders)
wmurdersTS <- ts(wmurders)
## pacf(wmurders - lag(wmurders,1))
## pacf(wmurdersTS - lag(wmurdersTS,1))
## d. Fit the model using R and examine the residuals. Is the model satisfactory?
## pacf(wmurdersTS - lag(wmurdersTS,1))
## pacf(wmurders - lag(wmurders,1))
## I received an error using the backshift operator format used previously but can assume that if this data is in a time series format the operator would effectively work in this situation.
## e. 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## f. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
plot(forecast(wmurders, h=3))
## g. 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.97autoplot(forecast(auto.arima(wmurders),h=3))
## The two plots give similar forecasts but some differences can be seen in the format and forecast range.- Consider austa, the total international visitors to Australia (in millions) for the period 1980-2015.
## a. Use auto.arima() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
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.57checkresiduals(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## The autoselection provided (0,1,1) as the appropriate model for the dataset
plot(forecast(aus, h=10))
## b. Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part a. Remove the MA term and plot again.
aus2 <- forecast(arima(austa, order=c(0,1,1)),h=10) ## Default drift is FALSE
plot(aus2)
## c. 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) ## include.drift not needed in this case (unused argument)
##plot(aus3)
#Error in arima(austa, order = c(2, 1, 3)) : non-stationary AR part from CSS
## d. 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))
## e. 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)
- For the usgdp series:
## a. if necessary, find a suitable Box-Cox transformation for the data
plot(BoxCox(usgdp, lambda=0), main="US GDP") ## Log transformation
usgdpT <- BoxCox(usgdp, lambda=0)
## b. fit a suitable ARIMA model to the transformed data using auto.arima()
auto.arima(usgdpT)## Series: usgdpT
## 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## c. try some other plausible models by experimenting with the orders chosen;
arima(usgdpT, order=c(0,0,1))##
## Call:
## arima(x = usgdpT, 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.87arima(usgdpT, order=c(1,1,1))##
## Call:
## arima(x = usgdpT, 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## d. choose what you think is the best model and check the residual diagnostics;
gdpA <- auto.arima(usgdpT)
checkresiduals(gdpA)
##
## 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## e. produce forecasts of your fitted model. Do the forecasts look reasonable?
forecast(gdpA, 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.664710autoplot(forecast(gdpA, h=25))
## f. compare the results with what you would obtain using ets() (with no transformation).
ets(usgdpT)## ETS(A,A,N)
##
## Call:
## ets(y = usgdpT)
##
## 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.2552etsGDP <- forecast(ets(usgdpT), h=25)
autoplot(etsGDP)
## The results are very similar in this case between the auto.arima() and ets() functions- Consider austourists, the quarterly visitor nights (in millions) spent by international tourists to Australia for the period 1999–2015.
## a. Describe the time plot.
autoplot(austourists)
## The plot shows a positive trend with very seasonal data
## b. What can you learn from the ACF graph?
acf(austourists)
## The acf graph gives information regarding autocorrelation seen in the data through the different lag values. In this case, many of the values lie outside of the critical region signaling the series is most likely not white noise. (±2/√T where T is the length of time).
## c. What can you learn from the PACF graph?
pacf(austourists)
## The pacf graph provides information regarding partial autocorrelation within the data. It focuses on the correlation of the residuals in the data. 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.
## d. Produce plots of the seasonally differenced data (1−B^4)Yt. What model do these graphs suggest?
## Lagged Differences function
plot(diff(austourists, lag=4, differences=1), xlab="Year", ylab="Total Night Visitors", main="International Tourists to Australia", col="dodgerblue")
## The graph suggests a model that captures the seasonality present within the data. We see a general trend but no strong overall trend can be identified.
## e. 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 also selects a seasonal model with (1−B^4)Yt
## f. Write the model in terms of the backshift operator, then without using the backshift operator.
## (1−B^4)Yt = BY(t-1) + e[i]- Consider usmelec, the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 – June 2013). In general there are two peaks per year: in mid-summer and mid-winter.
## 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="dodgerblue", 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 but for testing purposes we will use a log transformation 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") ## Log transformation
## 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))
## adf.test(diff(usmelec,12)) ## The aTSA pacakge conflicts with other packages in the R Markdown. Results pasted below:
## Augmented Dickey-Fuller Test
##alternative: stationary
##Type 1: no drift no trend
## lag ADF p.value
##[1,] 0 -11.21 0.01
##[2,] 1 -9.01 0.01
##[3,] 2 -6.43 0.01
##[4,] 3 -6.26 0.01
##[5,] 4 -5.58 0.01
##[6,] 5 -4.91 0.01
##Type 2: with drift no trend
## lag ADF p.value
##[1,] 0 -12.53 0.01
##[2,] 1 -10.30 0.01
##[3,] 2 -7.43 0.01
##[4,] 3 -7.35 0.01
##[5,] 4 -6.61 0.01
##[6,] 5 -5.92 0.01
##Type 3: with drift and trend
## lag ADF p.value
##[1,] 0 -12.63 0.01
##[2,] 1 -10.41 0.01
##[3,] 2 -7.54 0.01
##[4,] 3 -7.46 0.01
##[5,] 4 -6.72 0.01
##[6,] 5 -6.04 0.01
##----
##Note: in fact, p.value = 0.01 means p.value <= 0.01
## Log Transform
tsdisplay(diff(melLog, 1))
##adf.test(diff(melLog,12)) ## The aTSA pacakge conflicts with other packages in the R Markdown. Results pasted below:
##Augmented Dickey-Fuller Test
##alternative: stationary
##Type 1: no drift no trend
## lag ADF p.value
##[1,] 0 -10.35 0.01
##[2,] 1 -8.08 0.01
##[3,] 2 -5.90 0.01
##[4,] 3 -5.64 0.01
##[5,] 4 -5.02 0.01
##[6,] 5 -4.40 0.01
##Type 2: with drift no trend
## lag ADF p.value
##[1,] 0 -11.88 0.01
##[2,] 1 -9.48 0.01
##[3,] 2 -7.00 0.01
##[4,] 3 -6.78 0.01
##[5,] 4 -6.08 0.01
##[6,] 5 -5.43 0.01
##Type 3: with drift and trend
## lag ADF p.value
##[1,] 0 -12.25 0.01
##[2,] 1 -9.86 0.01
##[3,] 2 -7.33 0.01
##[4,] 3 -7.13 0.01
##[5,] 4 -6.40 0.01
##[6,] 5 -5.80 0.01
##----
##Note: in fact, p.value = 0.01 means p.value <= 0.01
## The alternate hypothesis of stationary is not met becuase the p-value is below 0.05.
## 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: 24tsdisplay(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 of the Augmented Dickey-Fuller Test
melA <- arima(usmelec, order=c(1,1,1))
tsdisplay(melA$residuals)
## adf.test(melA$residuals) ## The aTSA pacakge conflicts with other packages in the R Markdown. Results pasted below:
##Augmented Dickey-Fuller Test
##alternative: stationary
##Type 1: no drift no trend
## lag ADF p.value
##[1,] 0 -22.2 0.01
##[2,] 1 -17.8 0.01
##[3,] 2 -22.3 0.01
##[4,] 3 -26.0 0.01
##[5,] 4 -14.3 0.01
##[6,] 5 -14.2 0.01
##Type 2: with drift no trend
## lag ADF p.value
##[1,] 0 -22.2 0.01
##[2,] 1 -17.7 0.01
##[3,] 2 -22.3 0.01
##[4,] 3 -26.1 0.01
##[5,] 4 -14.3 0.01
##[6,] 5 -14.2 0.01
##Type 3: with drift and trend
## lag ADF p.value
##[1,] 0 -22.2 0.01
##[2,] 1 -17.7 0.01
##[3,] 2 -22.2 0.01
##[4,] 3 -26.0 0.01
##[5,] 4 -14.3 0.01
##[6,] 5 -14.2 0.01
##----
##Note: in fact, p.value = 0.01 means p.value <= 0.01
## f. Forecast the next 15 years of electricity generation by the U.S. electric industry. Get the latest figures from the EIA to check the accuracy of your forecasts.
## EIA US energy projections: https://www.eia.gov/analysis/projection-data.php
## The forecast is not overly accurate based on various sources from the above link but serves as a good example.
## g. Eventually, the prediction intervals are so wide that the forecasts are not particularly useful. How many years of forecasts do you think are sufficiently accurate to be usable?
## This will vary from dataset to dataset, but in this example, I believe that any forecasts past 5 years should be considered unpredictable with the varying prices, resources, and economy involved in energy consumption.- For the mcopper data:
## a. if necessary, find a suitable Box-Cox transformation for the data
autoplot(mcopper, col="dodgerblue") + ylab("Monthly Copper Prices") + xlab("Year")
mcT <- BoxCox(mcopper, lambda=0) ## Log transform
autoplot(mcT, col="dodgerblue") + ylab("Monthly Copper Prices") + xlab("Year")
## b. fit a suitable ARIMA model to the transformed data using auto.arima()
auto.arima(mcT)## Series: mcT
## 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) selected
## c. try some other plausible models by experimenting with the orders chosen
## Model 2
arima(mcT, order=c(0,0,1))##
## Call:
## arima(x = mcT, 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## Model 3
arima(mcT, order=c(1,1,1))##
## Call:
## arima(x = mcT, 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## d. choose what you think is the best model and check the residual diagnostics;
## Based on the results from the models, I would choose model 1 based on it's low AIC, AICc, and BIC values. In terms of AIC, 0 is an arbitrary number so the negative sign is not a concern.
## e. produce forecasts of your fitted model. Do the forecasts look reasonable?
mc <- arima(mcT)
mcF <- forecast(mc)
mcF## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Feb 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Mar 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Apr 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## May 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Jun 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Jul 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Aug 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Sep 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Oct 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Nov 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Dec 2007 6.730527 5.940282 7.520771 5.521952 7.939101
## Jan 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Feb 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Mar 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Apr 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## May 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Jun 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Jul 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Aug 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Sep 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Oct 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Nov 2008 6.730527 5.940282 7.520771 5.521952 7.939101
## Dec 2008 6.730527 5.940282 7.520771 5.521952 7.939101autoplot(mcF) + ylab("Monthly Copper Prices") + xlab("Year")
## f. compare the results with what you would obtain using ets() (with no transformation).
mc2 <- ets(mcT)
mcF2 <- forecast(mc2)
mcF2## 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.051922autoplot(mcF2) + ylab("Monthly Copper Prices") + xlab("Year")
## The results are fairly similar in the two models. A small difference can be observed in the trend of the forecasted models between the ETS and ARIMA models. The ETS shows a more negative trend compared to the ARIMA.- Choose one of the following seasonal time series: hsales, auscafe, qauselec, qcement, qgas.
## a. Do the data need transforming? If so, find a suitable transformation.
autoplot(hsales) + ylab("Monthly One-Family Home Sales") + xlab("Year")
homeT <- BoxCox(hsales, lambda=0) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
## b. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
## adf.test(homeT)
## According to the Augmented Dickey-Fuller Test, the data is not stationary after the log transform.
## Testing different lambda values
## Y^3
homeT <- BoxCox(hsales, lambda=3) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
## adf.test(homeT) ## The aTSA pacakge conflicts with other packages in the R Markdown. Results pasted below:
## y^.5
homeT <- BoxCox(hsales, lambda=.5) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
## adf.test(homeT)## The aTSA pacakge conflicts with other packages in the R Markdown. Results pasted below:
## Y^-2
homeT <- BoxCox(hsales, lambda=-2) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
## adf.test(homeT)## The aTSA pacakge conflicts with other packages in the R Markdown
## None of the lambda values tested provided stationary data transformation
## c. 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
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 autoselection performed best when compared to the other models. It produced an AIC value of 1630.92. This value was slightly lower than the other produced models.
## d. 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.
salesA <- auto.arima(hsales)
checkresiduals(salesA)
##
## 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: 24acf(hsales)
pacf(hsales)
## There is not white noise present in the model based on the acf values at different lag points. Many of the values lie outside of the critical region.
## e. Forecast the next 24 months of data using your preferred model.
salesF <- forecast(salesA, h=24)
autoplot(salesF)
## f. Compare the forecasts obtained using ets()
salesA2 <- ets(hsales)
salesF2 <- forecast(salesA2, h=24)
autoplot(salesF2)
## The forecast region of the ETS model is much wider than the ARIMA model.- 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?
salesSTLF <- stlf(hsales, method="arima")
salesSTLFf <- forecast(salesSTLF, h=24)
plot(salesSTLFf)
## The STLF() method seems to produce a more focuse forecast and reasonable results. I believe that this method is preffered in this situation.- For your retail time series (Exercise 5 above):
## a. develop an appropriate seasonal ARIMA model
## Ran across errors using the arima() function on the previously imported retail dataset
## retailTS <- ts(retail)
##auto.arima(retailTS)
## b. compare the forecasts with those you obtained in earlier chapters
## The retail dataset produced volatile graphs in previous examples. I believe the best ARIMA model for this situation would capture the seasonal components in the data while transforming the data to prevent autocorrelation.
## c. Obtain up-to-date retail data from the ABS website (Cat 8501.0, Table 11), and compare your forecasts with the actual numbers. How good were the forecasts from the various models?
## Source: https://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/8501.0Oct%202016?OpenDocument
##2015
##August –0.4 0.4 –2.2 –2.1 –0.6 –2.3 –2.9 –0.2 –0.8- Consider sheep, the sheep population of England and Wales from 1867–1939.
## a. Produce a time plot of the time series.
autoplot(sheep) + ylab("Sheep Population") + xlab("Year")
## b. Assume you decide to fit the following model:
## where εt is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?
## White noise series would be equivalent to ARIMA(0,0,1)
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## The autoselection model produces ARIMA(3,1,0)
## c. By examining the ACF and PACF of the differenced data, explain why this model is appropriate.
sheepA <- 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## Base dataset
acf(sheep)
pacf(sheep)
## This model would be appropriate because it would capture the seasonal trends evident in the data. The data does not show strong trends in either direction but the model was selected as the best for capturing 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?
sheepTS <- ts(sheep)
forecast(arima(sheepTS))## 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 similar from the ARIMA model but are subject to some constraints within R.- The annual bituminous coal production in the United States from 1920 to 1968 is in data set bicoal.
## a. Produce a time plot of the data.
autoplot(bicoal) + ylab("Bituminous Coal Production") + xlab("Year")
## b. You decide to fit the following model to the series:
## 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)?
## White noise series would be equivalent to ARIMA(0,0,1)
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 produces a ARIMA(4,0,0) with non-zero mean
## c. Explain why this model was chosen using the ACF and PACF.
coal <- auto.arima(bicoal)
coalF <- forecast(arima(bicoal))
## Base dataset
acf(bicoal)
pacf(bicoal)
## This model would be appropriate because it would capture the seasonal trends evident in the data. The data does not show strong trends in either direction but the model was selected as the best for capturing the data.
## d. 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
## e. Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?
coalTS <- ts(bicoal)
forecast(arima(coalTS))## 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 are subject to some constraints within R.- Before doing this exercise, you will need to install the Quandl package in R using
library(Quandl)
## a. 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.)
## I had issues with importing the data after account creation so I utilized the ibmclose dataset from the previous example
## b. Plot graphs of the data, and try to identify an appropriate ARIMA model.
autoplot(ibmclose) + ylab("Close Price") + xlab("Year")
ibmA <- auto.arima(ibmclose)
ibmA## Series: ibmclose
## ARIMA(0,1,0)
##
## sigma^2 estimated as 52.62: log likelihood=-1251.37
## AIC=2504.74 AICc=2504.75 BIC=2508.64## Autoselection chose ARIMA(0,1,0) as the appropriate model
## c. Do residual diagnostic checking of your ARIMA model. Are the residuals white noise?
checkresiduals(ibmA)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)
## Q* = 14.151, df = 10, p-value = 0.1662
##
## Model df: 0. Total lags used: 10## The residuals do show signs of white noise based on the ACF graph. Some values go beyond the large critical region
## d. Use your chosen ARIMA model to forecast the next four years.
forecast(arima(ibmclose))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 370 478.4688 370.6839 586.2538 313.626 643.3117
## 371 478.4688 370.6839 586.2538 313.626 643.3117
## 372 478.4688 370.6839 586.2538 313.626 643.3117
## 373 478.4688 370.6839 586.2538 313.626 643.3117
## 374 478.4688 370.6839 586.2538 313.626 643.3117
## 375 478.4688 370.6839 586.2538 313.626 643.3117
## 376 478.4688 370.6839 586.2538 313.626 643.3117
## 377 478.4688 370.6839 586.2538 313.626 643.3117
## 378 478.4688 370.6839 586.2538 313.626 643.3117
## 379 478.4688 370.6839 586.2538 313.626 643.3117plot(forecast(arima(ibmclose)))
## e. Now try to identify an appropriate ETS model.
ibm2 <- ets(ibmclose, model="ZZZ", lambda="auto")
ibm2## ETS(A,N,N)
##
## Call:
## ets(y = ibmclose, model = "ZZZ", lambda = "auto")
##
## Box-Cox transformation: lambda= 1.9999
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 109274.6793
##
## sigma: 3191.485
##
## AIC AICc BIC
## 8139.441 8139.507 8151.173autoplot(ibm2)
## The model selected is ETS(A,N,N)
## f. Do residual diagnostic checking of your ETS model. Are the residuals white noise?
checkresiduals(ibm2)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 16.112, df = 8, p-value = 0.04081
##
## Model df: 2. Total lags used: 10## The residual check does show signs of white noise in the ACF graph.
## g. Use your chosen ETS model to forecast the next four years.
ibm2F <- forecast.ets(ibm2)
ibm2F## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 370 356.9995 345.3475 368.2830 339.0173 374.1185
## 371 356.9995 340.4052 372.8561 331.2844 380.9830
## 372 356.9995 336.5634 376.3275 325.2259 386.1678
## 373 356.9995 333.2903 379.2294 320.0292 390.4853
## 374 356.9995 330.3797 381.7678 315.3799 394.2500
## 375 356.9995 327.7261 384.0482 311.1168 397.6229
## 376 356.9995 325.2667 386.1334 307.1442 400.6995
## 377 356.9995 322.9607 388.0642 303.3998 403.5421
## 378 356.9995 320.7798 389.8689 299.8405 406.1938
## 379 356.9995 318.7033 391.5682 296.4347 408.6860autoplot(ibm2F) + ylab("Stock Close") + xlab("Year")
## h. Which of the two models do you prefer?
## In this case, I prefer the ETS model forecast becuase the forecast range is smaller and acts as a better predictor.