Session 3 Script 5
Mary Keller
6/19/2021
library(quantmod)## Warning: package 'quantmod' 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##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Warning: package 'TTR' was built under R version 4.0.5## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoooptions('getSymbols.warning4.0'=FALSE)
library(tsbox)## Warning: package 'tsbox' was built under R version 4.0.5library(fpp3)## Warning: package 'fpp3' was built under R version 4.0.5## -- Attaching packages -------------------------------------------- fpp3 0.4.0 --## v tibble 3.0.4 v tsibble 1.0.1
## v dplyr 1.0.2 v tsibbledata 0.3.0
## v tidyr 1.1.2 v feasts 0.2.1
## v lubridate 1.7.9.2 v fable 0.3.0
## v ggplot2 3.3.3## Warning: package 'tsibble' was built under R version 4.0.5## Warning: package 'tsibbledata' was built under R version 4.0.5## Warning: package 'feasts' was built under R version 4.0.5## Warning: package 'fabletools' was built under R version 4.0.5## Warning: package 'fable' was built under R version 4.0.5## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks xts::first()
## x tsibble::index() masks zoo::index()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks xts::last()
## x tsibble::setdiff() masks base::setdiff()
## x tsibble::union() masks base::union()library(stringr)
library(tictoc)## Warning: package 'tictoc' was built under R version 4.0.5Quarterly US pork production data
Pork <- getSymbols('IPN311611T3PNQ',src='FRED',auto.assign = FALSE)
autoplot(Pork)
Data is quarter and we foucs on 1990 and 2019.
Pork <- Pork %>%
ts_tsibble() %>%
mutate(time = yearquarter(time)) %>%
filter_index('1990-1' ~ '2019-4')
Pork %>%
gg_tsdisplay(value)
myl <- Pork %>% features(value,guerrero)
myl## # A tibble: 1 x 1
## lambda_guerrero
## <dbl>
## 1 0.234Pork %>%
mutate(bcv = box_cox(value,myl)) %>%
gg_tsdisplay(bcv)
Guerrero : 0.2339248
Look at box cox vs log. They look the same.
Pork %>%
mutate(logv = log(value)) %>%
gg_tsdisplay(logv)
Pork %>%
mutate(bcv = box_cox(value,myl)) %>%
features(bcv,unitroot_nsdiffs)## # A tibble: 1 x 1
## nsdiffs
## <int>
## 1 1Since nsdiffs=1, take seasonal difference before checking ndiffs. ,4 for quarterly data.
Pork %>%
mutate(d4bcv = difference(box_cox(value,myl),4)) %>%
features(d4bcv,unitroot_ndiffs)## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 0Pork %>%
mutate(d4bcv = difference(box_cox(value,myl),4)) %>%
gg_tsdisplay(d4bcv,plot_type = 'partial')## Warning: Removed 4 row(s) containing missing values (geom_path).## Warning: Removed 4 rows containing missing values (geom_point).
The orders are hard to sort from the ACF and PACF
Lag 3 & 5 could be an interaction of seasonal and non-seasonal
Perhaps MA(2) and SMA(1)?
tic()
fit <- Pork %>% model(
arima002011 = ARIMA(box_cox(value,myl)~pdq(0,0,2)+PDQ(0,1,1)),
step = ARIMA(box_cox(value,myl)),
search = ARIMA(box_cox(value,myl),stepwise = FALSE,approximation = FALSE)
)Fit this mode : arima002011
toc()## 14.24 sec elapsedfit %>%
pivot_longer(everything(),names_to="Model Name",values_to="Orders")## # A mable: 3 x 2
## # Key: Model Name [3]
## `Model Name` Orders
## <chr> <model>
## 1 arima002011 <ARIMA(0,0,2)(0,1,1)[4]>
## 2 step <ARIMA(1,0,2)(2,1,1)[4] w/ drift>
## 3 search <ARIMA(3,0,0)(0,1,1)[4] w/ drift>glance(fit) %>%
arrange(AICc) %>%
select(.model:BIC)## # A tibble: 3 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 search 0.00426 147. -282. -281. -266.
## 2 step 0.00439 146. -277. -275. -255.
## 3 arima002011 0.00955 104. -201. -200. -190.Looks like AICC -281 stepwise is much better and ok with search.
fit %>%
select(arima002011) %>%
report()## Series: value
## Model: ARIMA(0,0,2)(0,1,1)[4]
## Transformation: box_cox(value, myl)
##
## Coefficients:
## ma1 ma2 sma1
## 0.9116 0.4303 -0.3029
## s.e. 0.1037 0.0756 0.0935
##
## sigma^2 estimated as 0.009547: log likelihood=104.36
## AIC=-200.71 AICc=-200.35 BIC=-189.77fit %>%
select(arima002011) %>%
gg_tsresiduals()
fit %>%
select(arima002011) %>%
augment() %>%
features(.innov,ljung_box,lag=20,dof=3)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arima002011 52.8 0.0000151fit %>%
select(step) %>%
report()## Series: value
## Model: ARIMA(1,0,2)(2,1,1)[4] w/ drift
## Transformation: box_cox(value, myl)
##
## Coefficients:
## ar1 ma1 ma2 sar1 sar2 sma1 constant
## 0.8056 0.1338 0.1993 0.0288 -0.0413 -0.9492 0.0112
## s.e. 0.0794 0.1088 0.1054 0.1168 0.1112 0.0964 0.0009
##
## sigma^2 estimated as 0.004388: log likelihood=146.4
## AIC=-276.81 AICc=-275.44 BIC=-254.92fit %>%
select(step) %>%
gg_tsresiduals()
fit %>%
select(step) %>%
augment() %>%
features(.innov,ljung_box,lag=20,dof=7)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 step 18.1 0.154fit %>%
select(search) %>%
report()## Series: value
## Model: ARIMA(3,0,0)(0,1,1)[4] w/ drift
## Transformation: box_cox(value, myl)
##
## Coefficients:
## ar1 ar2 ar3 sma1 constant
## 0.9592 0.0771 -0.206 -0.9477 0.0096
## s.e. 0.0918 0.1299 0.094 0.0791 0.0007
##
## sigma^2 estimated as 0.004259: log likelihood=146.99
## AIC=-281.98 AICc=-281.2 BIC=-265.57Decide with one with 5 paraments instead of 7 so with step.
fit %>%
select(step) %>%
gg_tsresiduals()
fit %>%
select(step) %>%
augment() %>%
features(.innov,ljung_box,lag=20,dof=5)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 step 18.1 0.258.60 for .init = 60
tic()
acc <- Pork %>%
stretch_tsibble(.init=60) %>%
model(arima102211=ARIMA(box_cox(value,myl) ~ 1+pdq(1,0,2)+PDQ(2,1,1)),
arima300011=ARIMA(box_cox(value,myl) ~ 1+pdq(3,0,0)+PDQ(0,1,1))
) %>%
forecast(h=12) %>%
group_by(.model,.id) %>%
mutate(h=row_number()) %>%
ungroup() %>%
accuracy(Pork, by = c('h','.model'), measures=list(rmse=RMSE))## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced## Warning in sqrt(diag(best$var.coef)): NaNs produced## Warning in log(s2): NaNs produced## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 12 observations are missing between 2019 Q3 and 2022 Q2toc()## 29.7 sec elapsedacc %>% print(n=Inf)## # A tibble: 24 x 4
## h .model .type rmse
## <int> <chr> <chr> <dbl>
## 1 1 arima102211 Test 2.29
## 2 1 arima300011 Test 2.24
## 3 2 arima102211 Test 3.38
## 4 2 arima300011 Test 3.38
## 5 3 arima102211 Test 4.16
## 6 3 arima300011 Test 4.20
## 7 4 arima102211 Test 4.65
## 8 4 arima300011 Test 4.72
## 9 5 arima102211 Test 5.12
## 10 5 arima300011 Test 5.20
## 11 6 arima102211 Test 5.46
## 12 6 arima300011 Test 5.53
## 13 7 arima102211 Test 5.74
## 14 7 arima300011 Test 5.80
## 15 8 arima102211 Test 5.92
## 16 8 arima300011 Test 5.96
## 17 9 arima102211 Test 6.06
## 18 9 arima300011 Test 6.08
## 19 10 arima102211 Test 6.15
## 20 10 arima300011 Test 6.17
## 21 11 arima102211 Test 6.28
## 22 11 arima300011 Test 6.28
## 23 12 arima102211 Test 6.46
## 24 12 arima300011 Test 6.42better accuaracy with arima300011.
If have to pick one, they are pretty close
acc %>%
ggplot(aes(x=h,y=rmse,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:12, minor_breaks = NULL)+
labs(y='RMSE',title='Pork: Point Prediction Accuracy')
pretty close call: search model is better at h=1, step model better at h>2
now STL+ARIMA
Pork %>%
model(STL(value)) %>%
components() %>%
autoplot()
Add myl transformation decomosition
Pork %>%
model(STL(box_cox(value,myl))) %>%
components() %>%
autoplot()
tic()
fit <- Pork %>%
model(
`stl+arima1` = decomposition_model(
STL(box_cox(value,myl)),
ARIMA(season_adjust~PDQ(0,0,0))),
`stl+arima2` = decomposition_model(
STL(box_cox(value,myl)),
ARIMA(season_adjust~PDQ(0,0,0),stepwise=FALSE,approximation = FALSE))
)Search for non seasons arina models for decomposition.
toc()## 1.11 sec elapsedfit %>%
select('stl+arima1') %>%
report()## Series: value
## Model: STL decomposition model
## Transformation: box_cox(value, myl)
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ARIMA(0,1,0) w/ drift
##
## Coefficients:
## constant
## 0.0155
## s.e. 0.0058
##
## sigma^2 estimated as 0.003967: log likelihood=157.98
## AIC=-311.97 AICc=-311.86 BIC=-306.44
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 0fit %>%
select('stl+arima1') %>%
gg_tsresiduals()## Warning: Removed 4 row(s) containing missing values (geom_path).## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 4 rows containing non-finite values (stat_bin).
fit %>%
select('stl+arima1') %>%
augment() %>%
features(.innov,ljung_box,lag=20,dof=0)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 stl+arima1 33.4 0.0308Pvalue is .03 which is lease than .05. fails residual test there is a white noise.
fit %>%
select('stl+arima2') %>%
report()## Series: value
## Model: STL decomposition model
## Transformation: box_cox(value, myl)
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ARIMA(0,1,4) w/ drift
##
## Coefficients:
## ma1 ma2 ma3 ma4 constant
## 0.0292 0.1997 0.0084 -0.3285 0.0156
## s.e. 0.0939 0.1035 0.0864 0.1699 0.0051
##
## sigma^2 estimated as 0.003741: log likelihood=163.15
## AIC=-314.3 AICc=-313.54 BIC=-297.73
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 0fit %>%
select('stl+arima2') %>%
gg_tsresiduals()## Warning: Removed 4 row(s) containing missing values (geom_path).## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 4 rows containing non-finite values (stat_bin).
fit %>%
select('stl+arima2') %>%
augment() %>%
features(.innov,ljung_box,lag=20,dof=5)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 stl+arima2 26.4 0.0338Pvalue is .03 which is lease than .05. fails residual test there is a white noise.
tic()
acc <- Pork %>%
stretch_tsibble(.init=60) %>%
model(arima300011=ARIMA(box_cox(value,myl) ~ 1+pdq(3,0,0)+PDQ(0,1,1)),
`stl+arima410`=decomposition_model(STL(box_cox(value,myl)),
ARIMA(season_adjust~1+pdq(4,1,0)+PDQ(0,0,0)))
) %>%
forecast(h=12) %>%
group_by(.model,.id) %>%
mutate(h=row_number()) %>%
ungroup() %>%
accuracy(Pork, by = c('h','.model'), measures=list(rmse=RMSE))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 12 observations are missing between 2019 Q3 and 2022 Q2toc()## 32.36 sec elapsedacc %>% print(n=Inf)## # A tibble: 24 x 4
## h .model .type rmse
## <int> <chr> <chr> <dbl>
## 1 1 arima300011 Test 2.24
## 2 1 stl+arima410 Test 2.29
## 3 2 arima300011 Test 3.38
## 4 2 stl+arima410 Test 3.48
## 5 3 arima300011 Test 4.20
## 6 3 stl+arima410 Test 4.30
## 7 4 arima300011 Test 4.72
## 8 4 stl+arima410 Test 4.78
## 9 5 arima300011 Test 5.20
## 10 5 stl+arima410 Test 5.15
## 11 6 arima300011 Test 5.53
## 12 6 stl+arima410 Test 5.47
## 13 7 arima300011 Test 5.80
## 14 7 stl+arima410 Test 5.86
## 15 8 arima300011 Test 5.96
## 16 8 stl+arima410 Test 6.20
## 17 9 arima300011 Test 6.08
## 18 9 stl+arima410 Test 6.40
## 19 10 arima300011 Test 6.17
## 20 10 stl+arima410 Test 6.51
## 21 11 arima300011 Test 6.28
## 22 11 stl+arima410 Test 6.77
## 23 12 arima300011 Test 6.42
## 24 12 stl+arima410 Test 7.20acc %>%
ggplot(aes(x=h,y=rmse,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:12, minor_breaks = NULL)+
labs(y='RMSE',title='Pork: Point Prediction Accuracy')
Looks like arima 300011 towards the end however, they switch /cross in 2 year since everything is in quarter.
Two ETS model with original and transform series. Turn seasonally off for seasonally adjust ed data.
Now ETS and STL+ETS:
tic()
fit <- Pork %>%
model(ets1 = ETS(value),
ets2 = ETS(box_cox(value,myl)),
`stl+ets` = decomposition_model(STL(box_cox(value,myl)),
ETS(season_adjust~season('N')))
)toc()## 1.55 sec elapsedfit %>%
select(ets1) %>%
report()## Series: value
## Model: ETS(M,Ad,M)
## Smoothing parameters:
## alpha = 0.9990932
## beta = 0.002378414
## gamma = 0.000103457
## phi = 0.9799766
##
## Initial states:
## l b s1 s2 s3 s4
## 55.11466 0.7307042 1.05932 0.9600011 0.9612917 1.019388
##
## sigma^2: 6e-04
##
## AIC AICc BIC
## 732.4824 734.5385 760.1893fit %>%
select(ets1) %>%
gg_tsresiduals()
fit %>% select(ets1) %>% augment() %>%
features(.innov,ljung_box,lag=20,dof=10)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ets1 29.4 0.00108Pvalue is 0.001 which is less than .05
fit %>% select(ets2) %>% report()## Series: value
## Model: ETS(A,A,A)
## Transformation: box_cox(value, myl)
## Smoothing parameters:
## alpha = 0.9998999
## beta = 0.0001000032
## gamma = 0.0001000173
##
## Initial states:
## l b s1 s2 s3 s4
## 6.73242 0.01674505 0.1686547 -0.1127978 -0.1113694 0.05551253
##
## sigma^2: 0.0046
##
## AIC AICc BIC
## -61.26461 -59.59794 -36.32845fit %>% select(ets2) %>% gg_tsresiduals()
fit %>% select(ets2) %>% augment() %>%
features(.innov,ljung_box,lag=20,dof=9)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ets2 27.8 0.00348fit %>% select('stl+ets') %>% report()## Series: value
## Model: STL decomposition model
## Transformation: box_cox(value, myl)
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ETS(A,A,N)
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.0001000049
##
## Initial states:
## l b
## 6.730403 0.01557049
##
## sigma^2: 0.004
##
## AIC AICc BIC
## -81.58664 -81.05093 -67.73322
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 0fit %>% select('stl+ets') %>% gg_tsresiduals()## Warning: Removed 4 row(s) containing missing values (geom_path).## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 4 rows containing non-finite values (stat_bin).
fit %>% select('stl+ets') %>% augment() %>%
features(.innov,ljung_box,lag=20,dof=2)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 stl+ets 33.4 0.0150Comment: None of the ETS models have white noise residuals, but thatis not usually a critical issue for those models.
Look at all the models.
make .init = 60
tic()
acc <- Pork %>%
stretch_tsibble(.init=60) %>%
model(`arima300011`=ARIMA(box_cox(value,myl) ~ 1+pdq(3,0,0)+PDQ(0,1,1)),
`stl+arima410`=decomposition_model(STL(box_cox(value,myl)),
ARIMA(season_adjust~1+pdq(4,1,0)+PDQ(0,0,0))),
`ets(MAdM)`=ETS(value~error('M')+trend('Ad')+season('M')),
`ets(bcAAA)`=ETS(box_cox(value,myl)~error('A')+trend('A')+season('A')),
`stl+ets(AAN)`=decomposition_model(STL(box_cox(value,myl)),
ETS(season_adjust~error('A')+trend('A')+season('N')))
) %>%
forecast(h=12) %>%
group_by(.model,.id) %>%
mutate(h=row_number()) %>%
ungroup() %>%
accuracy(Pork, by = c('h','.model'), measures=list(rmse=RMSE))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 12 observations are missing between 2019 Q3 and 2022 Q2toc()## 78.06 sec elapsedacc %>% print(n=Inf)## # A tibble: 60 x 4
## h .model .type rmse
## <int> <chr> <chr> <dbl>
## 1 1 arima300011 Test 2.24
## 2 1 ets(bcAAA) Test 2.17
## 3 1 ets(MAdM) Test 2.22
## 4 1 stl+arima410 Test 2.29
## 5 1 stl+ets(AAN) Test 2.21
## 6 2 arima300011 Test 3.38
## 7 2 ets(bcAAA) Test 3.33
## 8 2 ets(MAdM) Test 3.41
## 9 2 stl+arima410 Test 3.48
## 10 2 stl+ets(AAN) Test 3.39
## 11 3 arima300011 Test 4.20
## 12 3 ets(bcAAA) Test 4.16
## 13 3 ets(MAdM) Test 4.32
## 14 3 stl+arima410 Test 4.30
## 15 3 stl+ets(AAN) Test 4.18
## 16 4 arima300011 Test 4.72
## 17 4 ets(bcAAA) Test 4.65
## 18 4 ets(MAdM) Test 4.94
## 19 4 stl+arima410 Test 4.78
## 20 4 stl+ets(AAN) Test 4.63
## 21 5 arima300011 Test 5.20
## 22 5 ets(bcAAA) Test 5.10
## 23 5 ets(MAdM) Test 5.53
## 24 5 stl+arima410 Test 5.15
## 25 5 stl+ets(AAN) Test 5.08
## 26 6 arima300011 Test 5.53
## 27 6 ets(bcAAA) Test 5.47
## 28 6 ets(MAdM) Test 6.03
## 29 6 stl+arima410 Test 5.47
## 30 6 stl+ets(AAN) Test 5.45
## 31 7 arima300011 Test 5.80
## 32 7 ets(bcAAA) Test 5.83
## 33 7 ets(MAdM) Test 6.50
## 34 7 stl+arima410 Test 5.86
## 35 7 stl+ets(AAN) Test 5.79
## 36 8 arima300011 Test 5.96
## 37 8 ets(bcAAA) Test 6.11
## 38 8 ets(MAdM) Test 6.91
## 39 8 stl+arima410 Test 6.20
## 40 8 stl+ets(AAN) Test 6.05
## 41 9 arima300011 Test 6.08
## 42 9 ets(bcAAA) Test 6.29
## 43 9 ets(MAdM) Test 7.19
## 44 9 stl+arima410 Test 6.40
## 45 9 stl+ets(AAN) Test 6.20
## 46 10 arima300011 Test 6.17
## 47 10 ets(bcAAA) Test 6.52
## 48 10 ets(MAdM) Test 7.39
## 49 10 stl+arima410 Test 6.51
## 50 10 stl+ets(AAN) Test 6.29
## 51 11 arima300011 Test 6.28
## 52 11 ets(bcAAA) Test 6.89
## 53 11 ets(MAdM) Test 7.72
## 54 11 stl+arima410 Test 6.77
## 55 11 stl+ets(AAN) Test 6.57
## 56 12 arima300011 Test 6.42
## 57 12 ets(bcAAA) Test 7.35
## 58 12 ets(MAdM) Test 8.22
## 59 12 stl+arima410 Test 7.20
## 60 12 stl+ets(AAN) Test 7.03Graph it to see which one is the best.
acc %>%
ggplot(aes(x=h,y=rmse,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:12, minor_breaks = NULL)+
labs(y='RMSE',title='Pork: Point Prediction Accuracy')
# which model gives best forecast of mean? depends on h….
looks like ARIMA300011 was the best one. closes to h.
Electronics and Applicance store retail sales (monthly)
RSEAS <- getSymbols('MRTSSM443USN',src='FRED',auto.assign = FALSE)
RSEAS <- RSEAS %>%
ts_tsibble() %>%
mutate(time=yearmonth(time)) %>%
filter(year(time) <= 2019)
RSEAS %>% gg_tsdisplay(value)
myl <- RSEAS %>% features(value,guerrero)
myl## # A tibble: 1 x 1
## lambda_guerrero
## <dbl>
## 1 0.296Guerrero : 0.2964748
RSEAS %>% mutate(bcv=box_cox(value,myl)) %>% gg_tsdisplay(bcv)
RSEAS %>% mutate(logv=log(value)) %>% gg_tsdisplay(logv)
Let’s try some seasonal ARIMA models:
fit <- RSEAS %>%
model(arima1=ARIMA(box_cox(value,myl)),
arima2=ARIMA(log(value))
)
fit %>%
pivot_longer(everything(),names_to="Model Name",values_to="Orders")## # A mable: 2 x 2
## # Key: Model Name [2]
## `Model Name` Orders
## <chr> <model>
## 1 arima1 <ARIMA(0,1,1)(1,1,1)[12]>
## 2 arima2 <ARIMA(0,1,1)(1,1,1)[12]>glance(fit) %>%
select(.model:BIC)## # A tibble: 2 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 arima1 0.149 -150. 308. 308. 323.
## 2 arima2 0.000707 714. -1420. -1420. -1405.don’t even think about comparing AICc here – different scales!
fit %>%
select(arima1) %>%
report()## Series: value
## Model: ARIMA(0,1,1)(1,1,1)[12]
## Transformation: box_cox(value, myl)
##
## Coefficients:
## ma1 sar1 sma1
## -0.3811 0.1573 -0.5131
## s.e. 0.0543 0.1239 0.1059
##
## sigma^2 estimated as 0.1488: log likelihood=-150.18
## AIC=308.36 AICc=308.48 BIC=323.47fit %>%
select(arima1) %>%
gg_tsresiduals(lag_max=24)
fit %>%
select(arima1) %>%
augment() %>%
features(.innov,ljung_box,lag=24,dof=3)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arima1 54.0 0.0000996fit %>%
select(arima2) %>%
report()## Series: value
## Model: ARIMA(0,1,1)(1,1,1)[12]
## Transformation: log(value)
##
## Coefficients:
## ma1 sar1 sma1
## -0.3832 0.1831 -0.5722
## s.e. 0.0532 0.1124 0.0920
##
## sigma^2 estimated as 0.0007067: log likelihood=714.14
## AIC=-1420.28 AICc=-1420.15 BIC=-1405.17fit %>%
select(arima2) %>%
gg_tsresiduals(lag_max=24)
fit %>%
select(arima2) %>%
augment() %>%
features(.innov,ljung_box,lag=24,dof=3)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arima2 54.7 0.0000774Both stepwise models fail residual tests – not white noise
tic()
fit <- RSEAS %>%
model(arima1=ARIMA(box_cox(value,myl),stepwise = FALSE,approximation = FALSE),
arima2=ARIMA(log(value),stepwise = FALSE,approximation = FALSE)
)toc()## 136.37 sec elapsedfit %>%
pivot_longer(everything(),names_to="Model Name",values_to="Orders")## # A mable: 2 x 2
## # Key: Model Name [2]
## `Model Name` Orders
## <chr> <model>
## 1 arima1 <ARIMA(4,1,0)(0,1,1)[12]>
## 2 arima2 <ARIMA(4,1,0)(0,1,1)[12]>glance(fit) %>%
select(.model:BIC)## # A tibble: 2 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 arima1 0.147 -147. 306. 307. 329.
## 2 arima2 0.000698 717. -1422. -1422. -1400.don’t even think about comparing AICc here – different scales!
fit %>%
select(arima1) %>%
report()## Series: value
## Model: ARIMA(4,1,0)(0,1,1)[12]
## Transformation: box_cox(value, myl)
##
## Coefficients:
## ar1 ar2 ar3 ar4 sma1
## -0.3568 -0.1828 -0.0189 -0.1384 -0.3806
## s.e. 0.0558 0.0589 0.0586 0.0551 0.0556
##
## sigma^2 estimated as 0.147: log likelihood=-147.12
## AIC=306.24 AICc=306.5 BIC=328.9fit %>%
select(arima1) %>%
gg_tsresiduals(lag_max=24)
fit %>%
select(arima1) %>%
augment() %>%
features(.innov,ljung_box,lag=24,dof=5)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arima1 37.4 0.00699fit %>%
select(arima2) %>%
report()## Series: value
## Model: ARIMA(4,1,0)(0,1,1)[12]
## Transformation: log(value)
##
## Coefficients:
## ar1 ar2 ar3 ar4 sma1
## -0.3656 -0.1764 -0.0079 -0.1407 -0.4284
## s.e. 0.0557 0.0591 0.0589 0.0551 0.0557
##
## sigma^2 estimated as 0.0006983: log likelihood=717.1
## AIC=-1422.2 AICc=-1421.93 BIC=-1399.53fit %>%
select(arima2) %>%
gg_tsresiduals(lag_max=24)
fit %>%
select(arima2) %>%
augment() %>%
features(.innov,ljung_box,lag=24,dof=5)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arima2 37.9 0.00620Both models fail ljung-box test, but the signficant lags are at lag 10 or higher
and not at lags 12 or 24. Not bad….
fc <- fit %>%
forecast(h="2 years")
fc %>%
autoplot() +
facet_grid(.model ~ .)
fc %>%
autoplot(RSEAS) +
facet_grid(.model ~ .)
fc %>%
autoplot(filter_index(RSEAS,"2017-01"~.)) +
facet_grid(.model ~ .)
hilo(fc) %>%
arrange(time) %>%
print(n=Inf)## # A tsibble: 48 x 6 [1M]
## # Key: .model [2]
## .model time value .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 arima1 2020 Jan t(N(43, 0.15)) 6947. [6700.475, 7196.055]80
## 2 arima2 2020 Jan t(N(8.8, 7e-04)) 6961. [6726.440, 7197.794]80
## 3 arima1 2020 Feb t(N(42, 0.21)) 6481. [6202.455, 6763.597]80
## 4 arima2 2020 Feb t(N(8.8, 0.00098)) 6518. [6258.667, 6781.338]80
## 5 arima1 2020 Mar t(N(43, 0.26)) 7093. [6761.392, 7428.306]80
## 6 arima2 2020 Mar t(N(8.9, 0.0012)) 7114. [6797.676, 7435.322]80
## 7 arima1 2020 Apr t(N(42, 0.32)) 6372. [6029.517, 6718.992]80
## 8 arima2 2020 Apr t(N(8.8, 0.0015)) 6398. [6079.930, 6721.399]80
## 9 arima1 2020 May t(N(43, 0.36)) 6948. [6563.147, 7338.042]80
## 10 arima2 2020 May t(N(8.8, 0.0017)) 6962. [6596.798, 7334.636]80
## 11 arima1 2020 Jun t(N(43, 0.41)) 6910. [6500.100, 7326.667]80
## 12 arima2 2020 Jun t(N(8.8, 0.002)) 6935. [6545.892, 7333.613]80
## 13 arima1 2020 Jul t(N(43, 0.46)) 6983. [6544.849, 7429.311]80
## 14 arima2 2020 Jul t(N(8.9, 0.0022)) 6999. [6581.559, 7425.958]80
## 15 arima1 2020 Aug t(N(44, 0.51)) 7362. [6883.470, 7848.363]80
## 16 arima2 2020 Aug t(N(8.9, 0.0024)) 7375. [6913.025, 7848.094]80
## 17 arima1 2020 Sep t(N(43, 0.57)) 6957. [6474.140, 7448.335]80
## 18 arima2 2020 Sep t(N(8.8, 0.0027)) 6978. [6518.824, 7448.812]80
## 19 arima1 2020 Oct t(N(43, 0.62)) 6928. [6425.868, 7439.342]80
## 20 arima2 2020 Oct t(N(8.8, 0.0029)) 6944. [6467.763, 7433.417]80
## 21 arima1 2020 Nov t(N(47, 0.67)) 9080. [8448.586, 9723.918]80
## 22 arima2 2020 Nov t(N(9.1, 0.0032)) 9083. [8435.657, 9749.664]80
## 23 arima1 2020 Dec t(N(49, 0.72)) 10685. [9949.751, 11433.953]80
## 24 arima2 2020 Dec t(N(9.3, 0.0034)) 10686. [9895.583, 11500.170]80
## 25 arima1 2021 Jan t(N(43, 0.93)) 6804. [6196.310, 7426.347]80
## 26 arima2 2021 Jan t(N(8.8, 0.0044)) 6833. [6264.008, 7421.588]80
## 27 arima1 2021 Feb t(N(42, 1.1)) 6348. [5727.468, 6985.639]80
## 28 arima2 2021 Feb t(N(8.8, 0.005)) 6402. [5832.264, 6992.951]80
## 29 arima1 2021 Mar t(N(43, 1.2)) 6946. [6245.395, 7666.469]80
## 30 arima2 2021 Mar t(N(8.8, 0.0056)) 6984. [6326.328, 7666.343]80
## 31 arima1 2021 Apr t(N(42, 1.4)) 6240. [5553.853, 6948.191]80
## 32 arima2 2021 Apr t(N(8.7, 0.0063)) 6284. [5659.713, 6933.872]80
## 33 arima1 2021 May t(N(43, 1.5)) 6808. [6046.898, 7593.666]80
## 34 arima2 2021 May t(N(8.8, 0.0068)) 6840. [6131.939, 7577.352]80
## 35 arima1 2021 Jun t(N(43, 1.6)) 6772. [5979.777, 7589.614]80
## 36 arima2 2021 Jun t(N(8.8, 0.0074)) 6814. [6078.803, 7581.837]80
## 37 arima1 2021 Jul t(N(43, 1.7)) 6846. [6015.625, 7704.539]80
## 38 arima2 2021 Jul t(N(8.8, 0.008)) 6879. [6107.221, 7685.378]80
## 39 arima1 2021 Aug t(N(44, 1.9)) 7219. [6325.749, 8143.659]80
## 40 arima2 2021 Aug t(N(8.9, 0.0086)) 7249. [6407.319, 8130.356]80
## 41 arima1 2021 Sep t(N(43, 2)) 6821. [5933.178, 7741.874]80
## 42 arima2 2021 Sep t(N(8.8, 0.0092)) 6860. [6036.459, 7723.873]80
## 43 arima1 2021 Oct t(N(43, 2.1)) 6794. [5880.186, 7742.163]80
## 44 arima2 2021 Oct t(N(8.8, 0.0098)) 6828. [5983.126, 7716.116]80
## 45 arima1 2021 Nov t(N(47, 2.3)) 8917. [7775.949, 10100.001]80
## 46 arima2 2021 Nov t(N(9.1, 0.01)) 8933. [7795.347, 10130.503]80
## 47 arima1 2021 Dec t(N(49, 2.4)) 10502. [9184.090, 11867.515]80
## 48 arima2 2021 Dec t(N(9.3, 0.011)) 10511. [9135.560, 11961.462]80
## # ... with 1 more variable: `95%` <hilo>Cross validation: Does the type of transform matter?
tictoc::tic()
acc <- RSEAS %>%
stretch_tsibble(.init=225) %>%
model(arima1=ARIMA(box_cox(value,myl)~0+pdq(4,1,0)+PDQ(0,1,1)),
arima2=ARIMA(log(value)~0+pdq(4,1,0)+PDQ(0,1,1))
) %>%
forecast(h=24) %>%
group_by(.model,.id) %>%
mutate(h=row_number()) %>%
ungroup() %>%
accuracy(RSEAS, by = c('h','.model'), measures=list(rmse=RMSE))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 24 observations are missing between 2020 Jan and 2021 Dectictoc::toc()## 71.38 sec elapsedacc %>% print(n=Inf)## # A tibble: 48 x 4
## h .model .type rmse
## <int> <chr> <chr> <dbl>
## 1 1 arima1 Test 236.
## 2 1 arima2 Test 238.
## 3 2 arima1 Test 271.
## 4 2 arima2 Test 271.
## 5 3 arima1 Test 299.
## 6 3 arima2 Test 295.
## 7 4 arima1 Test 313.
## 8 4 arima2 Test 308.
## 9 5 arima1 Test 319.
## 10 5 arima2 Test 315.
## 11 6 arima1 Test 327.
## 12 6 arima2 Test 324.
## 13 7 arima1 Test 329.
## 14 7 arima2 Test 323.
## 15 8 arima1 Test 338.
## 16 8 arima2 Test 331.
## 17 9 arima1 Test 355.
## 18 9 arima2 Test 351.
## 19 10 arima1 Test 371.
## 20 10 arima2 Test 368.
## 21 11 arima1 Test 390.
## 22 11 arima2 Test 386.
## 23 12 arima1 Test 402.
## 24 12 arima2 Test 397.
## 25 13 arima1 Test 484.
## 26 13 arima2 Test 482.
## 27 14 arima1 Test 498.
## 28 14 arima2 Test 493.
## 29 15 arima1 Test 528.
## 30 15 arima2 Test 515.
## 31 16 arima1 Test 545.
## 32 16 arima2 Test 529.
## 33 17 arima1 Test 552.
## 34 17 arima2 Test 537.
## 35 18 arima1 Test 554.
## 36 18 arima2 Test 539.
## 37 19 arima1 Test 536.
## 38 19 arima2 Test 520.
## 39 20 arima1 Test 533.
## 40 20 arima2 Test 516.
## 41 21 arima1 Test 560.
## 42 21 arima2 Test 550.
## 43 22 arima1 Test 583.
## 44 22 arima2 Test 573.
## 45 23 arima1 Test 616.
## 46 23 arima2 Test 606.
## 47 24 arima1 Test 634.
## 48 24 arima2 Test 625.acc %>%
ggplot(aes(x=h,y=rmse,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:24, minor_breaks = NULL)+
labs(y='RMSE',title='RSEAS: Point Prediction Accuracy')
The log transform does slightly better almost across the board
Next, let’s look for ETS+ARIMA, ETS, and STL+ETS
tic()
fit <- RSEAS %>%
model(
arima=ARIMA(log(value)~0+pdq(4,1,0)+PDQ(0,1,1)),
`stl+arima`=decomposition_model(STL(log(value)),
ARIMA(season_adjust~PDQ(0,0,0),
stepwise = FALSE,approximation = FALSE)),
ets1=ETS(value),
ets2=ETS(log(value)),
`stl+ets`=decomposition_model(STL(log(value)),
ETS(season_adjust~season('N')))
)toc()
fit %>% select('stl+arima') %>% report()
fit %>% select('stl+arima') %>% gg_tsresiduals()
fit %>% select('stl+arima') %>% augment() %>%
features(.innov,ljung_box,lag=12,dof=5)
fit %>% select('stl+arima') %>% augment() %>%
features(.innov,ljung_box,lag=12,dof=5)fit %>% select('ets1') %>% report()## Series: value
## Model: ETS(M,A,M)
## Smoothing parameters:
## alpha = 0.5028789
## beta = 0.008467165
## gamma = 0.4704238
##
## Initial states:
## l b s1 s2 s3 s4 s5 s6
## 3950.324 51.6437 1.642281 1.070043 0.9564365 0.9257562 0.9457197 0.9685462
## s7 s8 s9 s10 s11 s12
## 0.9467981 0.8983164 0.8897136 0.9386781 0.8786271 0.939084
##
## sigma^2: 7e-04
##
## AIC AICc BIC
## 5505.508 5507.432 5570.398fit %>% select('ets1') %>% gg_tsresiduals()
fit %>% select('ets1') %>% augment() %>%
features(.innov,ljung_box,lag=24,dof=16)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ets1 46.7 0.000000174fit %>% select('ets2') %>% report()## Series: value
## Model: ETS(M,Ad,M)
## Transformation: log(value)
## Smoothing parameters:
## alpha = 0.5245922
## beta = 0.0001073937
## gamma = 0.4066682
## phi = 0.9788524
##
## Initial states:
## l b s1 s2 s3 s4 s5 s6
## 8.239854 0.01440058 1.059938 1.010434 0.9943403 0.9941961 0.9965518 0.9965459
## s7 s8 s9 s10 s11 s12
## 0.9938685 0.9896957 0.988435 0.9952858 0.9868694 0.9938393
##
## sigma^2: 0
##
## AIC AICc BIC
## -479.2819 -477.1242 -410.5739fit %>% select('ets2') %>% gg_tsresiduals()
fit %>% select('ets2') %>% augment() %>%
features(.innov,ljung_box,lag=24,dof=17)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 ets2 50.0 0.0000000143fit %>% select('stl+ets') %>% report()## Series: value
## Model: STL decomposition model
## Transformation: log(value)
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ETS(A,Ad,N)
## Smoothing parameters:
## alpha = 0.5962038
## beta = 0.0008418349
## phi = 0.9799997
##
## Initial states:
## l b
## 8.2294 0.01459199
##
## sigma^2: 5e-04
##
## AIC AICc BIC
## -615.4509 -615.1955 -592.5482
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 0fit %>% select('stl+ets') %>% gg_tsresiduals()## Warning: Removed 12 row(s) containing missing values (geom_path).## Warning: Removed 12 rows containing missing values (geom_point).## Warning: Removed 12 rows containing non-finite values (stat_bin).
fit %>% select('stl+ets') %>% augment() %>%
features(.innov,ljung_box,lag=24,dof=5)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 stl+ets 56.7 0.0000129fit %>% select('stl+ets') %>% augment() %>%
features(.innov,ljung_box,lag=12,dof=5)## # A tibble: 1 x 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 stl+ets 19.4 0.00707fc <- fit %>% forecast(h="2 years")
fc %>% autoplot() + facet_grid(.model ~ .)
fc %>% autoplot(RSEAS) + facet_grid(.model ~ .)
fc %>% autoplot(filter_index(RSEAS,"2017-01"~.)) + facet_grid(.model ~ .)
hilo(fc) %>% arrange(time) %>% print(n=Inf)## # A tsibble: 120 x 6 [1M]
## # Key: .model [5]
## .model time value .mean `80%`
## <chr> <mth> <dist> <dbl> <hilo>
## 1 arima 2020 Jan t(N(8.8, 7e-04)) 6961. [ 6726.440, 7197.794]80
## 2 stl+a~ 2020 Jan t(N(8.8, 0.00048)) 6824. [ 6634.293, 7016.367]80
## 3 ets1 2020 Jan N(7025, 33782) 7025. [ 6789.613, 7260.706]80
## 4 ets2 2020 Jan t(N(8.9, 0.00067)) 7028. [ 6795.818, 7262.224]80
## 5 stl+e~ 2020 Jan t(N(8.8, 0.00048)) 6840. [ 6648.429, 7033.274]80
## 6 arima 2020 Feb t(N(8.8, 0.00098)) 6518. [ 6258.667, 6781.338]80
## 7 stl+a~ 2020 Feb t(N(8.8, 0.00067)) 6684. [ 6463.236, 6907.231]80
## 8 ets1 2020 Feb N(6589, 37506) 6589. [ 6340.674, 6837.059]80
## 9 ets2 2020 Feb t(N(8.8, 0.00084)) 6650. [ 6404.391, 6899.696]80
## 10 stl+e~ 2020 Feb t(N(8.8, 0.00065)) 6723. [ 6504.836, 6943.784]80
## 11 arima 2020 Mar t(N(8.9, 0.0012)) 7114. [ 6797.676, 7435.322]80
## 12 stl+a~ 2020 Mar t(N(8.9, 0.00087)) 7007. [ 6744.482, 7272.925]80
## 13 ets1 2020 Mar N(7265, 55405) 7265. [ 6963.202, 7566.514]80
## 14 ets2 2020 Mar t(N(8.9, 0.0011)) 7335. [ 7033.098, 7642.202]80
## 15 stl+e~ 2020 Mar t(N(8.9, 0.00082)) 7075. [ 6817.854, 7335.946]80
## 16 arima 2020 Apr t(N(8.8, 0.0015)) 6398. [ 6079.930, 6721.399]80
## 17 stl+a~ 2020 Apr t(N(8.7, 0.0011)) 6296. [ 6026.941, 6569.209]80
## 18 ets1 2020 Apr N(6539, 53110) 6539. [ 6243.566, 6834.248]80
## 19 ets2 2020 Apr t(N(8.8, 0.0012)) 6600. [ 6309.171, 6896.888]80
## 20 stl+e~ 2020 Apr t(N(8.8, 0.00098)) 6319. [ 6066.748, 6574.774]80
## 21 arima 2020 May t(N(8.8, 0.0017)) 6962. [ 6596.798, 7334.636]80
## 22 stl+a~ 2020 May t(N(8.8, 0.0013)) 6728. [ 6421.562, 7038.990]80
## 23 ets1 2020 May N(7073, 72090) 7073. [ 6728.684, 7416.869]80
## 24 ets2 2020 May t(N(8.9, 0.0014)) 7131. [ 6791.115, 7478.304]80
## 25 stl+e~ 2020 May t(N(8.8, 0.0012)) 6806. [ 6512.017, 7104.067]80
## 26 arima 2020 Jun t(N(8.8, 0.002)) 6935. [ 6545.892, 7333.613]80
## 27 stl+a~ 2020 Jun t(N(8.8, 0.0015)) 6801. [ 6463.832, 7144.438]80
## 28 ets1 2020 Jun N(7004, 80795) 7004. [ 6639.992, 7368.540]80
## 29 ets2 2020 Jun t(N(8.9, 0.0016)) 7084. [ 6724.487, 7449.862]80
## 30 stl+e~ 2020 Jun t(N(8.8, 0.0013)) 6886. [ 6567.864, 7209.178]80
## 31 arima 2020 Jul t(N(8.9, 0.0022)) 6999. [ 6581.559, 7425.958]80
## 32 stl+a~ 2020 Jul t(N(8.8, 0.0018)) 6852. [ 6486.136, 7225.992]80
## 33 ets1 2020 Jul N(7069, 92928) 7069. [ 6678.578, 7459.916]80
## 34 ets2 2020 Jul t(N(8.9, 0.0018)) 7140. [ 6757.698, 7530.770]80
## 35 stl+e~ 2020 Jul t(N(8.8, 0.0015)) 6945. [ 6604.538, 7291.452]80
## 36 arima 2020 Aug t(N(8.9, 0.0024)) 7375. [ 6913.025, 7848.094]80
## 37 stl+a~ 2020 Aug t(N(8.9, 0.002)) 7154. [ 6747.160, 7570.427]80
## 38 ets1 2020 Aug N(7410, 114176) 7410. [ 6977.349, 7843.421]80
## 39 ets2 2020 Aug t(N(8.9, 0.002)) 7499. [ 7074.947, 7932.716]80
## 40 stl+e~ 2020 Aug t(N(8.9, 0.0017)) 7280. [ 6903.776, 7663.636]80
## 41 arima 2020 Sep t(N(8.8, 0.0027)) 6978. [ 6518.824, 7448.812]80
## 42 stl+a~ 2020 Sep t(N(8.8, 0.0023)) 6830. [ 6414.565, 7255.504]80
## 43 ets1 2020 Sep N(7017, 113547) 7017. [ 6585.239, 7448.920]80
## 44 ets2 2020 Sep t(N(8.9, 0.0022)) 7112. [ 6694.156, 7540.051]80
## 45 stl+e~ 2020 Sep t(N(8.8, 0.0018)) 6952. [ 6575.461, 7337.320]80
## 46 arima 2020 Oct t(N(8.8, 0.0029)) 6944. [ 6467.763, 7433.417]80
## 47 stl+a~ 2020 Oct t(N(8.8, 0.0026)) 6737. [ 6303.601, 7182.247]80
## 48 ets1 2020 Oct N(6994, 124244) 6994. [ 6542.033, 7445.483]80
## 49 ets2 2020 Oct t(N(8.9, 0.0023)) 7079. [ 6646.065, 7523.055]80
## 50 stl+e~ 2020 Oct t(N(8.8, 0.002)) 6875. [ 6485.430, 7273.024]80
## 51 arima 2020 Nov t(N(9.1, 0.0032)) 9083. [ 8435.657, 9749.664]80
## 52 stl+a~ 2020 Nov t(N(9.1, 0.0029)) 8822. [ 8221.987, 9437.674]80
## 53 ets1 2020 Nov N(9127, 231707) 9127. [ 8509.911, 9743.686]80
## 54 ets2 2020 Nov t(N(9.1, 0.0027)) 9228. [ 8624.445, 9847.521]80
## 55 stl+e~ 2020 Nov t(N(9.1, 0.0022)) 9022. [ 8489.917, 9566.627]80
## 56 arima 2020 Dec t(N(9.3, 0.0034)) 10686. [ 9895.583, 11500.170]80
## 57 stl+a~ 2020 Dec t(N(9.3, 0.0032)) 10822. [10047.037, 11620.263]80
## 58 ets1 2020 Dec N(10760, 350860) 10760. [10000.419, 11518.633]80
## 59 ets2 2020 Dec t(N(9.3, 0.003)) 10907. [10154.436, 11679.468]80
## 60 stl+e~ 2020 Dec t(N(9.3, 0.0023)) 11084. [10405.254, 11779.076]80
## 61 arima 2021 Jan t(N(8.8, 0.0044)) 6833. [ 6264.008, 7421.588]80
## 62 stl+a~ 2021 Jan t(N(8.8, 0.0036)) 6663. [ 6160.370, 7179.935]80
## 63 ets1 2021 Jan N(6944, 184902) 6944. [ 6393.404, 7495.545]80
## 64 ets2 2021 Jan t(N(8.9, 0.0033)) 7037. [ 6527.475, 7561.047]80
## 65 stl+e~ 2021 Jan t(N(8.8, 0.0025)) 6840. [ 6405.009, 7285.840]80
## 66 arima 2021 Feb t(N(8.8, 0.005)) 6402. [ 5832.264, 6992.951]80
## 67 stl+a~ 2021 Feb t(N(8.8, 0.0039)) 6537. [ 6020.724, 7069.960]80
## 68 ets1 2021 Feb N(6513, 173873) 6513. [ 5978.735, 7047.501]80
## 69 ets2 2021 Feb t(N(8.8, 0.0034)) 6659. [ 6166.809, 7165.771]80
## 70 stl+e~ 2021 Feb t(N(8.8, 0.0027)) 6723. [ 6281.616, 7176.316]80
## 71 arima 2021 Mar t(N(8.8, 0.0056)) 6984. [ 6326.328, 7666.343]80
## 72 stl+a~ 2021 Mar t(N(8.8, 0.0043)) 6867. [ 6299.640, 7452.922]80
## 73 ets1 2021 Mar N(7181, 225437) 7181. [ 6572.773, 7789.738]80
## 74 ets2 2021 Mar t(N(8.9, 0.0037)) 7345. [ 6782.003, 7925.433]80
## 75 stl+e~ 2021 Mar t(N(8.9, 0.0029)) 7075. [ 6596.416, 7567.517]80
## 76 arima 2021 Apr t(N(8.7, 0.0063)) 6284. [ 5659.713, 6933.872]80
## 77 stl+a~ 2021 Apr t(N(8.7, 0.0047)) 6120. [ 5592.675, 6665.980]80
## 78 ets1 2021 Apr N(6464, 194358) 6464. [ 5898.598, 7028.569]80
## 79 ets2 2021 Apr t(N(8.8, 0.0038)) 6609. [ 6096.105, 7138.198]80
## 80 stl+e~ 2021 Apr t(N(8.7, 0.003)) 6319. [ 5878.959, 6771.904]80
## 81 arima 2021 May t(N(8.8, 0.0068)) 6840. [ 6131.939, 7577.352]80
## 82 stl+a~ 2021 May t(N(8.8, 0.0051)) 6580. [ 5988.956, 7192.316]80
## 83 ets1 2021 May N(6991, 241509) 6991. [ 6361.423, 7621.022]80
## 84 ets2 2021 May t(N(8.9, 0.004)) 7141. [ 6568.803, 7731.675]80
## 85 stl+e~ 2021 May t(N(8.8, 0.0032)) 6806. [ 6318.909, 7307.543]80
## 86 arima 2021 Jun t(N(8.8, 0.0074)) 6814. [ 6078.803, 7581.837]80
## 87 stl+a~ 2021 Jun t(N(8.8, 0.0055)) 6644. [ 6023.647, 7288.602]80
## 88 ets1 2021 Jun N(6923, 251099) 6923. [ 6281.239, 7565.604]80
## 89 ets2 2021 Jun t(N(8.9, 0.0042)) 7093. [ 6512.379, 7692.779]80
## 90 stl+e~ 2021 Jun t(N(8.8, 0.0034)) 6886. [ 6380.513, 7407.325]80
## 91 arima 2021 Jul t(N(8.8, 0.008)) 6879. [ 6107.221, 7685.378]80
## 92 stl+a~ 2021 Jul t(N(8.8, 0.006)) 6689. [ 6039.492, 7363.309]80
## 93 ets1 2021 Jul N(6988, 270707) 6988. [ 6320.787, 7654.356]80
## 94 ets2 2021 Jul t(N(8.9, 0.0044)) 7150. [ 6551.402, 7768.183]80
## 95 stl+e~ 2021 Jul t(N(8.8, 0.0036)) 6946. [ 6422.708, 7484.473]80
## 96 arima 2021 Aug t(N(8.9, 0.0086)) 7249. [ 6407.319, 8130.356]80
## 97 stl+a~ 2021 Aug t(N(8.9, 0.0064)) 6998. [ 6293.833, 7732.082]80
## 98 ets1 2021 Aug N(7325, 314329) 7325. [ 6606.180, 8043.185]80
## 99 ets2 2021 Aug t(N(8.9, 0.0047)) 7509. [ 6864.470, 8176.242]80
## 100 stl+e~ 2021 Aug t(N(8.9, 0.0037)) 7281. [ 6719.828, 7859.626]80
## 101 arima 2021 Sep t(N(8.8, 0.0092)) 6860. [ 6036.459, 7723.873]80
## 102 stl+a~ 2021 Sep t(N(8.8, 0.0069)) 6671. [ 5974.871, 7396.542]80
## 103 ets1 2021 Sep N(6936, 3e+05) 6936. [ 6236.967, 7634.723]80
## 104 ets2 2021 Sep t(N(8.9, 0.0048)) 7122. [ 6501.871, 7763.316]80
## 105 stl+e~ 2021 Sep t(N(8.8, 0.0039)) 6954. [ 6405.507, 7519.066]80
## 106 arima 2021 Oct t(N(8.8, 0.0098)) 6828. [ 5983.126, 7716.116]80
## 107 stl+a~ 2021 Oct t(N(8.8, 0.0074)) 6584. [ 5873.207, 7326.865]80
## 108 ets1 2021 Oct N(6913, 311286) 6913. [ 6197.696, 7627.729]80
## 109 ets2 2021 Oct t(N(8.9, 0.005)) 7089. [ 6460.427, 7739.503]80
## 110 stl+e~ 2021 Oct t(N(8.8, 0.0041)) 6876. [ 6322.496, 7447.901]80
## 111 arima 2021 Nov t(N(9.1, 0.01)) 8933. [ 7795.347, 10130.503]80
## 112 stl+a~ 2021 Nov t(N(9.1, 0.008)) 8625. [ 7661.554, 9632.081]80
## 113 ets1 2021 Nov N(9021, 557880) 9021. [ 8063.721, 9978.137]80
## 114 ets2 2021 Nov t(N(9.1, 0.0055)) 9241. [ 8382.787, 10131.774]80
## 115 stl+e~ 2021 Nov t(N(9.1, 0.0043)) 9024. [ 8282.237, 9790.343]80
## 116 arima 2021 Dec t(N(9.3, 0.011)) 10511. [ 9135.560, 11961.462]80
## 117 stl+a~ 2021 Dec t(N(9.3, 0.0085)) 10576. [ 9355.957, 11854.161]80
## 118 ets1 2021 Dec N(10635, 814954) 10635. [ 9477.674, 11791.511]80
## 119 ets2 2021 Dec t(N(9.3, 0.0059)) 10922. [ 9871.791, 12014.317]80
## 120 stl+e~ 2021 Dec t(N(9.3, 0.0044)) 11086. [10157.040, 12047.429]80
## # ... with 1 more variable: `95%` <hilo>Cross validation: Which model is best?
tictoc::tic()
acc <- RSEAS %>%
stretch_tsibble(.init=225) %>%
model(
arima=ARIMA(log(value)~0+pdq(4,1,0)+PDQ(0,1,1)),
`stl+arima`=decomposition_model(STL(log(value)),
ARIMA(season_adjust~0+pdq(4,2,1)+PDQ(0,0,0))),
ets1=ETS(value~error('M')+trend('A')+season('M')),
ets2=ETS(log(value)~error('M')+trend('Ad')+season('M')),
`stl+ets`=decomposition_model(STL(log(value)),
ETS(season_adjust~error('A')+trend('Ad')+season('N')))
) %>%
forecast(h=24) %>%
group_by(.model,.id) %>%
mutate(h=row_number()) %>%
ungroup() %>%
accuracy(RSEAS, by = c('h','.model'), measures=list(rmse=RMSE))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 24 observations are missing between 2020 Jan and 2021 Dectictoc::toc()## 272.72 sec elapsedacc %>% print(n=Inf)## # A tibble: 120 x 4
## h .model .type rmse
## <int> <chr> <chr> <dbl>
## 1 1 arima Test 238.
## 2 1 ets1 Test 258.
## 3 1 ets2 Test 243.
## 4 1 stl+arima Test 348.
## 5 1 stl+ets Test 331.
## 6 2 arima Test 271.
## 7 2 ets1 Test 310.
## 8 2 ets2 Test 279.
## 9 2 stl+arima Test 352.
## 10 2 stl+ets Test 348.
## 11 3 arima Test 295.
## 12 3 ets1 Test 335.
## 13 3 ets2 Test 294.
## 14 3 stl+arima Test 380.
## 15 3 stl+ets Test 368.
## 16 4 arima Test 308.
## 17 4 ets1 Test 354.
## 18 4 ets2 Test 307.
## 19 4 stl+arima Test 359.
## 20 4 stl+ets Test 355.
## 21 5 arima Test 315.
## 22 5 ets1 Test 367.
## 23 5 ets2 Test 302.
## 24 5 stl+arima Test 372.
## 25 5 stl+ets Test 357.
## 26 6 arima Test 324.
## 27 6 ets1 Test 370.
## 28 6 ets2 Test 308.
## 29 6 stl+arima Test 375.
## 30 6 stl+ets Test 363.
## 31 7 arima Test 323.
## 32 7 ets1 Test 397.
## 33 7 ets2 Test 314.
## 34 7 stl+arima Test 387.
## 35 7 stl+ets Test 369.
## 36 8 arima Test 331.
## 37 8 ets1 Test 418.
## 38 8 ets2 Test 328.
## 39 8 stl+arima Test 384.
## 40 8 stl+ets Test 375.
## 41 9 arima Test 351.
## 42 9 ets1 Test 441.
## 43 9 ets2 Test 342.
## 44 9 stl+arima Test 440.
## 45 9 stl+ets Test 394.
## 46 10 arima Test 368.
## 47 10 ets1 Test 473.
## 48 10 ets2 Test 349.
## 49 10 stl+arima Test 407.
## 50 10 stl+ets Test 376.
## 51 11 arima Test 386.
## 52 11 ets1 Test 496.
## 53 11 ets2 Test 354.
## 54 11 stl+arima Test 425.
## 55 11 stl+ets Test 383.
## 56 12 arima Test 397.
## 57 12 ets1 Test 468.
## 58 12 ets2 Test 357.
## 59 12 stl+arima Test 406.
## 60 12 stl+ets Test 379.
## 61 13 arima Test 482.
## 62 13 ets1 Test 554.
## 63 13 ets2 Test 423.
## 64 13 stl+arima Test 578.
## 65 13 stl+ets Test 509.
## 66 14 arima Test 493.
## 67 14 ets1 Test 602.
## 68 14 ets2 Test 429.
## 69 14 stl+arima Test 569.
## 70 14 stl+ets Test 513.
## 71 15 arima Test 515.
## 72 15 ets1 Test 606.
## 73 15 ets2 Test 439.
## 74 15 stl+arima Test 583.
## 75 15 stl+ets Test 523.
## 76 16 arima Test 529.
## 77 16 ets1 Test 653.
## 78 16 ets2 Test 446.
## 79 16 stl+arima Test 561.
## 80 16 stl+ets Test 503.
## 81 17 arima Test 537.
## 82 17 ets1 Test 672.
## 83 17 ets2 Test 442.
## 84 17 stl+arima Test 577.
## 85 17 stl+ets Test 506.
## 86 18 arima Test 539.
## 87 18 ets1 Test 665.
## 88 18 ets2 Test 452.
## 89 18 stl+arima Test 565.
## 90 18 stl+ets Test 506.
## 91 19 arima Test 520.
## 92 19 ets1 Test 676.
## 93 19 ets2 Test 444.
## 94 19 stl+arima Test 571.
## 95 19 stl+ets Test 503.
## 96 20 arima Test 516.
## 97 20 ets1 Test 679.
## 98 20 ets2 Test 449.
## 99 20 stl+arima Test 569.
## 100 20 stl+ets Test 509.
## 101 21 arima Test 550.
## 102 21 ets1 Test 716.
## 103 21 ets2 Test 462.
## 104 21 stl+arima Test 656.
## 105 21 stl+ets Test 531.
## 106 22 arima Test 573.
## 107 22 ets1 Test 787.
## 108 22 ets2 Test 479.
## 109 22 stl+arima Test 615.
## 110 22 stl+ets Test 516.
## 111 23 arima Test 606.
## 112 23 ets1 Test 868.
## 113 23 ets2 Test 500.
## 114 23 stl+arima Test 640.
## 115 23 stl+ets Test 528.
## 116 24 arima Test 625.
## 117 24 ets1 Test 797.
## 118 24 ets2 Test 513.
## 119 24 stl+arima Test 638.
## 120 24 stl+ets Test 537.acc %>%
ggplot(aes(x=h,y=rmse,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:24, minor_breaks = NULL)+
labs(y='RMSE',title='RSEAS: Point Prediction Accuracy')
Let’s try distributional accuracy but be prepared to wait a while…
tictoc::tic()
acc <- RSEAS %>%
stretch_tsibble(.init=225) %>%
model(
arima=ARIMA(log(value)~0+pdq(4,1,0)+PDQ(0,1,1)),
`stl+arima`=decomposition_model(STL(log(value)),
ARIMA(season_adjust~0+pdq(4,2,1)+PDQ(0,0,0))),
ets1=ETS(value~error('M')+trend('A')+season('M')),
ets2=ETS(log(value)~error('M')+trend('Ad')+season('M')),
`stl+ets`=decomposition_model(STL(log(value)),
ETS(season_adjust~error('A')+trend('Ad')+season('N')))
) %>%
forecast(h=24) %>%
group_by(.model,.id) %>%
mutate(h=row_number()) %>%
ungroup() %>%
accuracy(RSEAS, by = c('h','.model'), measures=list(rmse=RMSE,crps=CRPS))## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 24 observations are missing between 2020 Jan and 2021 Dectictoc::toc()## 2695.11 sec elapsedacc %>% print(n=Inf)## # A tibble: 120 x 5
## h .model .type rmse crps
## <int> <chr> <chr> <dbl> <dbl>
## 1 1 arima Test 238. 129.
## 2 1 ets1 Test 258. 137.
## 3 1 ets2 Test 243. 130.
## 4 1 stl+arima Test 348. 183.
## 5 1 stl+ets Test 331. 170.
## 6 2 arima Test 271. 147.
## 7 2 ets1 Test 310. 164.
## 8 2 ets2 Test 279. 151.
## 9 2 stl+arima Test 352. 190.
## 10 2 stl+ets Test 348. 187.
## 11 3 arima Test 295. 164.
## 12 3 ets1 Test 335. 184.
## 13 3 ets2 Test 294. 162.
## 14 3 stl+arima Test 380. 209.
## 15 3 stl+ets Test 368. 200.
## 16 4 arima Test 308. 170.
## 17 4 ets1 Test 354. 196.
## 18 4 ets2 Test 307. 168.
## 19 4 stl+arima Test 359. 197.
## 20 4 stl+ets Test 355. 191.
## 21 5 arima Test 315. 176.
## 22 5 ets1 Test 367. 203.
## 23 5 ets2 Test 302. 168.
## 24 5 stl+arima Test 372. 208.
## 25 5 stl+ets Test 357. 195.
## 26 6 arima Test 324. 181.
## 27 6 ets1 Test 370. 207.
## 28 6 ets2 Test 308. 172.
## 29 6 stl+arima Test 375. 208.
## 30 6 stl+ets Test 363. 195.
## 31 7 arima Test 323. 185.
## 32 7 ets1 Test 397. 221.
## 33 7 ets2 Test 314. 178.
## 34 7 stl+arima Test 387. 215.
## 35 7 stl+ets Test 369. 200.
## 36 8 arima Test 331. 190.
## 37 8 ets1 Test 418. 232.
## 38 8 ets2 Test 328. 186.
## 39 8 stl+arima Test 384. 216.
## 40 8 stl+ets Test 375. 206.
## 41 9 arima Test 351. 198.
## 42 9 ets1 Test 441. 240.
## 43 9 ets2 Test 342. 192.
## 44 9 stl+arima Test 440. 237.
## 45 9 stl+ets Test 394. 215.
## 46 10 arima Test 368. 208.
## 47 10 ets1 Test 473. 251.
## 48 10 ets2 Test 349. 196.
## 49 10 stl+arima Test 407. 226.
## 50 10 stl+ets Test 376. 207.
## 51 11 arima Test 386. 216.
## 52 11 ets1 Test 496. 258.
## 53 11 ets2 Test 354. 198.
## 54 11 stl+arima Test 425. 233.
## 55 11 stl+ets Test 383. 208.
## 56 12 arima Test 397. 224.
## 57 12 ets1 Test 468. 257.
## 58 12 ets2 Test 357. 199.
## 59 12 stl+arima Test 406. 225.
## 60 12 stl+ets Test 379. 203.
## 61 13 arima Test 482. 264.
## 62 13 ets1 Test 554. 300.
## 63 13 ets2 Test 423. 233.
## 64 13 stl+arima Test 578. 297.
## 65 13 stl+ets Test 509. 260.
## 66 14 arima Test 493. 273.
## 67 14 ets1 Test 602. 319.
## 68 14 ets2 Test 429. 239.
## 69 14 stl+arima Test 569. 303.
## 70 14 stl+ets Test 513. 266.
## 71 15 arima Test 515. 289.
## 72 15 ets1 Test 606. 338.
## 73 15 ets2 Test 439. 248.
## 74 15 stl+arima Test 583. 318.
## 75 15 stl+ets Test 523. 275.
## 76 16 arima Test 529. 298.
## 77 16 ets1 Test 653. 362.
## 78 16 ets2 Test 446. 252.
## 79 16 stl+arima Test 561. 312.
## 80 16 stl+ets Test 503. 270.
## 81 17 arima Test 537. 302.
## 82 17 ets1 Test 672. 369.
## 83 17 ets2 Test 442. 251.
## 84 17 stl+arima Test 577. 322.
## 85 17 stl+ets Test 506. 272.
## 86 18 arima Test 539. 303.
## 87 18 ets1 Test 665. 372.
## 88 18 ets2 Test 452. 256.
## 89 18 stl+arima Test 565. 315.
## 90 18 stl+ets Test 506. 269.
## 91 19 arima Test 520. 300.
## 92 19 ets1 Test 676. 376.
## 93 19 ets2 Test 444. 254.
## 94 19 stl+arima Test 571. 319.
## 95 19 stl+ets Test 503. 267.
## 96 20 arima Test 516. 303.
## 97 20 ets1 Test 679. 381.
## 98 20 ets2 Test 449. 258.
## 99 20 stl+arima Test 569. 322.
## 100 20 stl+ets Test 509. 272.
## 101 21 arima Test 550. 317.
## 102 21 ets1 Test 716. 402.
## 103 21 ets2 Test 462. 264.
## 104 21 stl+arima Test 656. 355.
## 105 21 stl+ets Test 531. 284.
## 106 22 arima Test 573. 326.
## 107 22 ets1 Test 787. 416.
## 108 22 ets2 Test 479. 268.
## 109 22 stl+arima Test 615. 344.
## 110 22 stl+ets Test 516. 277.
## 111 23 arima Test 606. 337.
## 112 23 ets1 Test 868. 438.
## 113 23 ets2 Test 500. 273.
## 114 23 stl+arima Test 640. 351.
## 115 23 stl+ets Test 528. 279.
## 116 24 arima Test 625. 348.
## 117 24 ets1 Test 797. 440.
## 118 24 ets2 Test 513. 276.
## 119 24 stl+arima Test 638. 349.
## 120 24 stl+ets Test 537. 278.acc %>%
ggplot(aes(x=h,y=rmse,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:24, minor_breaks = NULL)+
labs(y='RMSE',title='RSEAS: Point Prediction Accuracy')
acc %>%
ggplot(aes(x=h,y=crps,color=.model))+
geom_point()+
guides(color=guide_legend(title='Forecast'))+
scale_x_continuous(breaks=1:24, minor_breaks = NULL)+
labs(y='CRPS',title='RSEAS: Distribution Accuracy')