ADAM ETS
First and foremost, ADAM implements ETS model, although in a more
flexible way than (Hyndman et al. 2008):
it supports different distributions for the error term, which are
regulated via distribution parameter. By default, the
additive error model relies on Normal distribution, while the
multiplicative error one assumes Inverse Gaussian. If you want to
reproduce the classical ETS, you would need to specify
distribution="dnorm". Here is an example of ADAM ETS(MMM)
with Normal distribution on AirPassengers data:
testModel <- adam(AirPassengers, "MMM", lags=c(1,12), distribution="dnorm",
h=12, holdout=TRUE)
summary(testModel)
#>
#> Model estimated using adam() function: ETS(MMM)
#> Response variable: AirPassengers
#> Distribution used in the estimation: Normal
#> Loss function type: likelihood; Loss function value: 467.4621
#> Coefficients:
#> Estimate Std. Error Lower 2.5% Upper 97.5%
#> alpha 0.7500 0.0899 0.5721 0.9278 *
#> beta 0.0055 0.0041 0.0000 0.0136
#> gamma 0.0000 0.0145 0.0000 0.0286
#>
#> Error standard deviation: 0.0352
#> Sample size: 132
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 128
#> Information criteria:
#> AIC AICc BIC BICc
#> 942.9242 943.2392 954.4554 955.2244
plot(forecast(testModel,h=12,interval="prediction"))
You might notice that the summary contains more than what is reported
by other smooth functions. This one also produces standard
errors for the estimated parameters based on Fisher Information
calculation. Note that this is computationally expensive, so if you have
a model with more than 30 variables, the calculation of standard errors
might take plenty of time. As for the default print()
method, it will produce a shorter summary from the model, without the
standard errors (similar to what es() does):
testModel
#> Time elapsed: 0.06 seconds
#> Model estimated using adam() function: ETS(MMM)
#> With backcasting initialisation
#> Distribution assumed in the model: Normal
#> Loss function type: likelihood; Loss function value: 467.4621
#> Persistence vector g:
#> alpha beta gamma
#> 0.7500 0.0055 0.0000
#>
#> Sample size: 132
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 128
#> Information criteria:
#> AIC AICc BIC BICc
#> 942.9242 943.2392 954.4554 955.2244
#>
#> Forecast errors:
#> ME: -4.77; MAE: 15.396; RMSE: 21.749
#> sCE: -21.807%; Asymmetry: -16.3%; sMAE: 5.865%; sMSE: 0.687%
#> MASE: 0.639; RMSSE: 0.694; rMAE: 0.203; rRMSE: 0.211Also, note that the prediction interval in case of multiplicative error models are approximate. It is advisable to use simulations instead (which is slower, but more accurate):
If you want to do the residuals diagnostics, then it is recommended
to use plot function, something like this (you can select,
which of the plots to produce):
By default ADAM will estimate models via maximising likelihood
function. But there is also a parameter loss, which allows
selecting from a list of already implemented loss functions (again, see
documentation for adam() for the full list) or using a
function written by a user. Here is how to do the latter on the example
of BJsales:
lossFunction <- function(actual, fitted, B){
return(sum(abs(actual-fitted)^3))
}
testModel <- adam(BJsales, "AAN", silent=FALSE, loss=lossFunction,
h=12, holdout=TRUE)
testModel
#> Time elapsed: 0.02 seconds
#> Model estimated using adam() function: ETS(AAN)
#> With backcasting initialisation
#> Distribution assumed in the model: Normal
#> Loss function type: custom; Loss function value: 599.299
#> Persistence vector g:
#> alpha beta
#> 1.000 0.227
#>
#> Sample size: 138
#> Number of estimated parameters: 2
#> Number of degrees of freedom: 136
#> Information criteria are unavailable for the chosen loss & distribution.
#>
#> Forecast errors:
#> ME: 3.015; MAE: 3.128; RMSE: 3.866
#> sCE: 15.915%; Asymmetry: 91.7%; sMAE: 1.376%; sMSE: 0.029%
#> MASE: 2.626; RMSSE: 2.52; rMAE: 1.009; rRMSE: 1.009
Note that you need to have parameters actual, fitted and B in the function, which correspond to the vector of actual values, vector of fitted values on each iteration and a vector of the optimised parameters.
loss and distribution parameters are
independent, so in the example above, we have assumed that the error
term follows Normal distribution, but we have estimated its parameters
using a non-conventional loss because we can. Some of distributions
assume that there is an additional parameter, which can either be
estimated or provided by user. These include Asymmetric Laplace
(distribution="dalaplace") with alpha,
Generalised Normal and Log-Generalised normal
(distribution=c("gnorm","dlgnorm")) with shape
and Student’s T (distribution="dt") with
nu:
The model selection in ADAM ETS relies on information criteria and
works correctly only for the loss="likelihood". There are
several options, how to select the model, see them in the description of
the function: ?adam(). The default one uses
branch-and-bound algorithm, similar to the one used in
es(), but only considers additive trend models (the
multiplicative trend ones are less stable and need more attention from a
forecaster):
testModel <- adam(AirPassengers, "ZXZ", lags=c(1,12), silent=FALSE,
h=12, holdout=TRUE)
#> Forming the pool of models based on... ANN , ANA , MNM , MAM , Estimation progress: 71 %86 %100 %... Done!
testModel
#> Time elapsed: 0.2 seconds
#> Model estimated using adam() function: ETS(MAM)
#> With backcasting initialisation
#> Distribution assumed in the model: Gamma
#> Loss function type: likelihood; Loss function value: 468.1117
#> Persistence vector g:
#> alpha beta gamma
#> 0.7976 0.0030 0.0090
#>
#> Sample size: 132
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 128
#> Information criteria:
#> AIC AICc BIC BICc
#> 944.2235 944.5384 955.7547 956.5236
#>
#> Forecast errors:
#> ME: 10.11; MAE: 20.814; RMSE: 26.223
#> sCE: 46.217%; Asymmetry: 66%; sMAE: 7.929%; sMSE: 0.998%
#> MASE: 0.864; RMSSE: 0.837; rMAE: 0.274; rRMSE: 0.255
Note that the function produces point forecasts if
h>0, but it won’t generate prediction interval. This is
why you need to use forecast() method (as shown in the
first example in this vignette).
Similarly to es(), function supports combination of
models, but it saves all the tested models in the output for a potential
reuse. Here how it works:
testModel <- adam(AirPassengers, "CXC", lags=c(1,12),
h=12, holdout=TRUE)
testForecast <- forecast(testModel,h=18,interval="semiparametric", level=c(0.9,0.95))
testForecast
#> Point forecast Lower bound (5%) Lower bound (2.5%) Upper bound (95%)
#> Jan 1960 413.0680 390.0118 385.7252 436.6783
#> Feb 1960 405.8957 377.3811 372.1223 435.2783
#> Mar 1960 468.3322 428.6131 421.3483 509.5205
#> Apr 1960 451.1823 409.8723 402.3457 494.1471
#> May 1960 453.9584 411.3825 403.6353 498.2832
#> Jun 1960 516.5757 466.0950 456.9305 569.2208
#> Jul 1960 570.8187 513.0840 502.6233 631.1197
#> Aug 1960 568.8944 510.2984 499.6933 630.1450
#> Sep 1960 498.8337 447.6481 438.3821 552.3291
#> Oct 1960 435.8123 390.7898 382.6427 482.8809
#> Nov 1960 380.1682 340.7948 333.6711 421.3356
#> Dec 1960 429.1452 383.8254 375.6356 476.5726
#> Jan 1961 437.0926 388.6635 379.9386 487.8906
#> Feb 1961 429.3893 379.3480 370.3634 482.0120
#> Mar 1961 495.3095 433.2238 422.1353 560.8545
#> Apr 1961 477.0476 414.3553 403.1995 543.4127
#> May 1961 479.8591 416.1189 404.7864 547.3763
#> Jun 1961 545.9096 473.0691 460.1234 623.0870
#> Upper bound (97.5%)
#> Jan 1960 441.3339
#> Feb 1960 441.1152
#> Mar 1960 517.7637
#> Apr 1960 502.7756
#> May 1960 507.1950
#> Jun 1960 579.8267
#> Jul 1960 643.2891
#> Aug 1960 642.5178
#> Sep 1960 563.1332
#> Oct 1960 492.3904
#> Nov 1960 429.6540
#> Dec 1960 486.1658
#> Jan 1961 498.1928
#> Feb 1961 492.7156
#> Mar 1961 574.2462
#> Apr 1961 557.0139
#> May 1961 561.2236
#> Jun 1961 638.9203
plot(testForecast)
Yes, now we support vectors for the levels in case you want to produce several. In fact, we also support side for prediction interval, so you can extract specific quantiles without a hustle:
forecast(testModel,h=18,interval="semiparametric", level=c(0.9,0.95,0.99), side="upper")
#> Point forecast Upper bound (90%) Upper bound (95%) Upper bound (99%)
#> Jan 1960 413.0680 431.3514 436.6783 446.7882
#> Feb 1960 405.8957 428.6126 435.2783 447.9667
#> Mar 1960 468.3322 500.1249 509.5205 527.4586
#> Apr 1960 451.1823 484.3213 494.1471 512.9326
#> May 1960 453.9584 488.1378 498.2832 517.6884
#> Jun 1960 516.5757 557.1531 569.2208 592.3215
#> Jul 1960 570.8187 617.2791 631.1197 657.6323
#> Aug 1960 568.8944 616.0766 630.1450 657.1042
#> Sep 1960 498.8337 540.0438 552.3291 575.8695
#> Oct 1960 435.8123 472.0687 482.8809 503.6017
#> Nov 1960 380.1682 411.8780 421.3356 439.4613
#> Dec 1960 429.1452 465.6685 476.5726 497.4792
#> Jan 1961 437.0926 476.1887 487.8906 510.3505
#> Feb 1961 429.3893 469.8635 482.0120 505.3562
#> Mar 1961 495.3095 545.6723 560.8545 590.0797
#> Apr 1961 477.0476 528.0054 543.4127 573.1076
#> May 1961 479.8591 531.6932 547.3763 577.6114
#> Jun 1961 545.9096 605.1560 623.0870 657.6600A brand new thing in the function is the possibility to use several
frequencies (double / triple / quadruple / … seasonal models). In order
to show how it works, we will generate an artificial time series,
inspired by half-hourly electricity demand using sim.gum()
function:
set.seed(41)
ordersGUM <- c(1,1,1)
lagsGUM <- c(1,48,336)
initialGUM1 <- -25381.7
initialGUM2 <- c(23955.09, 24248.75, 24848.54, 25012.63, 24634.14, 24548.22, 24544.63, 24572.77,
24498.33, 24250.94, 24545.44, 25005.92, 26164.65, 27038.55, 28262.16, 28619.83,
28892.19, 28575.07, 28837.87, 28695.12, 28623.02, 28679.42, 28682.16, 28683.40,
28647.97, 28374.42, 28261.56, 28199.69, 28341.69, 28314.12, 28252.46, 28491.20,
28647.98, 28761.28, 28560.11, 28059.95, 27719.22, 27530.23, 27315.47, 27028.83,
26933.75, 26961.91, 27372.44, 27362.18, 27271.31, 26365.97, 25570.88, 25058.01)
initialGUM3 <- c(23920.16, 23026.43, 22812.23, 23169.52, 23332.56, 23129.27, 22941.20, 22692.40,
22607.53, 22427.79, 22227.64, 22580.72, 23871.99, 25758.34, 28092.21, 30220.46,
31786.51, 32699.80, 33225.72, 33788.82, 33892.25, 34112.97, 34231.06, 34449.53,
34423.61, 34333.93, 34085.28, 33948.46, 33791.81, 33736.17, 33536.61, 33633.48,
33798.09, 33918.13, 33871.41, 33403.75, 32706.46, 31929.96, 31400.48, 30798.24,
29958.04, 30020.36, 29822.62, 30414.88, 30100.74, 29833.49, 28302.29, 26906.72,
26378.64, 25382.11, 25108.30, 25407.07, 25469.06, 25291.89, 25054.11, 24802.21,
24681.89, 24366.97, 24134.74, 24304.08, 25253.99, 26950.23, 29080.48, 31076.33,
32453.20, 33232.81, 33661.61, 33991.21, 34017.02, 34164.47, 34398.01, 34655.21,
34746.83, 34596.60, 34396.54, 34236.31, 34153.32, 34102.62, 33970.92, 34016.13,
34237.27, 34430.08, 34379.39, 33944.06, 33154.67, 32418.62, 31781.90, 31208.69,
30662.59, 30230.67, 30062.80, 30421.11, 30710.54, 30239.27, 28949.56, 27506.96,
26891.75, 25946.24, 25599.88, 25921.47, 26023.51, 25826.29, 25548.72, 25405.78,
25210.45, 25046.38, 24759.76, 24957.54, 25815.10, 27568.98, 29765.24, 31728.25,
32987.51, 33633.74, 34021.09, 34407.19, 34464.65, 34540.67, 34644.56, 34756.59,
34743.81, 34630.05, 34506.39, 34319.61, 34110.96, 33961.19, 33876.04, 33969.95,
34220.96, 34444.66, 34474.57, 34018.83, 33307.40, 32718.90, 32115.27, 31663.53,
30903.82, 31013.83, 31025.04, 31106.81, 30681.74, 30245.70, 29055.49, 27582.68,
26974.67, 25993.83, 25701.93, 25940.87, 26098.63, 25771.85, 25468.41, 25315.74,
25131.87, 24913.15, 24641.53, 24807.15, 25760.85, 27386.39, 29570.03, 31634.00,
32911.26, 33603.94, 34020.90, 34297.65, 34308.37, 34504.71, 34586.78, 34725.81,
34765.47, 34619.92, 34478.54, 34285.00, 34071.90, 33986.48, 33756.85, 33799.37,
33987.95, 34047.32, 33924.48, 33580.82, 32905.87, 32293.86, 31670.02, 31092.57,
30639.73, 30245.42, 30281.61, 30484.33, 30349.51, 29889.23, 28570.31, 27185.55,
26521.85, 25543.84, 25187.82, 25371.59, 25410.07, 25077.67, 24741.93, 24554.62,
24427.19, 24127.21, 23887.55, 24028.40, 24981.34, 26652.32, 28808.00, 30847.09,
32304.13, 33059.02, 33562.51, 33878.96, 33976.68, 34172.61, 34274.50, 34328.71,
34370.12, 34095.69, 33797.46, 33522.96, 33169.94, 32883.32, 32586.24, 32380.84,
32425.30, 32532.69, 32444.24, 32132.49, 31582.39, 30926.58, 30347.73, 29518.04,
29070.95, 28586.20, 28416.94, 28598.76, 28529.75, 28424.68, 27588.76, 26604.13,
26101.63, 25003.82, 24576.66, 24634.66, 24586.21, 24224.92, 23858.42, 23577.32,
23272.28, 22772.00, 22215.13, 21987.29, 21948.95, 22310.79, 22853.79, 24226.06,
25772.55, 27266.27, 28045.65, 28606.14, 28793.51, 28755.83, 28613.74, 28376.47,
27900.76, 27682.75, 27089.10, 26481.80, 26062.94, 25717.46, 25500.27, 25171.05,
25223.12, 25634.63, 26306.31, 26822.46, 26787.57, 26571.18, 26405.21, 26148.41,
25704.47, 25473.10, 25265.97, 26006.94, 26408.68, 26592.04, 26224.64, 25407.27,
25090.35, 23930.21, 23534.13, 23585.75, 23556.93, 23230.25, 22880.24, 22525.52,
22236.71, 21715.08, 21051.17, 20689.40, 20099.18, 19939.71, 19722.69, 20421.58,
21542.03, 22962.69, 23848.69, 24958.84, 25938.72, 26316.56, 26742.61, 26990.79,
27116.94, 27168.78, 26464.41, 25703.23, 25103.56, 24891.27, 24715.27, 24436.51,
24327.31, 24473.02, 24893.89, 25304.13, 25591.77, 25653.00, 25897.55, 25859.32,
25918.32, 25984.63, 26232.01, 26810.86, 27209.70, 26863.50, 25734.54, 24456.96)
y <- sim.gum(orders=ordersGUM, lags=lagsGUM, nsim=1, frequency=336, obs=3360,
measurement=rep(1,3), transition=diag(3), persistence=c(0.045,0.162,0.375),
initial=cbind(initialGUM1,initialGUM2,initialGUM3))$dataWe can then apply ADAM to this data:
testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
silent=FALSE, h=336, holdout=TRUE)
testModel
#> Time elapsed: 1.21 seconds
#> Model estimated using adam() function: ETS(MMdM)[48, 336]
#> With backcasting initialisation
#> Distribution assumed in the model: Gamma
#> Loss function type: likelihood; Loss function value: 19552.64
#> Persistence vector g:
#> alpha beta gamma1 gamma2
#> 0.0276 0.0000 0.1859 0.2400
#> Damping parameter: 0.4816
#> Sample size: 3024
#> Number of estimated parameters: 6
#> Number of degrees of freedom: 3018
#> Information criteria:
#> AIC AICc BIC BICc
#> 39117.28 39117.31 39153.37 39153.48
#>
#> Forecast errors:
#> ME: 80.315; MAE: 145.549; RMSE: 186.929
#> sCE: 89.037%; Asymmetry: 57.1%; sMAE: 0.48%; sMSE: 0.004%
#> MASE: 0.196; RMSSE: 0.182; rMAE: 0.023; rRMSE: 0.024
Note that the more lags you have, the more initial seasonal
components the function will need to estimate, which is a difficult
task. This is why we used initial="backcasting" in the
example above - this speeds up the estimation by reducing the number of
parameters to estimate. Still, the optimiser might not get close to the
optimal value, so we can help it. First, we can give more time for the
calculation, increasing the number of iterations via
maxeval (the default value is 40 iterations for each
estimated parameter, e.g. \(40 \times 5 =
200\) in our case):
testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
silent=FALSE, h=336, holdout=TRUE, maxeval=10000)
testModel
#> Time elapsed: 1.32 seconds
#> Model estimated using adam() function: ETS(MMdM)[48, 336]
#> With backcasting initialisation
#> Distribution assumed in the model: Gamma
#> Loss function type: likelihood; Loss function value: 19552.64
#> Persistence vector g:
#> alpha beta gamma1 gamma2
#> 0.0275 0.0000 0.1860 0.2397
#> Damping parameter: 0.4814
#> Sample size: 3024
#> Number of estimated parameters: 6
#> Number of degrees of freedom: 3018
#> Information criteria:
#> AIC AICc BIC BICc
#> 39117.28 39117.31 39153.37 39153.48
#>
#> Forecast errors:
#> ME: 80.314; MAE: 145.555; RMSE: 186.937
#> sCE: 89.036%; Asymmetry: 57.1%; sMAE: 0.48%; sMSE: 0.004%
#> MASE: 0.196; RMSSE: 0.182; rMAE: 0.023; rRMSE: 0.024
This will take more time, but will typically lead to more refined
parameters. You can control other parameters of the optimiser as well,
such as algorithm, xtol_rel,
print_level and others, which are explained in the
documentation for nloptr function from nloptr package (run
nloptr.print.options() for details). Second, we can give a
different set of initial parameters for the optimiser, have a look at
what the function saves:
and use this as a starting point for the reestimation (e.g. with a different algorithm):
testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
silent=FALSE, h=336, holdout=TRUE, B=testModel$B)
testModel
#> Time elapsed: 0.43 seconds
#> Model estimated using adam() function: ETS(MMdM)[48, 336]
#> With backcasting initialisation
#> Distribution assumed in the model: Gamma
#> Loss function type: likelihood; Loss function value: 19552.64
#> Persistence vector g:
#> alpha beta gamma1 gamma2
#> 0.0275 0.0000 0.1860 0.2397
#> Damping parameter: 0.7314
#> Sample size: 3024
#> Number of estimated parameters: 6
#> Number of degrees of freedom: 3018
#> Information criteria:
#> AIC AICc BIC BICc
#> 39117.28 39117.31 39153.37 39153.48
#>
#> Forecast errors:
#> ME: 80.314; MAE: 145.555; RMSE: 186.937
#> sCE: 89.036%; Asymmetry: 57.1%; sMAE: 0.48%; sMSE: 0.004%
#> MASE: 0.196; RMSSE: 0.182; rMAE: 0.023; rRMSE: 0.024
If you are ready to wait, you can change the initialisation to the
initial="optimal", which in our case will take much more
time because of the number of estimated parameters - 389 for the chosen
model. The estimation process in this case might take 20 - 30 times more
than in the example above.
In addition, you can specify some parts of the initial state vector or some parts of the persistence vector, here is an example:
testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
silent=TRUE, h=336, holdout=TRUE, persistence=list(beta=0.1))
testModel
#> Time elapsed: 1.01 seconds
#> Model estimated using adam() function: ETS(MMdM)[48, 336]
#> With backcasting initialisation
#> Distribution assumed in the model: Gamma
#> Loss function type: likelihood; Loss function value: 56395.13
#> Persistence vector g:
#> alpha beta gamma1 gamma2
#> 0.1003 0.1000 0.3081 0.3022
#> Damping parameter: 0.9503
#> Sample size: 3024
#> Number of estimated parameters: 5
#> Number of degrees of freedom: 3019
#> Information criteria:
#> AIC AICc BIC BICc
#> 112800.3 112800.3 112830.3 112830.4
#>
#> Forecast errors:
#> ME: -801836.573; MAE: 861868.38; RMSE: 11443964.786
#> sCE: -888922.037%; Asymmetry: 47.6%; sMAE: 2843.672%; sMSE: 14257063.817%
#> MASE: 1160.625; RMSSE: 11129.534; rMAE: 134.71; rRMSE: 1451.758The function also handles intermittent data (the data with zeroes) and the data with missing values. This is partially covered in the vignette on the oes() function. Here is a simple example:
testModel <- adam(rpois(120,0.5), "MNN", silent=FALSE, h=12, holdout=TRUE,
occurrence="odds-ratio")
testModel
#> Time elapsed: 0.03 seconds
#> Model estimated using adam() function: iETS(MNN)[O]
#> With backcasting initialisation
#> Occurrence model type: Odds ratio
#> Distribution assumed in the model: Mixture of Bernoulli and Gamma
#> Loss function type: likelihood; Loss function value: 23.4357
#> Persistence vector g:
#> alpha
#> 0.0306
#>
#> Sample size: 108
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 104
#> Information criteria:
#> AIC AICc BIC BICc
#> 192.3585 192.4728 203.0870 193.9903
#>
#> Forecast errors:
#> Asymmetry: -0.889%; sMSE: 28.566%; rRMSE: 0.8; sPIS: -67.368%; sCE: 98.476%
Finally, adam() is faster than es()
function, because its code is more efficient and it uses a different
optimisation algorithm with more finely tuned parameters by default.
Let’s compare:
adamModel <- adam(AirPassengers, "CCC",
h=12, holdout=TRUE)
esModel <- es(AirPassengers, "CCC",
h=12, holdout=TRUE)
"adam:"
#> [1] "adam:"
adamModel
#> Time elapsed: 0.83 seconds
#> Model estimated: ETS(CCC)
#> Loss function type: likelihood
#>
#> Number of models combined: 30
#> Sample size: 132
#> Average number of estimated parameters: 5.7052
#> Average number of degrees of freedom: 126.2948
#>
#> Forecast errors:
#> ME: 0.554; MAE: 16.037; RMSE: 22.028
#> sCE: 2.533%; sMAE: 6.109%; sMSE: 0.704%
#> MASE: 0.666; RMSSE: 0.703; rMAE: 0.211; rRMSE: 0.214
"es():"
#> [1] "es():"
esModel
#> Time elapsed: 0.78 seconds
#> Model estimated: ETS(CCC)
#> Loss function type: likelihood
#>
#> Number of models combined: 30
#> Sample size: 132
#> Average number of estimated parameters: 5.4538
#> Average number of degrees of freedom: 126.5462
#>
#> Forecast errors:
#> ME: 3.407; MAE: 16.501; RMSE: 22.884
#> sCE: 15.575%; sMAE: 6.286%; sMSE: 0.76%
#> MASE: 0.685; RMSSE: 0.73; rMAE: 0.217; rRMSE: 0.222
