Hyndman Exercises Chapter 7 & 8
Daniel Dychala
4/25/2021
knitr::opts_chunk$set(echo = TRUE)
library(seasonal)## Warning: package 'seasonal' was built under R version 4.0.5library(fpp2)## Warning: package 'fpp2' was built under R version 4.0.5## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## -- Attaching packages ---------------------------------------------- fpp2 2.4 --## v ggplot2 3.3.3 v fma 2.4
## v forecast 8.14 v expsmooth 2.3## Warning: package 'fma' was built under R version 4.0.5## Warning: package 'expsmooth' was built under R version 4.0.5## library(forecast)
library(knitr)
library(urca)
library(Quandl)## Warning: package 'Quandl' was built under R version 4.0.5## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericChapter 7 (7) Exercises
Question 1: Part A (7.1.A)
ses.pigs <- ses(pigs, h = 4)
summary(ses.pigs)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = pigs, h = 4)
##
## Smoothing parameters:
## alpha = 0.2971
##
## Initial states:
## l = 77260.0561
##
## sigma: 10308.58
##
## AIC AICc BIC
## 4462.955 4463.086 4472.665
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249 0.01282239
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995 98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995 98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995 98816.41 83958.37 113674.4 76092.99 121539.87.1.B
autoplot(ses.pigs)
#Find SE/SD for custom prediction interval
sd.pigs <- sd(ses.pigs$residuals)
pigs.upper <- ses.pigs$mean[1] + sd.pigs*1.96
pigs.lower <- ses.pigs$mean[1] - sd.pigs*1.96
pigs.int <- c(pigs.upper,pigs.lower)
#Extract R produced interval
pigs.og <- c(ses.pigs$upper[1,"95%"], ses.pigs$lower[1,"95%"])
#Compare the intervals
both.pigs <- rbind(pigs.og,pigs.int)
kable(both.pigs,col.names = c("Upper","Lower"))| Upper | Lower | |
|---|---|---|
| pigs.og | 119020.8 | 78611.97 |
| pigs.int | 118952.8 | 78679.97 |
7.2
custom.exp.fn <- function(y, alpha, level){
y_hat <- level
alpha_calc <- (1-alpha)
for(i in 1:length(y)){
y_hat <- alpha*y[i] + alpha_calc*y_hat
}
names(y_hat) <- "Custom Fn Prediction:"
return(y_hat)
}
custom.exp.fn(pigs,ses.pigs$model$par[1],ses.pigs$model$par[2])## Custom Fn Prediction:
## 98816.41paste("SES Prediction:",ses.pigs$mean[1])## [1] "SES Prediction: 98816.4061115907"7.3
custom.exp.fn1 <- function(pars= c(alpha, level), y){
y_hat <- pars[2]
alpha <- pars[1]
alpha_calc <- (1-alpha)
SSE <- 0
for(i in 1:length(y)){
error <- y[i]-y_hat
SSE <- SSE + error^2
y_hat <- alpha*y[i] + alpha_calc*y_hat
}
return(SSE)
}
optim.pigs <- optim(par = c(0.5,pigs[1]), y = pigs, fn = custom.exp.fn1)
paste("alpha:",optim.pigs$par[1],"level:",optim.pigs$par[2])## [1] "alpha: 0.299008094014243 level: 76379.2653476235"7.4
custom.exp.fn.final <- function(y){
step.1 <-function(pars= c(alpha, level), y){
y_hat <- pars[2]
alpha <- pars[1]
alpha_calc <- (1-alpha)
SSE <- 0
for(i in 1:length(y)){
error <- y[i]-y_hat
SSE <- SSE + error^2
y_hat <- alpha*y[i] + alpha_calc*y_hat
}
return(SSE)
}
step.2 <- optim(par = c(0.5,y[1]), y = y, fn = step.1)
y_hat <- step.2$par[2]
alpha <- step.2$par[1]
alpha_calc <- (1-alpha)
for(i in 1:length(y)){
y_hat <- alpha*y[i] + alpha_calc*y_hat
}
print(y_hat)
names(y_hat) <- "Prediction:"
return(list(
prediction_val = y_hat,
alpha = alpha,
level = step.2$par[2]
))
}
pigs.run <- custom.exp.fn.final(pigs)## [1] 98814.55paste("alpha:",pigs.run$alpha,"level:",pigs.run$level, "prediction:",pigs.run$prediction_val)## [1] "alpha: 0.299008094014243 level: 76379.2653476235 prediction: 98814.5529660019"7.5.A
autoplot(books)
7.5.B
books.h.ses <- ses(books[,'Hardcover'],h=4)
books.p.ses <- ses(books[,'Paperback'],h=4)
autoplot(books.h.ses, PI = FALSE) + autolayer(books.p.ses, PI = FALSE) 
7.5.C
h.75C <- paste("Hardcover RMSE:",sqrt(mean(books.h.ses$residuals^2)))
p.75C <- paste("Paperback RMSE:",sqrt(mean(books.p.ses$residuals^2)))
h.75C## [1] "Hardcover RMSE: 31.9310149844547"p.75C## [1] "Paperback RMSE: 33.637686782912"7.6.A
books.h.holt <- holt(books[,'Hardcover'],h=4)
books.p.holt <- holt(books[,'Paperback'],h=4)
autoplot(books.h.holt, PI = FALSE) + autolayer(books.p.holt, PI = FALSE)
7.6.B
h.76B <- paste("Hardcover RMSE:",sqrt(mean(books.h.holt$residuals^2)))
p.76B <- paste("Paperback RMSE:",sqrt(mean(books.p.holt$residuals^2)))
h.76B## [1] "Hardcover RMSE: 27.1935779818511"p.76B## [1] "Paperback RMSE: 31.1369230162347"7.6.C
paste(h.75C,"v.",h.76B)## [1] "Hardcover RMSE: 31.9310149844547 v. Hardcover RMSE: 27.1935779818511"paste(p.75C,"v.",p.76B)## [1] "Paperback RMSE: 33.637686782912 v. Paperback RMSE: 31.1369230162347"7.6.D
#Build holt Interval
sd.hb <- sd(books.h.holt$residuals)
hb.upper <- books.h.holt$mean[1] + sd.hb*1.96
hb.lower <- books.h.holt$mean[1] - sd.hb*1.96
hb.holt <- c(hb.upper,hb.lower)
#Auto produced holt interval
hb.holt.og <- c(books.h.holt$upper[1,"95%"], books.h.holt$lower[1,"95%"])
#Build SES Interval
sd.hb.ses <- sd(books.h.ses$residuals)
hb.ses.upper <- books.h.ses$mean[1] + sd.hb.ses*1.96
hb.ses.lower <- books.h.ses$mean[1] - sd.hb.ses*1.96
hb.ses <- c(hb.ses.upper,hb.ses.lower)
#Auto produced SES interval
hb.ses.og <- c(books.h.ses$upper[1,"95%"], books.h.ses$lower[1,"95%"])
#Combine
both.hb <- rbind(hb.holt.og,hb.holt,hb.ses.og,hb.ses)
kable(both.hb,col.names = c("Upper","Lower"))| Upper | Lower | |
|---|---|---|
| hb.holt.og | 307.4256 | 192.9222 |
| hb.holt | 304.3838 | 195.9640 |
| hb.ses.og | 304.3403 | 174.7799 |
| hb.ses | 300.5354 | 178.5848 |
#Build holt Interval
sd.pb <- sd(books.p.holt$residuals)
pb.upper <- books.p.holt$mean[1] + sd.pb*1.96
pb.lower <- books.p.holt$mean[1] - sd.pb*1.96
pb.holt <- c(pb.upper,pb.lower)
#Auto produced holt interval
pb.holt.og <- c(books.p.holt$upper[1,"95%"], books.p.holt$lower[1,"95%"])
#Build SES Interval
sd.pb.ses <- sd(books.p.ses$residuals)
pb.ses.upper <- books.p.ses$mean[1] + sd.pb.ses*1.96
pb.ses.lower <- books.p.ses$mean[1] - sd.pb.ses*1.96
pb.ses <- c(pb.ses.upper,pb.ses.lower)
#Auto produced SES interval
pb.ses.og <- c(books.p.ses$upper[1,"95%"], books.p.ses$lower[1,"95%"])
#Combine
both.pb <- rbind(pb.holt.og,pb.holt,pb.ses.og,pb.ses)
kable(both.pb,col.names = c("Upper","Lower"))| Upper | Lower | |
|---|---|---|
| pb.holt.og | 275.0205 | 143.9130 |
| pb.holt | 271.0945 | 147.8390 |
| pb.ses.og | 275.3523 | 138.8670 |
| pb.ses | 272.6230 | 141.5964 |
7.7
#Base
eggs.holt.1 <- holt(eggs, h = 100)
paste("Holt 1:",sqrt(mean(eggs.holt.1$residuals^2)))## [1] "Holt 1: 26.5821881569903"#Damped
eggs.holt.2 <- holt(eggs, h = 100, damped = TRUE)
paste("Holt 2:",sqrt(mean(eggs.holt.2$residuals^2)))## [1] "Holt 2: 26.5401852798325"#Box-Cox
eggs.holt.3 <- holt(eggs, h = 100, lambda = BoxCox.lambda(eggs))
paste("Holt 3:",sqrt(mean(eggs.holt.3$residuals^2)))## [1] "Holt 3: 1.03221710764543"#Plot
autoplot(eggs) +
autolayer(eggs.holt.1, series="Holt", PI=FALSE) +
autolayer(eggs.holt.2, series="Holt Damped", PI=FALSE) +
autolayer(eggs.holt.3, series="Holt Box-Cox", PI=FALSE) 
7.8.A
retail <- readxl::read_excel("retail.xlsx", skip = 1)
retail.ts <- ts(retail[,"A3349337W"],frequency=12, start=c(1982,4))
retail.ts.dc <- seas(retail.ts, x11="")
autoplot(retail.ts.dc)
7.8.A Response:
Multiplicative seasonality is necessary for this time series because the level of seasonality changes over time for this data.
7.8.B
retail.hw <- hw(retail.ts, "multiplicative", h=24)
retail.hw.damp <- hw(retail.ts, "multiplicative", h=24, damped = TRUE)
autoplot(retail.ts) +
autolayer(retail.hw, series="Holt-Winters", PI=FALSE) +
autolayer(retail.hw.damp, series="Holt-Winters Damped", PI=FALSE)
7.8.C
hw.acc <- accuracy(retail.hw)
hw.damp.acc <- accuracy(retail.hw.damp)
hw.acc## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.5571551 13.17456 9.918904 0.3439157 5.973236 0.4821535 0.1309425hw.damp.acc## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.4168871 13.15306 9.799695 0.2196695 5.827413 0.4763588 0.1060087.8.D
checkresiduals(retail.hw.damp)
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 150.7, df = 7, p-value < 2.2e-16
##
## Model df: 17. Total lags used: 247.8.E
retail.window.train<-window(retail.ts, end=c(2010,12))
retail.window.test<-window(retail.ts, start=c(2011, 1))
retail.window.naive<-snaive(retail.window.train, h=36)
retail.window.hw.damp <- hw(retail.window.train, seasonal="multiplicative", h=36, damped=TRUE)
snaive.acc <- accuracy(retail.window.naive, retail.window.test)
hw.damp.acc1 <- accuracy(retail.window.hw.damp, retail.window.test)
snaive.acc## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.460661 26.30758 21.23363 4.655690 12.762886 1.0000000 0.8070166
## Test set 14.094444 20.91121 17.30556 3.802915 4.911945 0.8150068 0.5265917
## Theil's U
## Training set NA
## Test set 0.6200865hw.damp.acc1## ME RMSE MAE MPE MAPE MASE
## Training set 0.1340865 13.14634 9.755892 0.05575759 6.322517 0.4594546
## Test set 15.1173579 23.10869 18.810903 3.99168849 5.241861 0.8859013
## ACF1 Theil's U
## Training set 0.09808033 NA
## Test set 0.63781919 0.6747655#Plot
autoplot(retail.ts) +
autolayer(retail.window.naive, series="Naive", PI=FALSE) +
autolayer(retail.window.hw.damp, series="Holt-Winters Damped", PI=FALSE)
7.9
retail.window.train.stl <- stlm(retail.window.train, lambda=BoxCox.lambda(retail.window.train),s.window=13, robust=TRUE, method="ets")
retail.window.for.stl <- forecast(retail.window.train.stl, lambda=BoxCox.lambda(retail.window.train), h = 36)## Warning in InvBoxCox(fcast$mean, lambda, biasadj, fcast): biasadj information
## not found, defaulting to FALSE.accuracy(retail.window.for.stl,retail.window.test)## ME RMSE MAE MPE MAPE MASE
## Training set 1.160974 12.03861 8.398332 0.8282614 5.354964 0.3955203
## Test set 17.744477 27.87913 21.913666 4.6529147 6.087432 1.0320262
## ACF1 Theil's U
## Training set 0.01707693 NA
## Test set 0.69510092 0.8203021#Plot
autoplot(retail.ts) +
autolayer(retail.window.naive, series="Naive", PI=FALSE) +
autolayer(retail.window.hw.damp, series="Holt-Winters Damped", PI=FALSE) +
autolayer(retail.window.for.stl, series="STL", PI=FALSE)
7.10.A
autoplot(ukcars)
7.10.B
ukcars.stl <- stl(ukcars,t.window=13, s.window="periodic", robust=TRUE)
autoplot(ukcars.stl)
7.10.C
ukcars.stlf <- stlf(ukcars,etsmodel="AAN", damped=TRUE, h = 8)
autoplot(ukcars.stlf)
7.10.D
ukcars.stlf.l <- stlf(ukcars,etsmodel="AAN", damped=FALSE, h = 8)
autoplot(ukcars.stlf.l)
7.10.E
ukcars.ets <- ets(ukcars)
ukcars.ets.for <- forecast(ukcars.ets,h=8)
autoplot(ukcars.ets.for)
7.10.F
accuracy(ukcars.stlf)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.551267 23.32113 18.48987 0.04121971 6.042764 0.602576 0.02262668accuracy(ukcars.stlf.l)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.3412727 23.295 18.1605 -0.5970778 5.98018 0.5918418 0.02103582accuracy(ukcars.ets.for)## ME RMSE MAE MPE MAPE MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
## ACF1
## Training set 0.025733347.10.G
autoplot(ukcars) +
autolayer(ukcars.stlf, series="stlf", PI=FALSE) +
autolayer(ukcars.stlf.l, series="stlf linear", PI=FALSE) +
autolayer(ukcars.ets.for, series="ets", PI=FALSE)
7.10.H
checkresiduals(ukcars.stlf.l)## Warning in checkresiduals(ukcars.stlf.l): The fitted degrees of freedom is based
## on the model used for the seasonally adjusted data.
##
## Ljung-Box test
##
## data: Residuals from STL + ETS(A,A,N)
## Q* = 22.061, df = 4, p-value = 0.0001949
##
## Model df: 4. Total lags used: 87.11.A
autoplot(visitors)
7.11.B
visitors.train<-window(visitors, end=c(2003,12))
visitors.test<-window(visitors, start=c(2004, 1))
visitors.hw <- hw(visitors.train,seasonal="multiplicative", h=24)
autoplot(visitors.hw)
7.11.C
Multiplicative is needed, as the seasonality component is changing over time.
7.11.D
visitors.ets <- ets(visitors.train)
visitors.ets.forecast<- forecast(visitors.ets, h = 24)
visitors.ets.add <- ets(visitors.train, additive.only = TRUE, lambda = BoxCox.lambda(visitors.train))
visitors.ets.add.forecast <- forecast(visitors.ets.add, h = 24)
visitors.snaive<-snaive(visitors.train, h = 24)
visitors.snaive<-snaive(visitors.train, h = 24)
visitors.stlf<-stlf(visitors.train, t.window = 13, s.window = "periodic", h = 24, robust = TRUE, method = "ets", lambda = BoxCox.lambda(visitors.train))
autoplot(visitors) +
autolayer(visitors.ets.forecast, series="ets", PI=FALSE) +
autolayer(visitors.ets.add.forecast, series="ets additive", PI=FALSE) +
autolayer(visitors.snaive, series="snaive", PI=FALSE)+
autolayer(visitors.stlf, series="stlf", PI=FALSE)
7.11.E
print("ets")## [1] "ets"accuracy(visitors.ets.forecast,visitors.test)## ME RMSE MAE MPE MAPE MASE
## Training set 1.232439 14.95912 10.80756 0.2153496 4.044462 0.4138209
## Test set 12.617119 22.90797 16.63487 2.5783148 3.521617 0.6369485
## ACF1 Theil's U
## Training set -0.002560158 NA
## Test set 0.518168412 0.3411748print("ets additive")## [1] "ets additive"accuracy(visitors.ets.add.forecast,visitors.test)## ME RMSE MAE MPE MAPE MASE
## Training set 1.168056 15.51547 11.22989 0.11241764 4.046645 0.4299918
## Test set 1.581019 20.22801 16.00030 -0.04770367 3.536397 0.6126507
## ACF1 Theil's U
## Training set -0.02809955 NA
## Test set 0.67007296 0.3064412print("snaive")## [1] "snaive"accuracy(visitors.snaive,visitors.test)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 16.59292 31.24792 26.11651 6.836282 10.18912 1.000000 0.6514168
## Test set 50.58125 59.02680 50.58125 11.730236 11.73024 1.936754 0.6653705
## Theil's U
## Training set NA
## Test set 0.9308711print("stlf")## [1] "stlf"accuracy(visitors.stlf,visitors.test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.02532068 15.39434 11.08552 -0.09928427 4.107846 0.4244642
## Test set -18.53595827 24.35689 20.18924 -4.57352057 4.909926 0.7730451
## ACF1 Theil's U
## Training set -0.02566624 NA
## Test set 0.12433192 0.3995668checkresiduals(visitors.ets.forecast)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,Ad,M)
## Q* = 22.713, df = 7, p-value = 0.001912
##
## Model df: 17. Total lags used: 247.11.F
fets <- function(y, h) {
forecast(ets(y), h = h)
}
fets_add <- function(y, h) {
forecast(ets(y,additive.only = TRUE, lambda = BoxCox.lambda(y)), h = h)
}
fstl <- function(y, h) {
forecast(stlf(y,t.window = 13, s.window = "periodic",robust = TRUE, method = "ets", lambda = BoxCox.lambda(y), h = h))
}
visitors.ets.csv<-tsCV(visitors,fets,h=1)
visitors.ets.add.cv<-tsCV(visitors,fets_add,h=1)
visitors.snaive.cv<-tsCV(visitors,snaive,h=1)
visitors.stlf.cv<-tsCV(visitors,fstl,h=1)
print("ets")## [1] "ets"sqrt(mean(visitors.ets.csv^2, na.rm = TRUE))## [1] 18.52985print("ets additive")## [1] "ets additive"sqrt(mean(visitors.ets.add.cv^2, na.rm = TRUE))## [1] 18.86439print("snaive")## [1] "snaive"sqrt(mean(visitors.snaive.cv^2, na.rm = TRUE))## [1] 32.78675print("stlf")## [1] "stlf"sqrt(mean(visitors.stlf.cv^2, na.rm = TRUE))## [1] 18.867047.12.A
qcement.ets.csv<-tsCV(qcement,fets,h=1)
qcement.snaive.cv<-tsCV(qcement,snaive,h=1)
qcement.ets.csv## Qtr1 Qtr2 Qtr3 Qtr4
## 1956 6.700000e-02 9.265000e-02 -2.442844e-02 -7.198045e-02
## 1957 7.500000e-02 9.175809e-03 -1.891133e-02 -3.235806e-02
## 1958 6.058900e-02 4.919003e-02 -1.878542e-03 -6.423717e-02
## 1959 7.523263e-02 1.086467e-02 1.571592e-02 -3.976007e-03
## 1960 1.082088e-02 3.148850e-02 6.303604e-03 2.540389e-02
## 1961 -1.569300e-02 -3.216928e-02 -6.815775e-02 -4.571899e-03
## 1962 4.543476e-02 1.206082e-02 4.557214e-03 -2.279211e-02
## 1963 1.909600e-02 4.028200e-02 3.200105e-02 6.035647e-02
## 1964 3.512984e-02 1.833127e-02 6.006993e-02 2.987929e-03
## 1965 -3.235374e-02 -3.992800e-02 1.989288e-02 -4.148989e-02
## 1966 3.109235e-03 -2.380817e-02 -2.543793e-02 -1.118874e-02
## 1967 5.741866e-02 -3.286375e-03 -3.352551e-02 3.687221e-03
## 1968 -4.154215e-02 4.276720e-02 -7.073400e-03 1.646154e-02
## 1969 -3.465909e-02 1.235001e-01 -4.140226e-03 8.169693e-03
## 1970 -4.347429e-02 3.218184e-02 1.217650e-03 -2.729623e-02
## 1971 1.720489e-02 -5.150032e-02 7.712976e-02 -6.040211e-02
## 1972 5.440756e-02 -8.166172e-03 5.787502e-03 -1.984992e-02
## 1973 1.653694e-03 -7.416339e-03 1.326106e-01 -7.849817e-02
## 1974 -2.279525e-02 -9.108580e-02 2.946654e-02 -9.323637e-02
## 1975 -9.185140e-03 -2.602798e-02 -5.019609e-02 -5.336273e-02
## 1976 2.699699e-02 6.380161e-03 -2.707303e-02 -6.387311e-02
## 1977 1.769847e-03 -2.056736e-03 1.604306e-02 -2.865398e-02
## 1978 -4.917760e-02 -1.392407e-03 -4.856133e-03 8.428459e-02
## 1979 -4.013761e-02 2.851826e-02 1.852664e-02 2.257385e-02
## 1980 1.789710e-04 9.042296e-03 1.472243e-02 8.902635e-03
## 1981 1.245945e-01 -4.748533e-02 1.173887e-01 -9.193956e-02
## 1982 2.601697e-02 -1.391883e-01 -1.751930e-01 -1.614772e-01
## 1983 -1.306019e-01 1.170459e-01 8.953425e-02 1.561802e-02
## 1984 6.627619e-02 5.409295e-02 6.198187e-02 -6.074919e-02
## 1985 6.390279e-03 1.373993e-01 -3.975074e-02 2.067906e-02
## 1986 -1.359085e-02 -3.041629e-02 -4.034562e-04 -4.223154e-02
## 1987 -7.239964e-02 -1.776867e-02 1.388662e-01 2.085564e-02
## 1988 3.330111e-02 9.755089e-02 4.568548e-02 8.369269e-02
## 1989 -1.171832e-01 1.610487e-01 -3.389748e-02 -8.002673e-02
## 1990 -1.650649e-01 -1.043661e-01 1.369842e-02 -8.155152e-02
## 1991 -9.767305e-02 -1.235748e-01 5.166749e-02 2.749059e-02
## 1992 5.985122e-02 -4.034956e-02 2.032241e-02 1.124560e-01
## 1993 3.802435e-02 -9.599193e-02 2.096743e-01 -9.954993e-02
## 1994 8.154217e-02 1.495824e-01 -1.345394e-01 4.275673e-02
## 1995 -1.338330e-01 -8.763383e-02 -2.128407e-01 1.440840e-01
## 1996 -9.895576e-02 2.800160e-02 -2.360591e-06 1.490874e-02
## 1997 9.813113e-02 -1.325704e-02 2.168898e-02 1.199172e-01
## 1998 -3.833019e-02 7.107875e-03 1.135933e-01 -5.235061e-02
## 1999 -6.147776e-03 -3.352216e-02 3.090990e-02 7.209592e-02
## 2000 5.208887e-02 -2.505475e-01 -2.846975e-01 2.668633e-02
## 2001 -2.076136e-02 -9.684253e-02 1.501796e-01 1.445685e-01
## 2002 9.107351e-02 3.841723e-03 -2.713982e-02 3.986384e-02
## 2003 -1.238251e-01 2.657085e-01 -5.784606e-02 9.731096e-02
## 2004 5.914216e-03 5.353615e-02 -1.133499e-01 4.269288e-02
## 2005 2.730024e-01 -1.526262e-01 -9.012514e-02 1.224639e-02
## 2006 -3.078828e-03 7.263933e-02 8.514199e-02 -2.093492e-02
## 2007 -4.586818e-02 8.309848e-02 7.198085e-02 -7.204070e-02
## 2008 8.016897e-02 -1.135926e-02 -2.132501e-01 -1.868232e-01
## 2009 -8.451826e-02 7.033035e-02 9.627761e-03 -1.100570e-01
## 2010 2.460655e-01 6.684919e-02 -1.522665e-01 -7.235882e-03
## 2011 -3.907263e-02 1.354478e-01 -4.901492e-02 -5.939948e-02
## 2012 -1.079507e-01 1.232096e-01 1.161515e-01 -1.378282e-01
## 2013 1.764193e-01 6.906144e-02 -1.303157e-02 -3.665988e-02
## 2014 NAqcement.snaive.cv## Qtr1 Qtr2 Qtr3 Qtr4
## 1956 0.067 0.096 0.105 0.064
## 1957 0.072 0.042 0.012 0.025
## 1958 0.016 0.043 0.055 0.019
## 1959 0.053 0.044 0.044 0.048
## 1960 0.025 0.063 0.047 0.067
## 1961 0.039 -0.011 -0.036 -0.051
## 1962 0.020 0.041 0.065 0.037
## 1963 0.017 0.052 0.081 0.123
## 1964 0.130 0.114 0.137 0.105
## 1965 0.070 0.020 -0.008 -0.053
## 1966 -0.013 -0.003 -0.045 -0.013
## 1967 0.037 0.059 0.049 0.056
## 1968 -0.025 0.022 0.046 0.056
## 1969 0.059 0.143 0.119 0.101
## 1970 0.102 0.039 0.052 0.009
## 1971 0.075 -0.030 0.065 0.012
## 1972 0.073 0.104 0.028 0.078
## 1973 0.023 0.039 0.171 0.036
## 1974 0.047 -0.026 -0.099 -0.063
## 1975 -0.055 0.002 -0.072 -0.012
## 1976 0.002 0.023 0.009 -0.008
## 1977 -0.015 -0.020 0.027 0.021
## 1978 -0.023 -0.009 -0.019 0.083
## 1979 0.041 0.065 0.061 0.027
## 1980 0.084 0.063 0.080 0.050
## 1981 0.187 0.093 0.202 0.055
## 1982 0.002 -0.058 -0.284 -0.251
## 1983 -0.417 -0.196 -0.044 0.080
## 1984 0.215 0.177 0.182 0.085
## 1985 0.083 0.175 0.060 0.124
## 1986 0.090 -0.063 0.006 -0.043
## 1987 -0.088 -0.067 0.054 0.086
## 1988 0.167 0.269 0.176 0.203
## 1989 0.094 0.202 0.103 -0.056
## 1990 -0.074 -0.314 -0.226 -0.222
## 1991 -0.204 -0.214 -0.171 -0.054
## 1992 0.060 0.095 0.071 0.161
## 1993 0.167 0.109 0.295 0.018
## 1994 0.087 0.314 -0.030 0.158
## 1995 -0.052 -0.229 -0.288 -0.100
## 1996 -0.110 -0.027 0.154 -0.015
## 1997 0.192 0.120 0.131 0.208
## 1998 0.076 0.130 0.237 0.018
## 1999 0.083 0.049 -0.040 0.098
## 2000 0.126 -0.107 -0.375 -0.281
## 2001 -0.353 -0.219 0.184 0.175
## 2002 0.275 0.351 0.142 0.102
## 2003 -0.100 0.197 0.112 0.132
## 2004 0.288 0.080 0.067 0.017
## 2005 0.301 0.033 0.108 0.047
## 2006 -0.203 0.087 0.186 0.113
## 2007 0.084 0.109 0.111 0.043
## 2008 0.196 0.076 -0.189 -0.220
## 2009 -0.398 -0.287 -0.100 -0.059
## 2010 0.241 0.169 0.023 0.151
## 2011 -0.128 0.005 0.094 0.012
## 2012 -0.050 -0.048 0.113 -0.018
## 2013 0.305 0.186 0.062 0.180
## 2014 NA7.12.B
print("ets")## [1] "ets"mean(qcement.ets.csv^2, na.rm = TRUE)## [1] 0.006959511print("snaive")## [1] "snaive"mean(qcement.snaive.cv^2, na.rm = TRUE)## [1] 0.017792437.13
runModels <- function(trainData, testData,h, windowT ){
model.ets <- forecast(ets(trainData), h = h)
model.snaive <- snaive(trainData, h = h)
model.stlf <- forecast(stlf(trainData,t.window = windowT, s.window = "periodic",robust = TRUE, lambda = BoxCox.lambda(trainData), h = h))
returnList <- list()
returnList[["ets"]] <- accuracy(model.ets,testData)
returnList[["snaive"]] <- accuracy(model.snaive,testData)
returnList[["stlf"]] <- accuracy(model.stlf,testData)
return(returnList)
}
getDataSets <-function(y){
y.freq <- frequency(y)
year_3 <- y.freq * 3
y.length <- length(y)
time(y)[y.length]
y.test <- window(y, start = time(y)[(y.length-year_3)+1])
y.train <- window(y, end = time(y)[(y.length-year_3)])
returnList <- list()
returnList[["train"]] <-y.train
returnList[["test"]] <-y.test
returnList[["frequency"]] <-y.freq
return(returnList)
}
runFull <- function(y){
all.sets <- getDataSets(y)
all.results <- runModels(all.sets[["train"]],all.sets[["test"]], all.sets[["frequency"]] * 3,all.sets[["frequency"]]+ 1)
return(all.results)
}
print("ausbeer")## [1] "ausbeer"print("bricksq")## [1] "bricksq"runFull(bricksq)## $ets
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.228708 21.96146 15.81586 0.2606171 3.912882 0.4299638 0.2038074
## Test set 37.682916 42.36696 37.68292 8.4157549 8.415755 1.0244329 0.6190546
## Theil's U
## Training set NA
## Test set 1.116608
##
## $snaive
## ME RMSE MAE MPE MAPE MASE
## Training set 6.438849 50.25482 36.78417 1.430237 8.999949 1.0000000
## Test set 15.333333 37.15284 33.66667 3.231487 7.549326 0.9152487
## ACF1 Theil's U
## Training set 0.81169827 NA
## Test set -0.09314867 0.9538702
##
## $stlf
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.391509 21.85179 15.23445 0.3066952 3.743250 0.4141577 0.1436315
## Test set 37.432254 42.22598 37.43225 8.3389919 8.338992 1.0176185 0.6226050
## Theil's U
## Training set NA
## Test set 1.111572print("dole")## [1] "dole"runFull(dole)## $ets
## ME RMSE MAE MPE MAPE MASE
## Training set 98.00253 16130.58 9427.497 0.5518769 6.20867 0.2995409
## Test set 221647.42595 279668.65 221647.426 32.5447637 32.54476 7.0424271
## ACF1 Theil's U
## Training set 0.5139536 NA
## Test set 0.9394267 11.29943
##
## $snaive
## ME RMSE MAE MPE MAPE MASE
## Training set 12606.58 56384.06 31473.16 3.350334 27.73651 1.000000
## Test set 160370.33 240830.92 186201.00 20.496197 27.38678 5.916184
## ACF1 Theil's U
## Training set 0.9781058 NA
## Test set 0.9208325 9.479785
##
## $stlf
## ME RMSE MAE MPE MAPE MASE
## Training set 283.3964 7930.906 4944.778 0.2303326 5.835319 0.1571109
## Test set 210956.3986 272047.968 211808.711 30.4291034 30.650066 6.7298206
## ACF1 Theil's U
## Training set 0.1264403 NA
## Test set 0.9390606 10.87011print("a10")## [1] "a10"runFull(a10)## $ets
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04378049 0.488989 0.3542285 0.2011915 3.993267 0.3611299
## Test set 1.62559331 2.562046 2.0101563 7.0011207 9.410027 2.0493197
## ACF1 Theil's U
## Training set -0.05238173 NA
## Test set 0.23062972 0.6929693
##
## $snaive
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.9577207 1.177528 0.9808895 10.86283 11.15767 1.000000 0.3779746
## Test set 4.3181585 5.180738 4.3306974 19.80542 19.89352 4.415071 0.6384822
## Theil's U
## Training set NA
## Test set 1.383765
##
## $stlf
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002177756 0.4553175 0.3313982 -0.1283507 3.847741 0.3378548
## Test set 1.176903363 2.1546770 1.8131163 4.6090383 8.631147 1.8484409
## ACF1 Theil's U
## Training set -0.1068058 NA
## Test set 0.1111434 0.5855374print("h02")## [1] "h02"runFull(h02)## $ets
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001532209 0.04653434 0.03463451 0.0008075938 4.575471 0.5850192
## Test set 0.023667588 0.07667744 0.06442193 1.7152319315 7.030603 1.0881653
## ACF1 Theil's U
## Training set 0.07982687 NA
## Test set -0.12074883 0.450176
##
## $snaive
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03880677 0.07122149 0.05920234 5.220203 8.156509 1.00000
## Test set -0.01479486 0.08551752 0.07161055 -1.308298 7.884556 1.20959
## ACF1 Theil's U
## Training set 0.40465289 NA
## Test set 0.02264221 0.5009677
##
## $stlf
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001690205 0.04522680 0.03309428 -0.1245389 4.428712 0.5590028
## Test set 0.029415081 0.07651069 0.06449846 2.1176176 7.063467 1.0894578
## ACF1 Theil's U
## Training set -0.0591871 NA
## Test set -0.1325979 0.4483587print("usmelec")## [1] "usmelec"runFull(usmelec)## $ets
## ME RMSE MAE MPE MAPE MASE
## Training set -0.07834851 7.447167 5.601963 -0.1168709 2.163495 0.6225637
## Test set -10.20080652 15.291359 13.123441 -3.1908161 3.926338 1.4584491
## ACF1 Theil's U
## Training set 0.2719249 NA
## Test set 0.4679496 0.4859495
##
## $snaive
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4.921564 11.58553 8.998217 2.000667 3.511648 1.000000 0.4841860
## Test set 5.475972 17.20037 13.494750 1.391437 3.767514 1.499714 0.2692544
## Theil's U
## Training set NA
## Test set 0.4765145
##
## $stlf
## ME RMSE MAE MPE MAPE MASE
## Training set -0.07431421 7.242232 5.336692 -0.03633071 2.076773 0.5930833
## Test set -24.40009972 27.993269 24.620052 -7.36385518 7.424595 2.7361034
## ACF1 Theil's U
## Training set 0.02701321 NA
## Test set 0.60890966 0.89613797.14.A
print("bicoal")## [1] "bicoal"accuracy(ets(bicoal))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.07050762 59.53 46.46042 -0.9663802 10.1218 0.9982543 0.03586774print("chicken")## [1] "chicken"accuracy(ets(chicken))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.116415 13.38863 9.53233 -4.930494 13.11671 0.9922021 0.02763414print("dole")## [1] "dole"accuracy(ets(dole))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 569.2395 16559.97 10000.7 0.5960691 5.910033 0.2424801 0.5556989print("usdeaths")## [1] "usdeaths"accuracy(ets(usdeaths))## ME RMSE MAE MPE MAPE MASE
## Training set -1.262752 264.2916 205.3357 -0.06186187 2.360075 0.4699475
## ACF1
## Training set -0.008409405print("lynx")## [1] "lynx"accuracy(ets(lynx))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 8.975647 1198.452 842.0649 -52.12968 101.3686 1.013488 0.3677583print("ibmclose")## [1] "ibmclose"accuracy(ets(ibmclose))## ME RMSE MAE MPE MAPE MASE
## Training set -0.2789829 7.243976 5.195358 -0.08443177 1.19398 0.9973353
## ACF1
## Training set 0.08562507print("eggs")## [1] "eggs"accuracy(ets(eggs))## ME RMSE MAE MPE MAPE MASE
## Training set -2.800547 26.58125 19.24163 -2.956833 10.02174 0.9491609
## ACF1
## Training set 0.027595517.14.B
autoplot(ibmclose) +autolayer(forecast(ets(ibmclose), h =24))
7.15
forecast(ets(usdeaths, model="MAM"), h =1)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1979 8233.107 7883.699 8582.515 7698.734 8767.481hw(usdeaths, seasonal="multiplicative", h=1)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1979 8217.64 7845.352 8589.928 7648.274 8787.0057.16
model.ets <- ets(usdeaths, model="ANN")
summary(model.ets)## ETS(A,N,N)
##
## Call:
## ets(y = usdeaths, model = "ANN")
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 9009.107
##
## sigma: 735.891
##
## AIC AICc BIC
## 1262.447 1262.800 1269.277
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3.206309 725.5983 633.6029 -0.3128852 7.28904 1.450113 0.02235815model.ets.for <- forecast(model.ets, h = 1)
alpha <- model.ets$par[1]
mysigma <- sqrt(model.ets$sigma2)
variation <- ((mysigma)*(1+(alpha^2)*(1)))
model.ets.for$mean[1]+variation## alpha
## 10711.57model.ets.for$upper[1,"95%"]## 95%
## 10682.267.17
Forecast = lt + (1-alpha) * lt-1 Variance = (sigma^2 * (1+(alpha^2)* (h - 1))) Prediction Interval Upper Bound = lt + (1-alpha) * lt-1 + (sigma^2 * (1+(alpha^2)* (h - 1)))
Chapter 8 (8) Exercises
8.1.A
All values fall within the 95% interval (indicated by the blue dotted lines), so it’s safe to say that the variation we see is just white noise.
8.1.B
Generally, this is the result of the central limit theorem, i.e. the fact that adding more observations narrows the width of confidence intervals, as more observations make us more sure of the “true” mean.
8.2
autoplot(ibmclose)
acf(ibmclose)
pacf(ibmclose)
8.3.A
boxcox.usnetelec <- BoxCox(usnetelec, BoxCox.lambda(usnetelec))
summary(ur.kpss(boxcox.usnetelec))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 1.4583
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(BoxCox(usnetelec,BoxCox.lambda(usnetelec)))
ggtsdisplay(diff(BoxCox(usnetelec, BoxCox.lambda(usnetelec)), differences = holdDiff))
8.3.B
boxcox.usgdp <- BoxCox(usgdp, BoxCox.lambda(usgdp))
summary(ur.kpss(boxcox.usgdp))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 4 lags.
##
## Value of test-statistic is: 4.8114
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(BoxCox(usgdp,BoxCox.lambda(usgdp)))
acf(diff(BoxCox(usgdp, BoxCox.lambda(usgdp)), differences = holdDiff))
ggtsdisplay(diff(BoxCox(usgdp, BoxCox.lambda(usgdp)), differences = 2))
8.3.C
boxcox.mcopper <- BoxCox(mcopper, BoxCox.lambda(mcopper))
summary(ur.kpss(boxcox.mcopper))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 6 lags.
##
## Value of test-statistic is: 6.2659
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(boxcox.mcopper)
ggtsdisplay(diff(boxcox.mcopper, differences = 2))
8.3.D
boxcox.enplanements <- BoxCox(enplanements, BoxCox.lambda(enplanements))
summary(ur.kpss(boxcox.enplanements))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 5 lags.
##
## Value of test-statistic is: 4.3785
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(boxcox.enplanements)
ggtsdisplay(diff(boxcox.enplanements, differences = holdDiff))
8.3.E
boxcox.visitors <- BoxCox(visitors, BoxCox.lambda(visitors))
summary(ur.kpss(boxcox.visitors))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 4 lags.
##
## Value of test-statistic is: 4.5233
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(boxcox.visitors)
ggtsdisplay(diff(boxcox.visitors, differences = holdDiff))
8.5
retaildata <- readxl::read_excel("retail.xlsx", skip = 1)
myts <- ts(retaildata[,"A3349337W"],frequency=12, start=c(1982,4))
ggtsdisplay(myts)
boxcox.myts <- BoxCox(myts, BoxCox.lambda(myts))
summary(ur.kpss(boxcox.myts))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 5 lags.
##
## Value of test-statistic is: 5.8558
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(boxcox.myts)
ggtsdisplay(diff(boxcox.myts, differences = holdDiff))
ggtsdisplay(boxcox.myts %>% diff(holdDiff) %>% diff(12))
8.6.A
ar.model <- function(phi){
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <-phi*y[i-1] + e[i]
}
return(y)
}
ar.model.6 <- ar.model(.6)
autoplot(ar.model.6)
8.6.B
ar.model.5 <- ar.model(.5)
ar.model.2 <- ar.model(.2)
ar.model.8 <- ar.model(.8)
autoplot(ar.model.6, series = ".6") + autolayer(ar.model.5,series =".5") + autolayer(ar.model.2,series =".2")+ autolayer(ar.model.8,series =".8")
8.6.C
ma.model <- function(theta){
y <- ts(numeric(100))
e <- rnorm(100)
e[1] <- 0
for(i in 2:100){
y[i] <-theta*e[i-1] + e[i]
}
return(y)
}
ma.model.6 <- ma.model(.6)8.6.D
ma.model.2 <- ma.model(.2)
ma.model.8 <- ma.model(.8)
autoplot(ma.model.6, series = ".6") + autolayer(ma.model.2,series =".2")+ autolayer(ma.model.8,series =".8")
8.6.E
arma.model <- function(theta, phi){
y <- ts(numeric(100))
e <- rnorm(100)
e[1] <- 0
for(i in 2:100){
y[i] <-phi*y[i-1]+ theta*e[i-1] + e[i]
}
return(y)
}
arma.model.hold <- arma.model(.6,.6)8.6.F
ar.model <- function(phi, phi2){
y <- ts(numeric(100))
e <- rnorm(100)
e[1] <- 0
for(i in 3:100){
y[i] <- phi*y[i-1]+ e[i] + phi2*y[i-2]
}
return(y)
}
ar.model.hold <- ar.model(-.8,.3)8.6.G
autoplot(arma.model.hold, series = "arma") + autolayer(ar.model.hold,series = "ar")
8.7.A
ggtsdisplay(wmurders)
boxcox.wmurders <- BoxCox(wmurders, BoxCox.lambda(wmurders))
summary(ur.kpss(boxcox.wmurders))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.6745
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739holdDiff <- ndiffs(boxcox.wmurders)
ggtsdisplay(diff(boxcox.wmurders, differences = holdDiff))
8.7.B
Yes, this model should include a constant. Without it, the model would follow a quadratic trend, which would throw off it’s accuracy/worsen its fit.
8.7.C
(1-beta)^2 * y_t
8.7.D
arima_wmurds_020 <- Arima(wmurders, order = c(0,2,0))
checkresiduals(arima_wmurds_020)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,0)
## Q* = 42.219, df = 10, p-value = 6.856e-06
##
## Model df: 0. Total lags used: 108.7.E
arima_wmurds_020.for <- forecast(arima_wmurds_020, h = 3)8.7.F
autoplot(arima_wmurds_020.for)
8.7.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.97summary(arima_wmurds_020)## Series: wmurders
## ARIMA(0,2,0)
##
## sigma^2 estimated as 0.1007: log likelihood=-14.38
## AIC=30.76 AICc=30.84 BIC=32.73
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0001657022 0.3115711 0.2248982 -0.02948252 6.457417 1.383019
## ACF1
## Training set -0.68301888.8.A
austa.auto.arima <- auto.arima(austa)
austa.auto.arima## 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.57austa.auto.arima.for <- forecast(austa.auto.arima,h = 10)
autoplot(austa.auto.arima.for)
8.8.B
arima_austa_011_nod <- Arima(austa, order = c(0,1,1), include.drift = FALSE)
arima_austa_011_nod## Series: austa
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## 0.4936
## s.e. 0.1265
##
## sigma^2 estimated as 0.04833: log likelihood=3.73
## AIC=-3.45 AICc=-3.08 BIC=-0.34arima_austa_011_nod.for <- forecast(arima_austa_011_nod, h = 10)
autoplot(arima_austa_011_nod.for)
arima_austa_010 <- Arima(austa, order = c(0,1,0), include.drift = FALSE)
arima_austa_010## Series: austa
## ARIMA(0,1,0)
##
## sigma^2 estimated as 0.06423: log likelihood=-1.62
## AIC=5.24 AICc=5.36 BIC=6.8arima_austa_010.for <- forecast(arima_austa_010, h = 10)
autoplot(arima_austa_010.for)
8.8.C
arima_austa_213 <- Arima(austa, order = c(2,1,3), include.drift = TRUE)
arima_austa_213## Series: austa
## ARIMA(2,1,3) with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 drift
## 1.7648 -0.8513 -1.7684 0.5571 0.2224 0.1725
## s.e. 0.1136 0.1182 0.2433 0.4351 0.2150 0.0051
##
## sigma^2 estimated as 0.0276: log likelihood=13.72
## AIC=-13.44 AICc=-9.29 BIC=-2.55arima_austa_213.for <- forecast(arima_austa_213, h = 10)
autoplot(arima_austa_213.for)
arima_austa_213_const <- Arima(austa-mean(austa), order = c(2,1,3), include.drift = TRUE)
arima_austa_213_const## Series: austa - mean(austa)
## ARIMA(2,1,3) with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 drift
## 1.7648 -0.8513 -1.7684 0.5571 0.2224 0.1725
## s.e. 0.1136 0.1182 0.2433 0.4351 0.2150 0.0051
##
## sigma^2 estimated as 0.0276: log likelihood=13.72
## AIC=-13.44 AICc=-9.29 BIC=-2.55arima_austa_213_const.for <- forecast(arima_austa_213_const, h = 10)
autoplot(arima_austa_213_const.for)
8.8.D
arima_austa_001 <- Arima(austa, order = c(0,0,1))
arima_austa_001## Series: austa
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 mean
## 1.0000 3.5515
## s.e. 0.0907 0.3090
##
## sigma^2 estimated as 0.9347: log likelihood=-50.64
## AIC=107.28 AICc=108.03 BIC=112.04arima_austa_001.for <- forecast(arima_austa_001, h = 10)
autoplot(arima_austa_001.for)
arima_austa_000 <- Arima(austa, order = c(0,0,0))
arima_austa_000## Series: austa
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 3.5414
## s.e. 0.2972
##
## sigma^2 estimated as 3.27: log likelihood=-71.9
## AIC=147.8 AICc=148.17 BIC=150.97arima_austa_000.for <- forecast(arima_austa_000, h = 10)
autoplot(arima_austa_000.for)
8.8.E
arima_austa_021 <- Arima(austa-mean(austa), order = c(0,2,1))
arima_austa_021## Series: austa - mean(austa)
## ARIMA(0,2,1)
##
## Coefficients:
## ma1
## -0.7263
## s.e. 0.2277
##
## sigma^2 estimated as 0.0397: log likelihood=6.74
## AIC=-9.48 AICc=-9.09 BIC=-6.43arima_austa_021.for <- forecast(arima_austa_021, h = 10)
autoplot(arima_austa_021.for)
8.9.A
boxcox.usgdp <- BoxCox(usgdp, BoxCox.lambda(usgdp))8.9.B
usgdp.auto.arima <- auto.arima(boxcox.usgdp)
usgdp.auto.arima## Series: boxcox.usgdp
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## 0.2795 0.1208 0.1829
## s.e. 0.0647 0.0648 0.0202
##
## sigma^2 estimated as 0.03518: log likelihood=61.56
## AIC=-115.11 AICc=-114.94 BIC=-101.268.9.C
arima_usgdp_001 <- Arima(boxcox.usgdp, order = c(0,0,1))
arima_usgdp_001.AIC <- AIC(arima_usgdp_001)
arima_usgdp_011 <- Arima(boxcox.usgdp, order = c(0,1,1))
arima_usgdp_011.AIC <- AIC(arima_usgdp_011)
arima_usgdp_111 <- Arima(boxcox.usgdp, order = c(1,1,1))
arima_usgdp_111.AIC <- AIC(arima_usgdp_111)
arima_usgdp_010 <- Arima(boxcox.usgdp, order = c(0,1,0))
arima_usgdp_010.AIC <- AIC(arima_usgdp_010)
arima_usgdp_221 <- Arima(boxcox.usgdp, order = c(2,2,1))
arima_usgdp_221.AIC <- AIC(arima_usgdp_221)8.9.D
allAic<- c(arima_usgdp_001.AIC,arima_usgdp_011.AIC,arima_usgdp_111.AIC,arima_usgdp_010.AIC,arima_usgdp_221.AIC)
names(allAic) <- c("001","011","111","010","200")
allAic## 001 011 111 010 200
## 1553.11358 -10.79300 -86.42836 53.03936 -109.15960checkresiduals(arima_usgdp_111)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)
## Q* = 16.914, df = 6, p-value = 0.009605
##
## Model df: 2. Total lags used: 88.9.E
arima_usgdp_111.for <- forecast(arima_usgdp_111, h = 12)
autoplot(arima_usgdp_111.for)
8.9.F
ets.usgdp <- ets(usgdp)
AIC(ets.usgdp)## [1] 3072.303ets.usgdp.for <- forecast(ets.usgdp, h = 12)
autoplot(ets.usgdp.for)
8.10.A
autoplot(austourists)
8.10.B
Acf(austourists)
8.10.C
Pacf(austourists)
8.10.D
ggtsdisplay(diff(austourists, lag=4, differences=1))
8.10.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.68.10.F
beta*Y(t-1) + e_t
Y(t-2) + e_t
8.11.A
ma.usmelec<-ma(usmelec, order=12)
autoplot(ma.usmelec)
8.11.B
boxcox.usmelec <- BoxCox(usmelec,BoxCox.lambda(usmelec))
autoplot(boxcox.usmelec)
8.11.C
ggtsdisplay(boxcox.usmelec)
boxcox.usmelec.diff <- diff(diff(boxcox.usmelec, lag = 12))
ggtsdisplay(boxcox.usmelec.diff)
8.11.D
arima_usmelec_001 <- Arima(boxcox.usmelec, order = c(0,0,1), seasonal=c(2,1,1))
arima_usmelec_001.AIC <- AIC(arima_usmelec_001)
arima_usmelec_011 <- Arima(boxcox.usmelec, order = c(0,1,1), seasonal=c(2,1,0))
arima_usmelec_011.AIC <- AIC(arima_usmelec_011)
arima_usmelec_111 <- Arima(boxcox.usmelec, order = c(1,1,1), seasonal=c(2,1,1))
arima_usmelec_111.AIC <- AIC(arima_usmelec_111)
arima_usmelec_010 <- Arima(boxcox.usmelec, order = c(0,1,0), seasonal=c(1,1,0))
arima_usgdp_010.AIC <- AIC(arima_usmelec_010)
arima_usmelec_221 <- Arima(boxcox.usmelec, order = c(2,1,1), seasonal=c(2,1,1))
arima_usmelec_221.AIC <- AIC(arima_usmelec_221)
allAic<- c(arima_usmelec_001.AIC,arima_usmelec_011.AIC,arima_usmelec_111.AIC,arima_usgdp_010.AIC,arima_usmelec_221.AIC)
names(allAic) <- c("001","011","111","010","221")
allAic## 001 011 111 010 221
## -4823.332 -4980.303 -5082.718 -4847.659 -5081.2538.11.E
summary(arima_usmelec_111)## Series: boxcox.usmelec
## ARIMA(1,1,1)(2,1,1)[12]
##
## Coefficients:
## ar1 ma1 sar1 sar2 sma1
## 0.3888 -0.8265 0.0403 -0.0958 -0.8471
## s.e. 0.0630 0.0375 0.0555 0.0531 0.0341
##
## sigma^2 estimated as 1.274e-06: log likelihood=2547.36
## AIC=-5082.72 AICc=-5082.54 BIC=-5057.76
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set -6.052081e-05 0.001107711 0.0008388277 -0.003631173 0.05030025
## MASE ACF1
## Training set 0.5592102 0.01175029checkresiduals(arima_usmelec_111)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(2,1,1)[12]
## Q* = 28.012, df = 19, p-value = 0.08319
##
## Model df: 5. Total lags used: 248.11.F
autoplot(forecast(arima_usmelec_111, h=(15*12)))
8.12.A
autoplot(mcopper)
boxcox.mcopper <- BoxCox(mcopper, BoxCox.lambda(mcopper))8.12.B
mcopper.auto.arima <- auto.arima(boxcox.mcopper )
mcopper.auto.arima## Series: boxcox.mcopper
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## 0.3720
## s.e. 0.0388
##
## sigma^2 estimated as 0.04997: log likelihood=45.05
## AIC=-86.1 AICc=-86.08 BIC=-77.438.12.C
arima_mcopper_001 <- Arima(boxcox.mcopper, order = c(0,0,1))
arima_mcopper_001.AIC <- AIC(arima_mcopper_001)
arima_mcopper_011 <- Arima(boxcox.mcopper, order = c(0,1,1))
arima_mcopper_011.AIC <- AIC(arima_mcopper_011)
arima_mcopper_111 <- Arima(boxcox.mcopper, order = c(1,1,1))
arima_mcopper_111.AIC <- AIC(arima_mcopper_111)
arima_mcopper_010 <- Arima(boxcox.mcopper, order = c(0,1,0))
arima_mcopper_010.AIC <- AIC(arima_mcopper_010)
arima_mcopper_221 <- Arima(boxcox.mcopper, order = c(2,2,1))
arima_mcopper_221.AIC <- AIC(arima_mcopper_221)
allAic<- c(arima_mcopper_001.AIC,arima_mcopper_011.AIC,arima_mcopper_111.AIC,arima_mcopper_010.AIC,arima_mcopper_221.AIC)
names(allAic) <- c("001","011","111","010","220")
allAic ## 001 011 111 010 220
## 1784.89050 -86.09690 -84.10448 -15.69837 -75.579408.12.D
checkresiduals(arima_mcopper_011)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)
## Q* = 22.913, df = 23, p-value = 0.4659
##
## Model df: 1. Total lags used: 248.12.E
autoplot(forecast(arima_mcopper_011, h = 24)) 
8.12.F
ets.mcopper <- ets(mcopper)
AIC(ets.mcopper)## [1] 8052.038ets.mcopper.for <- forecast(ets.mcopper, h = 24)
autoplot(ets.mcopper.for)
8.13.A
autoplot(qauselec)
boxcox.qauselec <- BoxCox(qauselec, BoxCox.lambda(mcopper))8.13.B
ggtsdisplay(boxcox.qauselec)
ggtsdisplay(diff(boxcox.qauselec, differences = 2) %>% diff(4))
8.13.C
arima_qauselec_001 <- Arima(boxcox.qauselec, order = c(0,0,1), seasonal=c(1,1,2))
arima_qauselec_001.AIC <- AIC(arima_qauselec_001)
arima_qauselec_011 <- Arima(boxcox.qauselec, order = c(0,1,1), seasonal=c(0,1,2))
arima_qauselec_011.AIC <- AIC(arima_qauselec_011)
arima_qauselec_111 <- Arima(boxcox.qauselec, order = c(1,1,1), seasonal=c(1,0,2))
arima_qauselec_111.AIC <- AIC(arima_qauselec_111)
arima_qauselec_010 <- Arima(boxcox.qauselec, order = c(0,1,0), seasonal=c(0,1,2))
arima_qauselec_010.AIC <- AIC(arima_qauselec_010)
arima_qauselec_221 <- Arima(boxcox.qauselec, order = c(2,2,1), seasonal=c(1,1,2))
arima_qauselec_221.AIC <- AIC(arima_qauselec_221)
allAic<- c(arima_qauselec_001.AIC,arima_qauselec_011.AIC,arima_qauselec_111.AIC,arima_qauselec_010.AIC,arima_qauselec_221.AIC)
names(allAic) <- c("001","011","111","010","220")
allAic ## 001 011 111 010 220
## -727.4932 -796.5163 -805.2229 -768.1753 -778.37828.13.D
checkresiduals(arima_qauselec_111)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,0,2)[4]
## Q* = 11.079, df = 3, p-value = 0.01131
##
## Model df: 5. Total lags used: 88.13.E
autoplot(forecast(arima_qauselec_111, h = 24))
8.13.F
ets.qauselec <- ets(boxcox.qauselec)
AIC(ets.qauselec)## [1] -281.0346ets.qauselec.for <- forecast(ets.qauselec, h = 24)
autoplot(ets.qauselec.for)
8.14
stlf.qauselec <- stl(qauselec, t.window=13, s.window="periodic")
stlf.qauselec.ses <- seasadj(stlf.qauselec)
stlf.qauselec.ses.arima <- stlf(stlf.qauselec.ses, method="arima")
autoplot(stlf.qauselec.ses.arima)
8.15.A
auto.arima(myts)## Series: myts
## ARIMA(1,0,1)(1,1,1)[12] with drift
##
## Coefficients:
## ar1 ma1 sar1 sma1 drift
## 0.9047 -0.2277 0.3904 -0.8288 0.7500
## s.e. 0.0268 0.0635 0.0819 0.0536 0.1447
##
## sigma^2 estimated as 169.2: log likelihood=-1471.51
## AIC=2955.02 AICc=2955.25 BIC=2978.498.16.A
autoplot(sheep)
8.16.B
Arima(3,1,0)
8.16.C
Acf(diff(sheep))
Pacf(diff(sheep))
8.16.D
nums <- c(1648,1665,1627,1791,1797)
phi1 <- .42
phi2 <- -.2
phi3 <- -.3
#1940
i <- 5
oneFor <- nums[i]+phi1*(nums[i]-nums[i-1])+phi2*(nums[i-1]-nums[i-2])+phi3*(nums[i-2]-nums[i-3])
oneFor## [1] 1778.12#1941
nums <-c(nums,oneFor)
i <- 6
twoFor <- nums[i]+phi1*(nums[i]-nums[i-1])+phi2*(nums[i-1]-nums[i-2])+phi3*(nums[i-2]-nums[i-3])
twoFor## [1] 1719.79#1942
nums <-c(nums,twoFor)
i <- 7
threeFor <- nums[i]+phi1*(nums[i]-nums[i-1])+phi2*(nums[i-1]-nums[i-2])+phi3*(nums[i-2]-nums[i-3])
threeFor## [1] 1697.2688.16.E
forecast(Arima(sheep, order=c(3,1,0)))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1940 1777.996 1687.454 1868.537 1639.524 1916.467
## 1941 1718.869 1561.544 1876.193 1478.261 1959.476
## 1942 1695.985 1494.151 1897.819 1387.306 2004.664
## 1943 1704.067 1482.975 1925.160 1365.935 2042.199
## 1944 1730.084 1499.953 1960.215 1378.129 2082.039
## 1945 1746.371 1508.361 1984.381 1382.366 2110.376
## 1946 1745.518 1495.743 1995.292 1363.520 2127.515
## 1947 1733.953 1468.206 1999.699 1327.529 2140.377
## 1948 1724.299 1442.092 2006.506 1292.701 2155.898
## 1949 1722.829 1426.874 2018.784 1270.204 2175.4538.17.A
autoplot(bicoal)
8.17.B
Arima(4,0,0)
8.17.C
Acf(bicoal)
Pacf(bicoal)
8.17.D
nums <- c(467,512,534,552,545)
phi1 <- .83
phi2 <- -.34
phi3 <- .55
phi4 <- -.38
#1964
i <- 5
oneFor <- 162 + phi1*nums[i]+phi2*nums[i-1]+phi3*nums[i-2] + phi4*nums[i-3]
oneFor## [1] 525.81#1965
nums <-c(nums,oneFor)
i <- 6
twoFor <- 162 + phi1*nums[i]+phi2*nums[i-1]+phi3*nums[i-2] + phi4*nums[i-3]
twoFor## [1] 513.8023#1965
nums <-c(nums,twoFor)
i <- 6
threeFor <- 162 + phi1*nums[i]+phi2*nums[i-1]+phi3*nums[i-2] + phi4*nums[i-3]
threeFor## [1] 513.80238.17.E
forecast(Arima(bicoal, order=c(4,0,0)))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1969 527.6291 459.8804 595.3779 424.0164 631.2419
## 1970 517.1923 429.0014 605.3832 382.3160 652.0686
## 1971 503.8051 412.4786 595.1315 364.1334 643.4768
## 1972 489.2909 390.4658 588.1160 338.1510 640.4309
## 1973 482.6041 379.6439 585.5643 325.1400 640.0682
## 1974 478.5776 375.5783 581.5770 321.0537 636.1016
## 1975 474.5657 371.4760 577.6554 316.9036 632.2278
## 1976 474.4001 371.2205 577.5796 316.6006 632.1996
## 1977 475.9463 372.5126 579.3801 317.7581 634.1346
## 1978 476.5973 372.9772 580.2174 318.1240 635.07068.18.A
pall <- Quandl("JOHNMATT/PALL", type="ts", start_date="1982-01-01")## Warning in Quandl("JOHNMATT/PALL", type = "ts", start_date = "1982-01-01"): Type
## 'ts' does not support frequency 365. Returning zoo.pall <- pall$`New York 9:30`
pall <- apply.monthly(as.xts(pall), mean,na.rm=TRUE)
autoplot(pall)
8.18.B
boxcox.pall <- BoxCox(pall, BoxCox.lambda(pall))
auto.arima.pal <- auto.arima(boxcox.pall)
auto.arima.pal## Series: boxcox.pall
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.2926 0.0364
## s.e. 0.0508 0.0187
##
## sigma^2 estimated as 0.07264: log likelihood=-36.23
## AIC=78.47 AICc=78.54 BIC=908.18.C
checkresiduals(auto.arima.pal)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1) with drift
## Q* = 8.8348, df = 8, p-value = 0.3564
##
## Model df: 2. Total lags used: 108.18.D
autoplot(forecast(boxcox.pall, h= 36))
8.18.E
ets.pal<- ets(boxcox.pall)
ets.pal## ETS(M,Ad,N)
##
## Call:
## ets(y = boxcox.pall)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.2001
## phi = 0.8026
##
## Initial states:
## l = 7.215
## b = 0.0535
##
## sigma: 0.0229
##
## AIC AICc BIC
## 1122.000 1122.248 1145.0798.18.F
checkresiduals(ets.pal)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,Ad,N)
## Q* = 14.359, df = 5, p-value = 0.01348
##
## Model df: 5. Total lags used: 108.18.G
autoplot(forecast(ets.pal, h = 36))
8.18.H
AIC(ets.pal)## [1] 1122AIC(auto.arima.pal)## [1] 78.46904