library(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.5 v fma 2.4
## v forecast 8.15 v expsmooth 2.3## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'forecast' was built under R version 4.0.5## 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(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.0 v forcats 0.5.1
## v purrr 0.3.4## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'purrr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'stringr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(tseries)## Warning: package 'tseries' was built under R version 4.0.5library(Quandl)## Warning: package 'Quandl' was built under R version 4.0.5## Loading required package: xts## Warning: package 'xts' was built under R version 4.0.5## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.0.5##
## 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, lastpigs <- fma::pigs
tsdisplay(pigs)
plot(pigs)
sespigs = ses(pigs, h=4)sespigs$model## 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.665autoplot(sespigs)
sespigs$upper[1, "95%"]## 95%
## 119020.8sespigs$lower[1, "95%"]## 95%
## 78611.97s = sd(sespigs$residuals)
sespigs$mean[1] + 1.96*s## [1] 118952.8sespigs$mean[1] - 1.96*s## [1] 78679.97#q2
SES <- function(y, alpha, l0){
y_hat <- l0
for(index in 1:length(y)){
y_hat <- alpha*y[index] + (1 - alpha)*y_hat
}
cat("Next Observation of SES",
as.character(y_hat),
sep = "\n")
}
alpha <- sespigs$model$par[1]
l0 <- sespigs$model$par[2]
SES(pigs, alpha = alpha, l0 = l0)## Next Observation of SES
## 98816.4061115907writeLines(paste(
"Next observation of SES", as.character(sespigs$mean[1])
))## Next observation of SES 98816.4061115907sespigs = ses(pigs, h=4)
sespigs## 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.8The point observation is the same from SES function and my own function
#Q3
SES <- function(pars = c(alpha, l0), y){
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)
}
# compare ses and SES using pigs data
opt_SES_pigs <- optim(par = c(0.5, pigs[1]), y = pigs, fn = SES)
writeLines(paste(
"Optimal parameters",
"\n",
as.character(opt_SES_pigs$par[1]),
", ",
as.character(opt_SES_pigs$par[2]),
sep = ""
))## Optimal parameters
## 0.299008094014243, 76379.2653476235writeLines(paste(
"Parameters",
"\n",
as.character(sespigs$model$par[1]),
", ",
as.character(sespigs$model$par[2]),
sep = ""))## Parameters
## 0.297148833372095, 77260.0561458528The alpha value is a little larger using the optimal ses function and the the los were close as well.
#Q4
SES <- function(init_pars, data){
fc_next <- 0
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
}
fc_next <<- y_hat
return(SSE)
}
.
optim_pars <- optim(par = init_pars, data = data, fn = SSE)
return(list(
Next_observation_forecast = fc_next,
alpha = optim_pars$par[1],
l0 = optim_pars$par[2]
))
}
print("Next obsfc by ses")## [1] "Next obsfc by ses"sespigs$mean[1]## [1] 98816.41print("alpha by ses")## [1] "alpha by ses"sespigs$model$par[1]## alpha
## 0.2971488print("l0 by ses")## [1] "l0 by ses"sespigs$model$par[2]## l
## 77260.06#Q5
autoplot(books)
sespaper= ses(books[,"Paperback"], h=4)
seshard = ses(books[,"Hardcover"], h=4)
autoplot(sespaper)
autoplot(seshard)
spaper = sqrt(mean(sespaper$residuals^2))
shard = sqrt(mean(seshard$residuals^2))#Q6
#6A
holtpaper = holt(books[,"Paperback"], h=4)
holthard = holt(books[,"Hardcover"], h=4)#6B
hpaper = sqrt(mean(holtpaper$residuals^2))
hhard = sqrt(mean(holthard$residuals^2))#6C
autoplot(holtpaper)
autoplot(sespaper)
autoplot(holthard)
autoplot(seshard)
#Holt’s method is better because it captures the upward trend better
#6D
holtpaper$upper[1, "95%"]## 95%
## 275.0205holtpaper$lower[1, "95%"]## 95%
## 143.913holtpaper$mean[1] +1.96*hpaper## [1] 270.4951holtpaper$mean[1] -1.96*hpaper## [1] 148.4384holthard$upper[1, "95%"]## 95%
## 307.4256holthard$lower[1, "95%"]## 95%
## 192.9222holthard$mean[1] + 1.96*hhard## [1] 303.4733holthard$mean[1] - 1.96*hhard## [1] 196.8745#7
## Time-Series [1:94] from 1900 to 1993: 277 315 315 321 315 ...## Time Series:
## Start = 1900
## End = 1905
## Frequency = 1
## [1] 276.79 315.42 314.87 321.25 314.54 317.92
## RMSE when using holt function## [1] 26.58219## RMSE when using holt function with damped option## [1] 26.54019## RMSE when using holt function with Box-Cox transformation## [1] 1.032217## RMSE when using holt function with damped option and Box-Cox transformation## [1] 1.039187. The best model was the Box-Cox transformation with Holt’s linear method. It gave pint forecasts and prediction intervals witha damping effect. Holt’s method went into the negatiive the damped method keeps from going to 0 but doesnt follow the overal downward trend.
#8
## [1] 14.72762## [1] 14.94306
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 42.932, df = 7, p-value = 3.437e-07
##
## Model df: 17. Total lags used: 24
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4556121 8.681456 6.24903 0.2040939 3.151257 0.3916228
## Test set 94.7346169 111.911266 94.73462 24.2839784 24.283978 5.9369594
## ACF1 Theil's U
## Training set -0.01331859 NA
## Test set 0.60960299 1.90013
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03021223 9.107356 6.553533 0.001995484 3.293399 0.4107058
## Test set 78.34068365 94.806617 78.340684 19.945024968 19.945025 4.9095618
## ACF1 Theil's U
## Training set 0.02752875 NA
## Test set 0.52802701 1.613903#9
fc_stl_ets_retail_train <- ts_retail_train %>%
stlm(
s.window = 13,
robust = TRUE,
method = "ets",
lambda = BoxCox.lambda(ts_retail_train)
) %>%
forecast(
h = 36,
lambda = BoxCox.lambda(ts_retail_train)
)## Warning in InvBoxCox(fcast$mean, lambda, biasadj, fcast): biasadj information
## not found, defaulting to FALSE.autoplot(fc_stl_ets_retail_train)
accuracy(fc_stl_ets_retail_train, ts_retail_test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.6782982 8.583559 5.918078 -0.3254076 2.913104 0.3708823
## Test set 82.1015276 98.384220 82.101528 21.0189982 21.018998 5.1452516
## ACF1 Theil's U
## Training set 0.02704667 NA
## Test set 0.52161725 1.679783fc_stl_ets_retail_train_without_tr <-
ts_retail_train %>%
stlm(
s.window = 13,
robust = TRUE,
method = "ets"
) %>%
forecast(h = 36)
autoplot(fc_stl_ets_retail_train_without_tr)
accuracy(fc_stl_ets_retail_train_without_tr,
ts_retail_test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.5020795 10.00826 6.851597 -0.391432 3.489759 0.4293853
## Test set 74.2529959 91.04491 74.252996 18.837766 18.837766 4.6533890
## ACF1 Theil's U
## Training set 0.09741266 NA
## Test set 0.48917501 1.549271# I couldn't expect it because when I also used transformation, the accuracy of training set was better. In fact, the actual values in forecast horizon increased exponentially. Without using transformation, the forecast may reflect that bigger values have bigger variation. ETS forecasting after STL decomposition 'without' Box-Cox transformation yielded better result than when ETS(A, Ad, M) or ETS(A, A, M) was used. #10
#10A
autoplot(ukcars)
#TThere is a downward trend until around 1982, then trends upwards until 2000, where it dips and then plateus.#10B
stlcars = stl(ukcars, s.window=4, robust=TRUE)
autoplot(stlcars)
seas.stlcars = seasadj(stlcars)
autoplot(seas.stlcars)
#10C
add.damp.stl = stlf(seas.stlcars, etsmodel= "AAN", damped=TRUE, h=8)
autoplot(add.damp.stl)
#10D
holtfc = stlf(seas.stlcars, etsmodel="AAN", damped=FALSE, h=8)
autoplot(holtfc)
#10E
etsfc = ets(ukcars)
summary(etsfc)## 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
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
## ACF1
## Training set 0.02573334#10F
#Accuracy of additive damped trend method
accuracy(add.damp.stl)## ME RMSE MAE MPE MAPE MASE
## Training set 1.854493 21.0918 16.59365 0.3029367 5.320469 0.5380042
## ACF1
## Training set -0.005454679accuracy(holtfc)## ME RMSE MAE MPE MAPE MASE
## Training set -0.354125 21.1403 16.47804 -0.4358911 5.319257 0.5342557
## ACF1
## Training set 0.004332171accuracy(etsfc)## ME RMSE MAE MPE MAPE MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
## ACF1
## Training set 0.02573334#10G
autoplot(add.damp.stl)
autoplot(holtfc)
autoplot(forecast(etsfc, h=8))
#2 are seasonally adjusted#10H
checkresiduals(add.damp.stl)## Warning in checkresiduals(add.damp.stl): 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,Ad,N)
## Q* = 48.348, df = 3, p-value = 1.796e-10
##
## Model df: 5. Total lags used: 8checkresiduals(holtfc)## Warning in checkresiduals(holtfc): 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* = 45.837, df = 4, p-value = 2.663e-09
##
## Model df: 4. Total lags used: 8checkresiduals(forecast(etsfc, h=8))
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,A)
## Q* = 18.324, df = 3, p-value = 0.0003772
##
## Model df: 6. Total lags used: 9#11
autoplot(visitors)
Upward trend line, looks like the seasonality increases propotional to the level #11B
test = window(visitors, start = c(2003,4))
train = window(visitors, end = c(2004, 3))
hwmult = hw(train, seasonal="multiplicative", h=24)
autoplot(hwmult)
#11C Multiplicative is needed because the seasonal variation in the data increases as the level of the series increases.
#11D
etsfc <- train %>% ets() %>% forecast(h=24)
autoplot(etsfc)
etsadd.BCfc <- train %>% ets(lambda = BoxCox.lambda(train), additive.only=TRUE) %>% forecast(h=24)
autoplot(etsadd.BCfc)
snaive = snaive(train, h=24)
autoplot(snaive)
BCstl.ets <- train %>% stlm(
lambda = BoxCox.lambda(train),
s.window = 13,
robust=TRUE,
method="ets") %>% forecast(h=24)
autoplot(BCstl.ets)
#11E
accuracy(hwmult, test)## ME RMSE MAE MPE MAPE MASE
## Training set -1.557919 14.70228 11.17152 -1.019484 4.452697 0.4268183
## Test set 15.090935 20.85700 16.04326 3.122066 3.359691 0.6129476
## ACF1 Theil's U
## Training set 0.09876254 NA
## Test set 0.21025230 0.2869292accuracy(etsfc, test)## ME RMSE MAE MPE MAPE MASE
## Training set 1.220635 14.90665 10.78433 0.2197614 4.018479 0.4120252
## Test set 31.260025 36.53961 31.26002 6.6962946 6.696295 1.1943180
## ACF1 Theil's U
## Training set -0.01726088 NA
## Test set 0.39874182 0.5090901accuracy(etsadd.BCfc, test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.373716 15.55474 11.35863 -0.1677915 4.070118 0.4339668
## Test set 12.013636 19.11730 15.31311 2.3827550 3.229648 0.5850514
## ACF1 Theil's U
## Training set -0.05004163 NA
## Test set 0.23875318 0.2521544accuracy(snaive, test)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 16.78326 31.24472 26.17395 6.837756 10.14381 1.000000 0.6512811
## Test set 48.30000 55.23646 48.30000 11.417380 11.41738 1.845346 0.6747630
## Theil's U
## Training set NA
## Test set 0.7596856accuracy(BCstl.ets, test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.8379729 13.88885 9.778233 -0.2528935 3.536527 0.3735864
## Test set 9.4686802 17.54139 13.069371 1.8664731 2.747569 0.4993274
## ACF1 Theil's U
## Training set -0.04339136 NA
## Test set -0.04401040 0.2405589#12
str(ausbeer)## Time-Series [1:218] from 1956 to 2010: 284 213 227 308 262 228 236 320 272 233 ...head(ausbeer)## Qtr1 Qtr2 Qtr3 Qtr4
## 1956 284 213 227 308
## 1957 262 228autoplot(ausbeer)
ausbeer_train <- subset(
ausbeer, end = length(ausbeer) - 12
)
ausbeer_test <- subset(
ausbeer, start = length(ausbeer) - 11
)
ets_ausbeer_train <- forecast(
ets(ausbeer_train), h = 12
)
snaive_ausbeer_train <- snaive(ausbeer_train, h = 12)
stlf_ausbeer_train <- stlf(
ausbeer_train,
h = 12,
s.window = 5,
robust = TRUE,
lambda = BoxCox.lambda(ausbeer_train))
accuracy(ets_ausbeer_train, ausbeer_test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.3466095 15.781973 11.979955 -0.064073435 2.864311 0.7567076
## Test set 0.1272571 9.620828 8.919347 0.009981598 2.132836 0.5633859
## ACF1 Theil's U
## Training set -0.1942276 NA
## Test set 0.3763918 0.1783972accuracy(snaive_ausbeer_train, ausbeer_test)## ME RMSE MAE MPE MAPE MASE
## Training set 3.336634 19.67005 15.83168 0.9044009 3.771370 1.0000000
## Test set -2.916667 10.80509 9.75000 -0.6505927 2.338012 0.6158537
## ACF1 Theil's U
## Training set 0.0009690785 NA
## Test set 0.3254581810 0.1981463accuracy(stlf_ausbeer_train, ausbeer_test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.5704833 13.61609 9.816091 0.1279204 2.334194 0.6200283
## Test set -4.0433841 10.11436 8.429747 -0.9670308 2.051915 0.5324605
## ACF1 Theil's U
## Training set -0.1131945 NA
## Test set 0.2865992 0.1689574# Without RMSE, all the other errors show that the best model is STL + ETS(M, Ad, N).
autoplot(stlf_ausbeer_train) +
autolayer(ausbeer_test)
# The forecasts are similar to real data.
forecast.models <- function(y, h){
m <- frequency(y)
y_train <- subset(
y, end = length(y) - m*h
)
y_test <- subset(
y, start = length(y) - m*h + 1
)
# try each model and forecast.
ets_y_train <- forecast(
ets(y_train), h = m*h
)
snaive_y_train <- snaive(y_train, h = m*h)
stlf_y_train <- stlf(
y_train,
h = m*h,
s.window = m + 1,
robust = TRUE
)
stlf_y_train_with_BoxCox <- stlf(
y_train,
h = m*h,
s.window = m + 1,
robust = TRUE,
lambda = BoxCox.lambda(y_train))
# combine forecasts to models variable
models <- list(ets_y_train,
snaive_y_train,
stlf_y_train,
stlf_y_train_with_BoxCox)
names(models) <- c("ets", "snaive", "stlf", "stlf_with_BoxCox")
# get accuracy for each model using test set
acc_ets <- accuracy(ets_y_train, y_test)
acc_snaive <- accuracy(snaive_y_train, y_test)
acc_stlf <- accuracy(stlf_y_train, y_test)
acc_stlf_with_BoxCox <- accuracy(stlf_y_train_with_BoxCox, y_test)
# combine accuracies to accuracies variable.
accuracies <- list(acc_ets,
acc_snaive,
acc_stlf,
acc_stlf_with_BoxCox)
names(accuracies) <- c("acc_ets", "acc_snaive", "acc_stlf", "acc_stlf_with_BoxCox")
# return output values
output <- list(y_train, y_test, m, models, accuracies)
names(output) <- c("train", "test", "m", "models", "accuracies")
return(output)
}
# bricksq data case
fc_bricksq <- forecast.models(bricksq, 3)
fc_bricksq$accuracies## $acc_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
##
## $acc_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
##
## $acc_stlf
## ME RMSE MAE MPE MAPE MASE
## Training set 1.318736 21.49082 14.41766 0.2855898 3.501287 0.3919528
## Test set 44.266931 50.20039 45.89301 10.0933739 10.487099 1.2476294
## ACF1 Theil's U
## Training set 0.1502759 NA
## Test set 0.1053316 1.31732
##
## $acc_stlf_with_BoxCox
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.42724 21.40466 14.30375 0.3185472 3.500221 0.3888560 0.1625098
## Test set 34.69481 40.13618 36.34649 7.7322091 8.132133 0.9881014 0.4746242
## Theil's U
## Training set NA
## Test set 1.04561# All errors show that the best model is seasonal naive method.
autoplot(fc_bricksq$models$snaive) +
autolayer(fc_bricksq$test)
# In about an year the forecasts were similar to real data, but after that the real data increased exponentially while the trend of forecasts didn't change. But real data were still in the 80% prediction interval.
fc_dole <- forecast.models(dole, 3)
fc_dole$accuracies## $acc_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
##
## $acc_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
##
## $acc_stlf
## ME RMSE MAE MPE MAPE MASE
## Training set -96.4811 7801.958 4530.06 -0.2719573 5.116149 0.1439341
## Test set 328005.9874 407547.190 328005.99 48.6435815 48.643581 10.4217690
## ACF1 Theil's U
## Training set 0.08037493 NA
## Test set 0.93373958 16.57033
##
## $acc_stlf_with_BoxCox
## ME RMSE MAE MPE MAPE MASE
## Training set 145.1468 6688.083 3659.404 0.1464869 3.82385 0.1162706
## Test set 205385.6547 268135.992 207317.238 29.3540384 29.85051 6.5871126
## ACF1 Theil's U
## Training set -0.09446256 NA
## Test set 0.93778748 10.68587# All errors show best model is seasonal naive method.
autoplot(fc_dole$models$snaive) +
autolayer(fc_dole$test)
# The forecasts were completely wrong. Even the best model couldn't predict such change.
fc_a10 <- forecast.models(a10, 3)
fc_a10$accuracies## $acc_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
##
## $acc_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
##
## $acc_stlf
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06704027 0.6809938 0.4297557 0.3200743 4.732989 0.4381286
## Test set 1.83324369 3.1118106 2.4570607 7.0963340 11.223584 2.5049311
## ACF1 Theil's U
## Training set -0.01970509 NA
## Test set 0.23918322 0.8089697
##
## $acc_stlf_with_BoxCox
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01686093 0.4244591 0.3098511 0.003515845 3.594686 0.3158879
## Test set 1.25869582 2.2242431 1.8428354 5.032385108 8.701050 1.8787390
## ACF1 Theil's U
## Training set -0.1780966 NA
## Test set 0.1502463 0.5964321# All errors show that the best model is STL + ETS(A, A, N) with BoxCox transformation model.
autoplot(fc_a10$models$stlf_with_BoxCox) +
autolayer(fc_a10$test)
#The best model could follow the general trend, but a little short of predicting more fastly increasing trend.
fc_h02 <- forecast.models(h02, 3)
fc_h02$accuracies## $acc_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
##
## $acc_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
##
## $acc_stlf
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001987465 0.04609626 0.03396307 -0.1760405 4.751272 0.5736778
## Test set 0.046520958 0.09174622 0.07544541 3.5302072 8.001263 1.2743653
## ACF1 Theil's U
## Training set 0.01953927 NA
## Test set 0.05661406 0.5085785
##
## $acc_stlf_with_BoxCox
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003628017 0.04230875 0.03068561 0.02682146 4.159581 0.5183175
## Test set 0.031706735 0.07574975 0.06381807 2.51742657 6.903027 1.0779653
## ACF1 Theil's U
## Training set -0.09600393 NA
## Test set -0.25297574 0.4385025# All errors show that the best model is STL + ETS(A, Ad, N) method.
autoplot(fc_h02$models$stlf_with_BoxCox) +
autolayer(fc_h02$test)
# The forecasts were similar to real data for the most part.
fc_usmelec <- forecast.models(usmelec, 3)
fc_usmelec$accuracies## $acc_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
##
## $acc_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
##
## $acc_stlf
## ME RMSE MAE MPE MAPE MASE
## Training set -0.07064178 6.659696 4.907696 -0.07729202 1.910234 0.5454076
## Test set -11.41131201 17.776488 15.918880 -3.66455560 4.779167 1.7691150
## ACF1 Theil's U
## Training set 0.1126383 NA
## Test set 0.5099758 0.5639709
##
## $acc_stlf_with_BoxCox
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1351692 6.474779 4.761143 -0.05026415 1.848092 0.5291207
## Test set -19.9969503 24.406877 20.951753 -6.06977735 6.315723 2.3284338
## ACF1 Theil's U
## Training set 0.08255205 NA
## Test set 0.60877348 0.7777035# Most of errors show the best model is ETS(M, A, M) method.
autoplot(fc_usmelec$models$ets) +
autolayer(fc_usmelec$test)
# Real data was within the prediction interval for the most part. #Chapter 8 #2
ggtsdisplay(ibmclose)
# ACF plot shows that the autocorrelation values are bigger than critical value and decrease slowly. Also, r1 is large(near to 1) and positive. It means that the IBM stock is, predictable using lagged values. PACF plot shows that there is a strong correlation between IBM stock data and their 1 lagged values. It means that IBM stock data can be predicted by 1 lagged values. To get stationary data, IBM stock data need differencing. Differencing can make non-staionary data stationary.transform.plot <- function(x){
lam <- BoxCox.lambda(x)
transformed <- BoxCox(x, lambda = lam)
tsdisplay(transformed)
return(paste("Box-Cox Lambda = ", round(lam, 4)))
}
data("usnetelec")
transform.plot(usnetelec)
## [1] "Box-Cox Lambda = 0.5168"#The series needs one differencing in order to be stationary
data("usgdp")
transform.plot(usgdp)
## [1] "Box-Cox Lambda = 0.3664"# The series needs one differencing in order to be stationary
data("mcopper")
transform.plot(mcopper)
## [1] "Box-Cox Lambda = 0.1919"#The series needs two steps of differencing in order to be stationary
data("enplanements")
transform.plot(enplanements)
## [1] "Box-Cox Lambda = -0.2269"#The series needs one differencing and seasonal differencing in order to be stationary
data("visitors")
transform.plot(visitors)
## [1] "Box-Cox Lambda = 0.2775"## [1] "Box-Cox Lambda = 0.2775"#5a
autoplot(austa)
fc_austa_autoarima <- forecast(
auto.arima(austa), h = 10
)
fc_austa_autoarima$model## 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# ARIMA(0, 1, 1) with drift model was chosen.
checkresiduals(fc_austa_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1) with drift
## Q* = 2.297, df = 5, p-value = 0.8067
##
## Model df: 2. Total lags used: 7# The residuals are like white noise.
autoplot(fc_austa_autoarima)
# 5b
fc_austa_arima.0.1.1 <- forecast(
Arima(austa, order = c(0, 1, 1)), h = 10
)
autoplot(fc_austa_arima.0.1.1)
fc_austa_arima.0.1.0 <- forecast(
Arima(austa, order = c(0, 1, 0)), h = 10
)
autoplot(fc_austa_arima.0.1.0)
# the forecasts of both models are like the result of naive forecast. Increasing trend isn't reflected in the forecasts.
fc_austa_arima.0.1.1$upper - fc_austa_arima.0.1.0$upper## Time Series:
## Start = 2016
## End = 2025
## Frequency = 1
## fc_austa_arima.0.1.1$upper.80% fc_austa_arima.0.1.1$upper.95%
## 2016 0.1138117 0.09101057
## 2017 0.2039418 0.22885273
## 2018 0.2527212 0.30345426
## 2019 0.2886706 0.35843413
## 2020 0.3180942 0.40343373
## 2021 0.3434784 0.44225553
## 2022 0.3660754 0.47681466
## 2023 0.3866116 0.50822207
## 2024 0.4055495 0.53718500
## 2025 0.4232034 0.56418442fc_austa_arima.0.1.0$lower - fc_austa_arima.0.1.1$lower## Time Series:
## Start = 2016
## End = 2025
## Frequency = 1
## fc_austa_arima.0.1.0$lower.80% fc_austa_arima.0.1.0$lower.95%
## 2016 -0.199956384 -0.22275751
## 2017 -0.109826244 -0.08491535
## 2018 -0.061046922 -0.01031382
## 2019 -0.025097519 0.04466605
## 2020 0.004326137 0.08966565
## 2021 0.029710347 0.12848745
## 2022 0.052307350 0.16304658
## 2023 0.072843552 0.19445399
## 2024 0.091781392 0.22341692
## 2025 0.109435361 0.25041634# But prediction interval of ARIMA(0, 1, 1) model was generally larger than the one of ARIMA(0, 1, 0) model. I think that it is because of one more error term in ARIMA(0, 1, 1) model.#5c
fc_austa_arima.2.1.3.drift <- forecast(
Arima(austa, order = c(2, 1, 3), include.drift = TRUE),
h = 10
)
autoplot(fc_austa_arima.2.1.3.drift)
# The forecasts are increasing, but the speed of the increase is decreasing.
drift_austa <- fc_austa_arima.2.1.3.drift$model$coef[6]
fc_austa_arima.2.1.3.nodrift <- fc_austa_arima.2.1.3.drift$mean - drift_austa*seq_len(10)
autoplot(fc_austa_arima.2.1.3.drift) +
autolayer(fc_austa_arima.2.1.3.nodrift)
#5d
fc_austa_arima.0.0.1.const <- forecast(
Arima(
austa, order = c(0, 0, 1), include.constant = TRUE
),
h = 10
)
autoplot(fc_austa_arima.0.0.1.const)
fc_austa_arima.0.0.0.const <- forecast(
Arima(austa, order = c(0, 0, 0), include.constant = TRUE),
h = 10
)
autoplot(fc_austa_arima.0.0.0.const)
# All of the forecasts are the mean of the data history. It is like the result of mean method.#5e
fc_austa_arima.0.2.1 <- forecast(
Arima(austa, order = c(0, 2, 1)),
h = 10
)
autoplot(fc_austa_arima.0.2.1)
# the forecasts show increasing trend. PI is being larger for the farther future forecast.#6a
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
plot(y[i])
#6B
ar1 <- function(phi1){
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <- phi1*y[i-1] + e[i]
}
return(y)
}
autoplot(ar1(0.3), series = "0.3") +
geom_line(size = 1, colour = "red") +
autolayer(y, series = "0.6", size = 1) +
autolayer(ar1(0.9), size = 1, series = "0.9") +
ylab("AR(1) models") +
guides(colour = guide_legend(title = "Phi1"))
#The variation in y increases as phi increases.#6C
ma1 <- function(theta1){
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <- theta1*e[i-1] + e[i]
}
return(y)
}#6D
autoplot(ma1(0.3), series = "0.3") +
geom_line(size = 1, colour = "red") +
autolayer(y, series = "0.6", size = 1) +
autolayer(ma1(0.9), size = 1, series = "0.9") +
ylab("MA(1) models") +
guides(colour = guide_legend(title = "Theta1"))
#as theta increases, variation of y increases, just like phi. #6E
MYarima11 <- ts(numeric(50))
e <- rnorm(50)
for(i in 2:50){
MYarima11[i] <- 0.6*MYarima11[i-1] + 0.6*e[i-1] + e[i]
}#6F
MYarima2 <- ts(numeric(50))
e <- rnorm(50)
for(i in 3:50){
MYarima2[i] <- -0.8*MYarima2[i-1] + 0.3*MYarima2[i-2] + e[i]
}#6G
autoplot(MYarima11, series = "ARMA(1, 1)") +
autolayer(MYarima2, series = "AR(2)") +
ylab("y") +
guides(colour = guide_legend(title = "ARIMA Method"))
#7
# a
autoplot(wmurders)
# It looked like the data don't need seasonal differencing or Box-Cox transformation.
autoplot(diff(wmurders))
# It looked like 1 more differencing would be needed to make the data stationary. Differenced data slowly go to minus infinity.
ndiffs(wmurders)## [1] 2# ndiffs function shows that the data need 2 differencing.
autoplot(diff(wmurders, differences = 2))
kpss.test(diff(wmurders, differences = 2))## Warning in kpss.test(diff(wmurders, differences = 2)): p-value greater than
## printed p-value##
## KPSS Test for Level Stationarity
##
## data: diff(wmurders, differences = 2)
## KPSS Level = 0.045793, Truncation lag parameter = 3, p-value = 0.1# twice differencing made the data stationary.
diff(wmurders, differences = 2) %>% ggtsdisplay()
# PACF is decaying. And there are significant spikes at lag 1, and 2 in the ACF, but none beyond lag 2. If the data can be modelled by ARIMA(0, 2, q) or ARIMA(p, 2, 0), I'm going to model the data by ARIMA(0, 2, 2).
# b
# ARIMA model of the data includes twice differencing. If there is a constant in the model, twice integrated contant will yield quadratic trend, which is dangerous for forecasting. Therefore I won't include a constant in the model.
# c
# (1 - B)^2*yt = (1 + theta1*B + theta2*B^2)*et
# d
wmurders_arima.0.2.2 <- Arima(wmurders,
order = c(0, 2, 2))
checkresiduals(wmurders_arima.0.2.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,2)
## Q* = 11.764, df = 8, p-value = 0.1621
##
## Model df: 2. Total lags used: 10# The residuals of the model can be thought of as white noise series. A little sorry that they aren't normally distributed. But it is satisfactory to get them.
# e
fc_wmurders_arima.0.2.2 <- forecast(
wmurders_arima.0.2.2, h = 3
)
# forecasts by Arima function
fc_wmurders_arima.0.2.2$mean## Time Series:
## Start = 2005
## End = 2007
## Frequency = 1
## [1] 2.480525 2.374890 2.269256# get forecasts by manual calculation
fc_wmurders_arima.0.2.2$model## Series: wmurders
## ARIMA(0,2,2)
##
## Coefficients:
## ma1 ma2
## -1.0181 0.1470
## s.e. 0.1220 0.1156
##
## sigma^2 estimated as 0.04702: log likelihood=6.03
## AIC=-6.06 AICc=-5.57 BIC=-0.15# formula
# (1 - B)^2*yt = (1 - 1.0181*B + 0.1470*B^2)*et
# yt = 2yt-1 - yt-2 + et - 1.0181*et-1 + 0.1470*et-2
years <- length(wmurders)
e <- fc_wmurders_arima.0.2.2$residuals
fc1 <- 2*wmurders[years] - wmurders[years - 1] - 1.0181*e[years] + 0.1470*e[years - 1]
fc2 <- 2*fc1 - wmurders[years] + 0.1470*e[years]
fc3 <- 2*fc2 - fc1
# forecasts by manual calculation
c(fc1, fc2, fc3)## [1] 2.480523 2.374887 2.269252# the forecasts are almost similar to the ones got by Arima function.
# f
autoplot(fc_wmurders_arima.0.2.2)
# g
fc_wmurders_autoarima <- forecast(
auto.arima(wmurders), h = 3
)
# Without RMSE, all errors show that ARIMA(0, 2, 2) is better than ARIMA(1, 2, 1).
accuracy(fc_wmurders_arima.0.2.2)## ME RMSE MAE MPE MAPE MASE
## Training set -0.0113461 0.2088162 0.1525773 -0.2403396 4.331729 0.9382785
## ACF1
## Training set -0.05094066accuracy(fc_wmurders_autoarima)## ME RMSE MAE MPE MAPE MASE
## Training set -0.01065956 0.2072523 0.1528734 -0.2149476 4.335214 0.9400996
## ACF1
## Training set 0.02176343# try using auto.arima function with stepwise and approximation options false.
fc_wmurders_autoarima2 <- forecast(
auto.arima(wmurders, stepwise = FALSE, approximation = FALSE),
h = 3
)
# It is ARIMA(0, 2, 3) model.
accuracy(fc_wmurders_autoarima2)## ME RMSE MAE MPE MAPE MASE
## Training set -0.01336585 0.2016929 0.1531053 -0.3332051 4.387024 0.9415259
## ACF1
## Training set -0.03193856# In this case, some errors were better while others were worse. I'll check residuals and ACF, PACF plots.
ggtsdisplay(diff(wmurders, differences = 2))
# It looked like that the data are similar to ARIMA(0, 2, 2) rather than ARIMA(0, 2, 3).
checkresiduals(fc_wmurders_arima.0.2.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,2)
## Q* = 11.764, df = 8, p-value = 0.1621
##
## Model df: 2. Total lags used: 10checkresiduals(fc_wmurders_autoarima2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,3)
## Q* = 10.706, df = 7, p-value = 0.152
##
## Model df: 3. Total lags used: 10# almost similar residuals.
# Therefore I'll choose ARIMA(0, 2, 2).#8
usgdparima = auto.arima(usgdp)
usgdparima## Series: usgdp
## ARIMA(2,2,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.1228 0.3106 -0.5835 -0.3669
## s.e. 0.2869 0.0872 0.3004 0.2862
##
## sigma^2 estimated as 1604: log likelihood=-1199.57
## AIC=2409.13 AICc=2409.39 BIC=2426.43usgdpets = ets(usgdp)
usgdpets## ETS(A,A,N)
##
## Call:
## ets(y = usgdp)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.278
##
## Initial states:
## l = 1557.4589
## b = 18.6862
##
## sigma: 41.8895
##
## AIC AICc BIC
## 3072.303 3072.562 3089.643checkresiduals(usgdparima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,2,2)
## Q* = 8.6247, df = 4, p-value = 0.0712
##
## Model df: 4. Total lags used: 8usgdparima %>% forecast(h=8) %>% autoplot()
usgdpets %>% forecast(h=8) %>% autoplot
#10
autoplot(austourists)
#A
acf(austourists)
#B
pacf(austourists)
#C
ggtsdisplay(diff(austourists, lag=4))
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#11
mausmelec = ma(usmelec, order=12, centre = TRUE)
autoplot(mausmelec)
ndiffs(usmelec)## [1] 1elecarima = auto.arima(usmelec)
elecarima## 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.22checkresiduals(elecarima)
##
## 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: 24elecfc = forecast(elecarima, h= 12*15)
autoplot(elecfc)
#13
autoplot(qauselec)
# The data need Box-Cox transformation to make the variations evenly over time.
lambda_qauselec <- BoxCox.lambda(qauselec)# b
nsdiffs(qauselec)## [1] 1ndiffs(qauselec)## [1] 1# The data need 1 seasonal differencing.
kpss.test(diff(qauselec, lag = 4))##
## KPSS Test for Level Stationarity
##
## data: diff(qauselec, lag = 4)
## KPSS Level = 0.39382, Truncation lag parameter = 4, p-value = 0.07982# The data doesn't need first differencing.
# c
ggtsdisplay(diff(
BoxCox(qauselec, lambda_qauselec), lag = 4
))
ggtsdisplay(diff(diff(
BoxCox(qauselec, lambda_qauselec), lag = 4
)))
qauselec_arima0.0.1.0.1.1.4 <- Arima(
qauselec, lambda = lambda_qauselec,
order = c(0, 0, 1), seasonal = c(0, 1, 1)
)
qauselec_arima0.0.1.0.1.1.4## Series: qauselec
## ARIMA(0,0,1)(0,1,1)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ma1 sma1
## 0.6495 0.2738
## s.e. 0.0746 0.0642
##
## sigma^2 estimated as 0.03639: log likelihood=51.5
## AIC=-97.01 AICc=-96.89 BIC=-86.91qauselec_arima0.1.1.0.1.1.4 <- Arima(
qauselec, lambda = lambda_qauselec,
order = c(0, 1, 1), seasonal = c(0, 1, 1)
)
qauselec_arima0.1.1.0.1.1.4## Series: qauselec
## ARIMA(0,1,1)(0,1,1)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ma1 sma1
## -0.4976 -0.6910
## s.e. 0.0772 0.0418
##
## sigma^2 estimated as 0.01435: log likelihood=149.29
## AIC=-292.59 AICc=-292.47 BIC=-282.5qauselec_arima0.1.1.0.1.2.4 <- Arima(
qauselec, lambda = lambda_qauselec,
order = c(0, 1, 1), seasonal = c(0, 1, 2)
)
qauselec_arima0.1.1.0.1.2.4## Series: qauselec
## ARIMA(0,1,1)(0,1,2)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ma1 sma1 sma2
## -0.5037 -0.7993 0.1249
## s.e. 0.0710 0.0870 0.0864
##
## sigma^2 estimated as 0.01425: log likelihood=150.36
## AIC=-292.73 AICc=-292.53 BIC=-279.28qauselec_autoarima <- auto.arima(
qauselec, lambda = lambda_qauselec
)
qauselec_autoarima## Series: qauselec
## ARIMA(1,1,1)(1,1,2)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ar1 ma1 sar1 sma1 sma2
## 0.2523 -0.6905 0.8878 -1.6954 0.7641
## s.e. 0.1562 0.1279 0.0864 0.0983 0.0772
##
## sigma^2 estimated as 0.01345: log likelihood=156.42
## AIC=-300.84 AICc=-300.43 BIC=-280.67# According to AIC values, ARIMA(1, 1, 1)(1, 1, 2)[4] with Box-Cox transformation is the best model# d
checkresiduals(qauselec_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,1,2)[4]
## Q* = 10.605, df = 3, p-value = 0.01407
##
## Model df: 5. Total lags used: 8# The residuals aren't like white noise.
# try using other models.
checkresiduals(qauselec_arima0.0.1.0.1.1.4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,1)(0,1,1)[4]
## Q* = 84.595, df = 6, p-value = 4.441e-16
##
## Model df: 2. Total lags used: 8checkresiduals(qauselec_arima0.1.1.0.1.1.4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[4]
## Q* = 15.102, df = 6, p-value = 0.01948
##
## Model df: 2. Total lags used: 8checkresiduals(qauselec_arima0.1.1.0.1.2.4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,2)[4]
## Q* = 13.684, df = 5, p-value = 0.01775
##
## Model df: 3. Total lags used: 8# The residuals don't resemble white noise so I'm going to use the best model.e
fc_qauselec_autoarima <- forecast(
qauselec_autoarima, h = 8
)
autoplot(fc_qauselec_autoarima)
# The forecasts are reasonable.
# f. Compare the forecasts obtained using ets().
fc_qauselec_ets <- forecast(
ets(qauselec), h = 8
)
autoplot(fc_qauselec_ets)
# These forecasts also are reasonable.
#17 
## [1] 525.8100 513.8023 499.6705## Time Series:
## Start = 1969
## End = 1971
## Frequency = 1
## [1] 526.2057 514.0658 500.0111## [1] 526.2057 514.0658 500.0111