#packages needed
library(readr)
library(ggplot2)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(fpp2)
## Loading required package: fma
##
## Attaching package: 'fma'
## The following object is masked from 'package:GGally':
##
## pigs
## Loading required package: expsmooth
library(seasonal)
library(apricom)
library(psych)
##
## Attaching package: 'psych'
## The following object is masked from 'package:seasonal':
##
## outlier
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:GGally':
##
## nasa
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# CHAPTER 7
#Q1
#a
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.4
#break into 4 quarters
pigsSES2 <- ses(pigs, h=4)
#b.
pigsSES$upper[1,"95%"]
## 95%
## 119020.8
pigsSES$lower[1,"95%"]
## 95%
## 78611.97
res <- sd(pigsSES$residuals)
res
## [1] 10273.69
pigsSES$mean[1] + 1.96*res
## [1] 118952.8
pigsSES$mean[1] - 1.96*res
## [1] 78679.97
#2
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.2
forecast(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 forecast seems to be slightly different forecast when using the forecast function
#3
SES <- function(pars = c(alpha, l0), y){
# change the first argument as vector of alpha and l0, rather than separate alpha and l0 because optim function wants to take a function that requires vector as its first argument as fn argument.
error <- 0
SSE <- 0
alpha <- pars[1]
l0 <- pars[2]
y_hat <- l0
for(index in 1:length(y)){
error <- y[index] - y_hat
SSE <- SSE + error^2
y_hat <- alpha*y[index] + (1 - alpha)*y_hat
}
return(SSE)
}
opt_SES_pigs <- optim(par = c(0.5, pigs[1]), y = pigs, fn = SES)
writeLines(paste(
"Optimal parameters for the result of SES function: ",
"\n",
as.character(opt_SES_pigs$par[1]),
", ",
as.character(opt_SES_pigs$par[2]),
sep = ""
))
## Optimal parameters for the result of SES function:
## 0.299008094014243, 76379.2653476235
writeLines(paste(
"Parameters got from the result of ses function: ",
"\n",
as.character(pigsSES2$model$par[1]),
", ",
as.character(pigsSES2$model$par[2]),
sep = ""
))
## Parameters got from the result of ses function:
## 0.297148833372095, 77260.0561458528
# when comparing the 2, you can see that they are almost the same no matter the function being used.
#4
# modify SES function to find the optimal values of alpha and l0, and produce the next observation forecast.
SES <- function(init_pars, data){
# init_pars is c(init_alpha, init_l0)
# make next observation forecast variable
fc_next <- 0
# make SSE function to get SSE if parameters and data are given.
SSE <- function(pars, data){
error <- 0
SSE <- 0
alpha <- pars[1]
l0 <- pars[2]
y_hat <- l0
for(index in 1:length(data)){
error <- data[index] - y_hat
SSE <- SSE + error^2
y_hat <- alpha*data[index] + (1 - alpha)*y_hat
}
# use superassignment to make forecast value possible to use outside SSE function.
fc_next <<- y_hat
return(SSE)
}
# use optim function to get optimal values of alpha and l0.
optim_pars <- optim(par = init_pars, data = data, fn = SSE)
# return results
return(list(
Next_observation_forecast = fc_next,
alpha = optim_pars$par[1],
l0 = optim_pars$par[2]
))
}
# compare the result using pigs data
SES(c(0.5, pigs[1]), pigs)
## $Next_observation_forecast
## [1] 98814.63
##
## $alpha
## [1] 0.2990081
##
## $l0
## [1] 76379.27
print("Next observation forecast by ses function")
## [1] "Next observation forecast by ses function"
pigsSES2$mean[1]
## [1] 98816.41
print("alpha calculated by ses function")
## [1] "alpha calculated by ses function"
pigsSES2$model$par[1]
## alpha
## 0.2971488
print("l0 calculated by ses function")
## [1] "l0 calculated by ses function"
pigsSES2$model$par[2]
## l
## 77260.06
# compare the result using ausbeer data
ses_ausbeer <- ses(ausbeer, h = 1)
SES(c(0.5, ausbeer[1]), ausbeer)
## $Next_observation_forecast
## [1] 421.3166
##
## $alpha
## [1] 0.1488036
##
## $l0
## [1] 259.6585
print("Next observation forecast by ses function")
## [1] "Next observation forecast by ses function"
ses_ausbeer$mean[1]
## [1] 421.3136
print("alpha calculated by ses function")
## [1] "alpha calculated by ses function"
ses_ausbeer$model$par[1]
## alpha
## 0.1489004
print("l0 calculated by ses function")
## [1] "l0 calculated by ses function"
ses_ausbeer$model$par[2]
## l
## 259.6592
# SES function worked similar to ses function.
#Question 5
#a
plot(books, ylab="Book sales", xlab="Number of Days", main="Book Sales over Time")
#both types of books show sales that are seasonal and actually tend to act as inverses of each other.
#b.
books <- data.frame(books)
hardTS <- ts(books$Hardcover)
paperTS <- ts(books$Paperback)
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
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
autoplot(hF)
autoplot(pF)
#C
round(forecast::accuracy(pF),2)[2]
## [1] 33.64
round(forecast::accuracy(hF),2)[2]
## [1] 31.93
## Hardcover RMSE = 33.64
## Paperback RMSE = 31.93
#6
#a
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.3334
h2 <- holt(paperTS, h=4)
h2
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 209.4668 166.6035 252.3301 143.9130 275.0205
## 32 210.7177 167.8544 253.5811 145.1640 276.2715
## 33 211.9687 169.1054 254.8320 146.4149 277.5225
## 34 213.2197 170.3564 256.0830 147.6659 278.7735
#b
## 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.
#Looking at the data you can see that the holts method producted the best overal results with the lowest RMSE value
#d
## Upper Hardcover
h1$mean[1] + round(forecast::accuracy(h1),2)*1.96
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 249.8995 303.4663 295.5675 246.0383 274.0075 251.5263 250.1151
## Lower Hardcover
h1$mean[1] - round(forecast::accuracy(h1),2)*1.96
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 250.4483 196.8815 204.7803 254.3095 226.3403 248.8215 250.2327
## Upper Paperback
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
#Q7
## 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.6195
lambda <- 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")
#Q8
install.packages("readxl")
## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)
library(readxl)
retailXL <- read_excel("retail.xlsx")
## New names:
## * `` -> ...1
#a.
retailTS <- ts(retailXL[,20], frequency=12, start=c(1982,4))
## Warning in data.matrix(data): NAs introduced by coercion
autoplot(retailTS)
#multiplicative seasonality will be necessary for this series becasue the auto graph has an increasing trend
#b.
## Damped = FALSE
holtFit1 <- hw(retailTS, seasonal="multiplicative", damped=FALSE, h=30)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Missing values encountered. Using longest contiguous portion of time series
autoplot(retailTS) + autolayer(holtFit1) + ylab("Retail Sales") + xlab("Year")
## Damped = TRUE
holtFit2 <- hw(retailTS, seasonal="multiplicative", damped=TRUE, h=30)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Missing values encountered. Using longest contiguous portion of time series
autoplot(retailTS) + autolayer(holtFit2) + ylab("Retail Sales") + xlab("Year")
#c.
holtFit1$mean[1]
## [1] 392.1741
holtFit2$mean[1]
## [1] 394.1059
round(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.07
round(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
#when looking at the values above the damped version has better results. Along with that it also has slightly lower rmse
#d.
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.
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
#rmse of 27.31
#9.
lambdaTS <- BoxCox.lambda(retailTS)
bc1 <- BoxCox(retailTS, lambdaTS)
bc1
## Jan Feb Mar Apr May Jun Jul Aug
## 1982 NA 4.171794 4.165599 4.122013 4.140355
## 1983 4.285306 4.246614 4.168702 4.154121 4.111036 4.126364 4.111036 4.164562
## 1984 4.213951 4.288901 4.220778 4.261562 4.207068 4.244726 4.191123 4.205091
## 1985 4.349467 4.308405 4.188100 4.252250 4.258779 4.277161 4.239040 4.280791
## 1986 4.480012 4.402869 4.290692 4.355293 4.387039 4.415323 4.361903 4.363547
## 1987 4.575353 4.458994 4.394990 4.402869 4.410674 4.410674 4.396572 4.483584
## 1988 4.574073 4.553985 4.477144 4.535421 4.499822 4.515066 4.494209 4.508863
## 1989 4.578543 4.565707 4.461193 4.504702 4.474266 4.494209 4.562467 4.563116
## 1990 4.728110 4.767675 4.655198 4.704320 4.703223 4.704320 4.691054 4.679271
## 1991 4.742381 4.681528 4.577906 4.638757 4.643487 4.628016 4.617743 4.634596
## 1992 4.760013 4.754865 4.682092 4.757443 4.726509 4.729708 4.663875 4.636382
## 1993 4.698266 4.632209 4.528014 4.481443 4.546079 4.492800 4.538098 4.576630
## 1994 4.737651 4.805454 4.646430 4.774259 4.647018 4.655198 4.654616 4.689939
## 1995 4.743429 4.776275 4.669039 4.756928 4.756928 4.751762 4.717380 4.803507
## 1996 5.006971 5.036951 4.952566 4.998147 4.975867 4.957441 4.921857 4.917619
## 1997 5.000074 5.011158 4.877582 4.970701 4.939840 4.954601 4.903046 4.939840
## 1998 5.012296 4.972692 4.871767 4.901748 4.782293 4.755897 4.766657 4.793214
## 1999 4.819909 4.787274 4.660415 4.743952 4.707604 4.722762 4.703772 4.673033
## 2000 4.741857 4.697713 4.618958 4.692168 4.723835 4.703772 4.736068 4.875797
## 2001 4.986887 4.968306 4.883358 4.984930 4.844387 4.808849 4.835526 4.860012
## 2002 4.987668 5.027698 4.829881 4.941080 4.933194 4.901748 4.888214 4.949710
## 2003 5.029186 5.014946 4.828463 4.918469 5.020972 5.032153 4.989619 5.046102
## 2004 5.174021 5.153496 5.031041 5.075427 5.157005 5.147720 5.118223 5.180235
## 2005 5.198904 5.136374 5.043550 5.097615 5.109839 5.097615 5.104768 5.204934
## 2006 5.319718 5.234730 5.109839 5.190692 5.177444 5.170270 5.156687 5.243950
## 2007 5.243377 5.195872 5.085208 5.183941 5.174021 5.174332 5.185788 5.266558
## 2008 5.309463 5.237911 5.209428 5.256462 5.248807 5.264605 5.232118 5.273222
## 2009 5.385715 5.353505 5.210622 5.305751 5.316576 5.331136 5.285571 5.378206
## 2010 5.552045 5.524641 5.376501 5.427775 5.487176 5.476025 5.484616 5.561627
## 2011 5.622297 5.526064 5.381606 5.432571 5.410180 5.369153 5.364960 5.351497
## 2012 5.502996 5.426629 5.292083 5.366936 5.429148 5.402185 5.409946 5.465595
## 2013 5.550273 5.415778 5.306813 5.435072 5.430519 5.416242 5.411116 5.413683
## 2014 5.562793
## Sep Oct Nov Dec
## 1982 4.122013 4.137147 4.191123 4.192129
## 1983 4.169734 4.160401 4.115444 4.145676
## 1984 4.220778 4.210024 4.280791 4.290692
## 1985 4.311904 4.311904 4.366008 4.416095
## 1986 4.400513 4.379817 4.426072 4.457525
## 1987 4.471378 4.489267 4.498422 4.504702
## 1988 4.503310 4.481443 4.510246 4.513004
## 1989 4.584889 4.581087 4.624404 4.641717
## 1990 4.696053 4.570865 4.646430 4.647018
## 1991 4.659258 4.624404 4.674171 4.678706
## 1992 4.634000 4.647605 4.677008 4.667894
## 1993 4.605496 4.642898 4.694390 4.698818
## 1994 4.693835 4.701574 4.654616 4.650533
## 1995 4.817038 4.815599 4.921012 4.934444
## 1996 4.937354 4.900448 4.967906 4.944378
## 1997 4.948073 4.957036 4.975867 4.952566
## 1998 4.722226 4.715759 4.781793 4.751762
## 1999 4.716300 4.738704 4.750725 4.744475
## 2000 4.873560 4.927755 4.993506 4.894790
## 2001 4.855447 4.857733 4.929432 4.869069
## 2002 4.922703 4.912932 4.986104 4.927335
## 2003 5.058746 5.041722 5.116219 5.070847
## 2004 5.131791 5.152856 5.178686 5.115549
## 2005 5.142877 5.158913 5.255332 5.236756
## 2006 5.225992 5.211815 5.178065 5.179306
## 2007 5.266558 5.235310 5.268229 5.253635
## 2008 5.263486 5.212708 5.295588 5.249091
## 2009 5.383059 5.397920 5.458125 5.430291
## 2010 5.522196 5.542150 5.549287 5.504874
## 2011 5.344935 5.374793 5.373080 5.336529
## 2012 5.436433 5.442079 5.433254 5.444999
## 2013 5.455030 5.418098 5.465595 5.486324
## 2014
lambdaTS <- BoxCox.lambda(testRetail)
bc1 <- BoxCox(testRetail, lambdaTS)
bc1
## Jan Feb Mar Apr May Jun Jul Aug
## 2011 5061.528 4285.012 3340.600 3646.792 3508.894 3269.857 3246.380 3172.157
## 2012 4117.653 3609.668 2864.726 3257.422 3625.360 3460.951 3507.481 3860.340
## 2013 4468.125 3542.860 2937.995 3662.534 3633.928 3545.695 3514.548 3530.111
## 2014 4565.935
## Sep Oct Nov Dec
## 2011 3136.611 3301.703 3292.000 3091.674
## 2012 3671.129 3707.010 3651.083 3725.711
## 2013 3790.669 3557.042 3860.340 4000.857
## 2014
retailSTL <- 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 performs poorly, the previous test was much better.
#10
#a.
autoplot(ukcars) + ylab("Vehicle Production") + xlab("Year")
#shows an overall positive trend in vechile production, with some sesaonality.
#b
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.309325330
summary(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.
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.
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.
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
#the ets function chose ETS (ANA) as the correct version
#f.
## 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)
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## ME RMSE MAE MPE MAPE
## Training set NaN NaN NaN NaN NaN
#g.
## 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
#overal the results are similr. i picked model 1 becasue the ets function picked the best function by itself which removes user error.
#h
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.071449002
summary(mForecast1$residuals)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -74.4850 -17.8438 2.2317 -0.3413 14.8230 53.6357
#11.
#a.
autoplot(visitors) + ylab("Short-term Overseas Visitors") + xlab("Year")
#the plot has a positive trend with slight seasonality.
#b.
## Testing set
testVisitor <- window(visitors, start=2004)
## Forecasting one year (12 months)
hw(testVisitor, seasonal="multiplicative", damped=FALSE, h=12)
## Warning in ets(x, "MAM", alpha = alpha, beta = beta, gamma = gamma, phi = phi, :
## Seasonal component could not be estimated
## 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.
#the reason we need multiplicative seasonality is because the trend is increasing. this method would capture all trends that additive would miss.
#d.
## 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.8476
summary(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.
#Looking at the data it seems that the seasonal naive model is the best. The sample size is small though which can limit the accuracy. it does not pass the residual test but if the data set was larger it might.
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: 3
summary(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.
tsCV(tsCV(visitors, forecastfunction, h = 1, window = NULL, xreg = NULL,initial = 0))
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1985 NA NA NA NA NA NA NA NA
## 1986 NA NA NA NA NA NA NA NA NA NA NA NA
## 1987 NA NA NA NA NA NA NA NA NA NA NA NA
## 1988 NA NA NA NA NA NA NA NA NA NA NA NA
## 1989 NA NA NA NA NA NA NA NA NA NA NA NA
## 1990 NA NA NA NA NA NA NA NA NA NA NA NA
## 1991 NA NA NA NA NA NA NA NA NA NA NA NA
## 1992 NA NA NA NA NA NA NA NA NA NA NA NA
## 1993 NA NA NA NA NA NA NA NA NA NA NA NA
## 1994 NA NA NA NA NA NA NA NA NA NA NA NA
## 1995 NA NA NA NA NA NA NA NA NA NA NA NA
## 1996 NA NA NA NA NA NA NA NA NA NA NA NA
## 1997 NA NA NA NA NA NA NA NA NA NA NA NA
## 1998 NA NA NA NA NA NA NA NA NA NA NA NA
## 1999 NA NA NA NA NA NA NA NA NA NA NA NA
## 2000 NA NA NA NA NA NA NA NA NA NA NA NA
## 2001 NA NA NA NA NA NA NA NA NA NA NA NA
## 2002 NA NA NA NA NA NA NA NA NA NA NA NA
## 2003 NA NA NA NA NA NA NA NA NA NA NA NA
## 2004 NA NA NA NA NA NA NA NA NA NA NA NA
## 2005 NA NA NA NA
#12.
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.
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.614039
f2 <- 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.
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.04
round(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.61
summary(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 MASE
## Training set 0.0006237224 0.07829922 0.05603263 -0.08782872 3.606505 0.5485176
## ACF1
## Training set -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.614039
summary(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 ACF1
## Training set 0.03419651 0.1338564 0.1021528 2.338509 6.768085 1 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 becasue of missing values in the data set. this could definately cause issues down the line for modeling.
#13.
## 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.07
sn1 <- 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.12
st1 <- 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 looks like the best fit
#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.02
sn2 <- 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.3
st2 <- 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
#STLF() still looks like the best option
#c.
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.47
sn3 <- 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.74
st3 <- 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
#this dole test set has interesting results. if negative value for ME are helpful then the stlf is the best. otherwise the seasonal naive method shows the best results.
#d
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.08
sn4 <- 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.14
st4 <- 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
#ets is the best option
#e.
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.33
sn5 <- 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.19
st5 <- 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
# ets is slightly better.
#f.
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.15
sn6 <- 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.2
st6 <- 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
#stlf is the best option
#14.
#a.
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
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
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
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
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
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.
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.
#the ets function wasnt great for the usdeaths data set. this is probably becasue of the number of factors in the set. also there isnt a clear pattern that exist.
usdeaths
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1973 9007 8106 8928 9137 10017 10826 11317 10744 9713 9938 9161 8927
## 1974 7750 6981 8038 8422 8714 9512 10120 9823 8743 9129 8710 8680
## 1975 8162 7306 8124 7870 9387 9556 10093 9620 8285 8433 8160 8034
## 1976 7717 7461 7776 7925 8634 8945 10078 9179 8037 8488 7874 8647
## 1977 7792 6957 7726 8106 8890 9299 10625 9302 8314 8850 8265 8796
## 1978 7836 6892 7791 8129 9115 9434 10484 9827 9110 9070 8633 9240
#15
# The 3 letter code in an ets model describes the classifcation and notification. There are 4 different inputs possible. N = none, A = Additive , M= multiplicative, and z= automatic. For this case the ETS MAM utlizes the multiplicative and additive process. this matches the process used within the holts winters multiplicative met hod.
#16
#For the ETS(ANN) model the traditonal paramaters is 0<A<1 but the amissable region is 0<a<2. so for this case o^2[1+a^2(h-1)] which will the the same forecast variance that can be found in the ets ANN model.
#17
#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^σh
#############################
#CHAPTER 8 PROBLEMS
#Q1
#a. The figures show a correlation between a series of lags in the time series. It is a correlation coefficient therefore 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. The CV region is dependent on the numbers used in the sample size. For this case the series of numbers increased rapidly between the series. As the series number increase the critical region decreases. This indicates that less is required for th correlation to leave the CV range
#2.
autoplot(ibmclose) + ylab("Close Price") + xlab("Day")
plot(acf(ibmclose))
plot(pacf(ibmclose))
#the first plot shows a strong trend in the timeseries for auto correlation. This can be seen by the pattern when looking at the correlation as lag increases. the correlation slowly goes down as lag goes up.
# the patrial acf plot shows the outliers at the first lag position that has a correlation of 1. any values past this point are still critical by are highly affected by correlation
#3.
#a.
plot(BoxCox(usnetelec, lambda=0), main="US Net Electricity Generation")
#b.
plot(BoxCox(usgdp, lambda=1/2), main="US GDP")
#c.
plot(BoxCox(usgdp, lambda=1/3), main="Monthly Copper Prices")
#d.
plot(BoxCox(enplanements, lambda=0), main="Monthly US Domestic Enplanements")
#e.
plot(BoxCox(visitors, lambda=0), main="Monthly Australian Overseas Visitors")
#4.
enp <- plot(BoxCox(usgdp, lambda=0), main="Monthly US Domestic Enplanements")
enp
## NULL
acf(enplanements)
acf(BoxCox(enplanements, lambda=0))
pacf(enplanements)
# no big differences in the plots, next we will tweak the lambda values.
acf(enplanements)
acf(BoxCox(enplanements, lambda=1/3))
pacf(BoxCox(enplanements, lambda=1/3))
#if we add in a backshift of 2 we see the following
#pacf(enplanements - lag(enplanements,2))
#Q5
## Importing the retail dataset
retail <- read_excel("retail.xlsx")
## New names:
## * `` -> ...1
retailTS <- ts(retail)
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
## Warning in data.matrix(data): NAs introduced by coercion
retailTransform <- BoxCox(retailTS, lambda = 0) ## Log transformation
#6.
#a
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
#b
autoplot(y)
y2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y2[i] <- 0.1*y2[i-1] + e[i]
plot(y2)
#c.
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
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)
#as theta goes up the plots will become less condesne and volitile. when we hit 1.5 the graph becomes lost.
#e.
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.3001
## s.e. 0.1223
##
## sigma^2 estimated as 1.495: log likelihood = -162, aic = 328
yAr <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
yAr[i] <- 0.6*yAr[i-1] + 0.6*e[i-1] + e[i]
plot(yAr)
arima(yAr)
##
## Call:
## arima(x = yAr)
##
## Coefficients:
## intercept
## 0.5319
## s.e. 0.1608
##
## sigma^2 estimated as 2.587: log likelihood = -189.42, aic = 382.83
#f.
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
plot(yAr)
plot(yAr2)
#its clear that the second graph shows higher volitility
#7.
#a.
autoplot(wmurders) + ylab("Female Murder Rate") + xlab("Year")
summary(wmurders)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.233 2.609 3.463 3.393 4.055 4.492
wm1 <- 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
#ets picked MNN
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
## arima picks (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.
#a value that is constant could possibly casuse unintentional changes ot the ARIMA model for the data set.
#c.
wmAuto <- auto.arima(wmurders)
wmAuto
## Series: wmurders
## ARIMA(1,2,1)
##
## Coefficients:
## ar1 ma1
## -0.2434 -0.8261
## s.e. 0.1553 0.1143
##
## sigma^2 estimated as 0.04632: log likelihood=6.44
## AIC=-6.88 AICc=-6.39 BIC=-0.97
acf(wmurders)
pacf(wmurders)
wmurdersTS <- ts(wmurders)
#pacf(wmurders - lag(wmurders,1))
# pacf(wmurders - lag(wmurders,1))
# pacf(wmurdersTS - lag(wmurdersTS,1))
# d.
# pacf(wmurdersTS - lag(wmurdersTS,1))
# pacf(wmurders - lag(wmurders,1))
#I was unable to get the backshift operator to format correctly but based on my knowledge i can assum that the data given is set in a time series format, which would let the operator work as should.
#e
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.
plot(forecast(wmurders, h=3))
#g.
auto.arima(wmurders)
## Series: wmurders
## ARIMA(1,2,1)
##
## Coefficients:
## ar1 ma1
## -0.2434 -0.8261
## s.e. 0.1553 0.1143
##
## sigma^2 estimated as 0.04632: log likelihood=6.44
## AIC=-6.88 AICc=-6.39 BIC=-0.97
autoplot(forecast(auto.arima(wmurders),h=3))
# the two plots give similar forecasts but have slight differences that could be seen
#8.
#a.
aus <- auto.arima(austa)
aus
## Series: austa
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.3006 0.1735
## s.e. 0.1647 0.0390
##
## sigma^2 estimated as 0.03376: log likelihood=10.62
## AIC=-15.24 AICc=-14.46 BIC=-10.57
checkresiduals(aus)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1) with drift
## Q* = 2.297, df = 5, p-value = 0.8067
##
## Model df: 2. Total lags used: 7
#autoselection picked 0,1,1
plot(forecast(aus, h=10))
#b.
aus2 <- forecast(arima(austa, order=c(0,1,1)),h=10) ## Default drift is FALSE
plot(aus2)
#c.
##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.
aus4 <- forecast(arima(austa, order=c(0,0,1)),h=10)
plot(aus4)
## Removing MA term
plot(forecast(arima(austa),h=10))
#e.
aus5 <- forecast(arima(austa, order=c(0,2,1)),h=10)
plot(aus5)
#9.
#a.
plot(BoxCox(usgdp, lambda=0), main="US GDP") ## Log transformation
usgdpT <- BoxCox(usgdp, lambda=0)
#b.
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.
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.87
arima(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.
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.
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.664710
autoplot(forecast(gdpA, h=25))
#f.
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.2552
etsGDP <- forecast(ets(usgdpT), h=25)
autoplot(etsGDP)
#overal results are similar between the two functions
#10
#a.
autoplot(austourists)
#shows postive trend with very seasonal data.
#b.
acf(austourists)
#this graph above gives us info regarding auto correlation seen in the data through the different values. for this issue, alof of value fall outside of the critical region which tells me that it is not white noise.
#c.
pacf(austourists)
#the pacf graph show info on the partial auto correlation in the data set. it will focus on the residual in the data. it shows that there is not much white noice in the data.
#d.
plot(diff(austourists, lag=4, differences=1), xlab="Year", ylab="Total Night Visitors", main="International Tourists to Australia", col="dodgerblue")
#the graph suggest a model that shows the seasonality present within the data. overal there is a general trend but nothing to crazy
#e.
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 season model
#f.
## (1−B^4)Yt = BY(t-1) + e[i]
#11.
#a.
mel1 <- ma(usmelec, order=12)
plot(mel1, col="dodgerblue", main="Electricity Monthly Total Net Generation", xlab="Year", ylab="Electricity Generation (billions kWh")
#b.
#the data seems to need no transformation but for the testing purpose we will use a log transformation to examine.
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.
#data as small seasonal patterns
tsdisplay(diff(usmelec, 1))
tsdisplay(diff(melLog, 1))
#d.
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 auto selection model is better than AIC
#e.
autoMel <- auto.arima(usmelec)
checkresiduals(autoMel)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,1)[12] with drift
## Q* = 42.725, df = 19, p-value = 0.001413
##
## Model df: 5. Total lags used: 24
tsdisplay(autoMel$residuals)
melARIMAauto <- arima(usmelec)
autoMelF <- forecast(melARIMAauto)
autoMelF
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2013 259.5576 171.3822 347.733 124.705 394.4102
## Aug 2013 259.5576 171.3822 347.733 124.705 394.4102
## Sep 2013 259.5576 171.3822 347.733 124.705 394.4102
## Oct 2013 259.5576 171.3822 347.733 124.705 394.4102
## Nov 2013 259.5576 171.3822 347.733 124.705 394.4102
## Dec 2013 259.5576 171.3822 347.733 124.705 394.4102
## Jan 2014 259.5576 171.3822 347.733 124.705 394.4102
## Feb 2014 259.5576 171.3822 347.733 124.705 394.4102
## Mar 2014 259.5576 171.3822 347.733 124.705 394.4102
## Apr 2014 259.5576 171.3822 347.733 124.705 394.4102
## May 2014 259.5576 171.3822 347.733 124.705 394.4102
## Jun 2014 259.5576 171.3822 347.733 124.705 394.4102
## Jul 2014 259.5576 171.3822 347.733 124.705 394.4102
## Aug 2014 259.5576 171.3822 347.733 124.705 394.4102
## Sep 2014 259.5576 171.3822 347.733 124.705 394.4102
## Oct 2014 259.5576 171.3822 347.733 124.705 394.4102
## Nov 2014 259.5576 171.3822 347.733 124.705 394.4102
## Dec 2014 259.5576 171.3822 347.733 124.705 394.4102
## Jan 2015 259.5576 171.3822 347.733 124.705 394.4102
## Feb 2015 259.5576 171.3822 347.733 124.705 394.4102
## Mar 2015 259.5576 171.3822 347.733 124.705 394.4102
## Apr 2015 259.5576 171.3822 347.733 124.705 394.4102
## May 2015 259.5576 171.3822 347.733 124.705 394.4102
## Jun 2015 259.5576 171.3822 347.733 124.705 394.4102
melA <- arima(usmelec, order=c(1,1,1))
tsdisplay(melA$residuals)
#12
#a.
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.
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
#c.
#model2
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.
#based on the results from the models i would have picked model 1 since it has the lowest AIC and BIC values.
#e
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.939101
autoplot(mcF) + ylab("Monthly Copper Prices") + xlab("Year")
#f
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.051922
autoplot(mcF2) + ylab("Monthly Copper Prices") + xlab("Year")
#the models are pretty much the same. the forecasted models are slightly diffrent but the ets shows a more negative trend compared to the arima
#13
#a.
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.
#adf.test(homeT)
homeT <- BoxCox(hsales, lambda=3) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
homeT <- BoxCox(hsales, lambda=.5) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
homeT <- BoxCox(hsales, lambda=-2) ## Log Transform
autoplot(homeT) + ylab("Monthly One-Family Home Sales") + xlab("Year")
#c.
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
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
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 auto selection was the best method. the AIC was lower than the others.
#d.
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: 24
acf(hsales)
pacf(hsales)
#there is no white noise. Most values are out of the creitical region
#e.
salesF <- forecast(salesA, h=24)
autoplot(salesF)
#f.
salesA2 <- ets(hsales)
salesF2 <- forecast(salesA2, h=24)
autoplot(salesF2)
# the forecast is wider via ETS than the ARIMA
#14
salesSTLF <- stlf(hsales, method="arima")
salesSTLFf <- forecast(salesSTLF, h=24)
plot(salesSTLFf)
#stlf looks to produce a more focused forecast with good results.
#15
#a.
#retailTS<- ts(retail)
#auto.arima(retailTS)
#b.
#the retail has more volatile graphs then the previous. I think the ARIMA model was best for this data. It caputures the seasonal components in the data. The goal would be to tranform data to prevent autocorrelation.
#c.
## 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
#16.
#a.
autoplot(sheep) + ylab("Sheep Population") + xlab("Year")
#b.
#white noise series would be similar 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
#auto selection model picked ARIMA 3,1,0
#c.
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
acf(sheep)
pacf(sheep)
#the model is appropriate because it would capture the seasonality that is in the data. the data is pretty week in trends in either direction but the model was best for captureing data.
#d.
## 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.
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
#results are pretty close but the arima function is limited by constraints in R.
#17
#a.
autoplot(bicoal) + ylab("Bituminous Coal Production") + xlab("Year")
#b.
## 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
#auto selection picked 4,0,0
#c.
coal <- auto.arima(bicoal)
coalF <- forecast(arima(bicoal))
acf(bicoal)
pacf(bicoal)
# this model captures the seasonal trends nor shows strong trends in either directions.
#d.
## 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.
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
#results are similar
#18
install.packages("Quandl")
## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)
library(Quandl)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
#a
y<-Quandl(raw())
#had issues with quandl so i used the data from IBM close from previous problem
#b.
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
#auto arima picked 0,1,0
#c.
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 show signs of white noise and some values go past the critical range
#d.
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.3117
plot(forecast(arima(ibmclose)))
#e,
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.173
autoplot(ibm2)
#the ets is ANN
#f
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
#residuals show white noise
#g.
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.6860
autoplot(ibm2F) + ylab("Stock Close") + xlab("Year")
#h for this case i would pick the ets model since the forecase range is smaller and looks like a better predictor.