From 9ac843fc51e96468175a5f98f34411ad80017eb1 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Mon, 1 Apr 2024 14:50:02 +0000 Subject: [PATCH] version 4.0.1 --- DESCRIPTION | 15 +- MD5 | 80 +++---- NAMESPACE | 1 + NEWS | 19 ++ R/RcppExports.R | 28 ++- R/adam-es.R | 82 ++++--- R/adam-sma.R | 52 ++--- R/adam.R | 400 +++++++++++++++++-------------- R/adamGeneral.R | 31 ++- R/autoadam.R | 5 + R/cma.R | 2 +- R/msdecompose.R | 2 + R/smooth-package.R | 1 - build/partial.rdb | Bin 61 -> 61 bytes build/vignette.rds | Bin 492 -> 492 bytes inst/doc/adam.html | 518 ++++++++++++++++++++--------------------- inst/doc/ces.html | 8 +- inst/doc/es.html | 160 ++++++------- inst/doc/gum.html | 34 +-- inst/doc/oes.html | 122 +++++----- inst/doc/simulate.html | 80 +++---- inst/doc/sma.html | 20 +- inst/doc/smooth.Rmd | 2 +- inst/doc/smooth.html | 14 +- inst/doc/ssarima.html | 48 ++-- man/adam.Rd | 12 +- man/cma.Rd | 9 +- man/es.Rd | 2 +- man/gum.Rd | 7 +- man/msdecompose.Rd | 8 + man/sim.gum.Rd | 7 +- man/sma.Rd | 7 +- man/smooth.Rd | 8 +- src/RcppExports.cpp | 66 ++++-- src/adamGeneral.cpp | 178 +++++++++++--- src/adamGeneral.h | 1 - src/adamRefitter.cpp | 36 +-- src/adamSimulator.cpp | 18 +- src/ssGeneral.cpp | 17 -- src/ssGeneral.h | 18 +- vignettes/smooth.Rmd | 2 +- 41 files changed, 1181 insertions(+), 939 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3922e57..5b3020d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,9 @@ Package: smooth Type: Package Title: Forecasting Using State Space Models -Version: 4.0.0 -Date: 2023-09-15 -Authors@R: person("Ivan", "Svetunkov", email = "ivan@svetunkov.ru", role = c("aut", "cre"), +Version: 4.0.1 +Date: 2024-04-01 +Authors@R: person("Ivan", "Svetunkov", email = "ivan@svetunkov.com", role = c("aut", "cre"), comment="Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK") URL: https://github.com/config-i1/smooth BugReports: https://github.com/config-i1/smooth/issues @@ -24,12 +24,13 @@ LinkingTo: Rcpp, RcppArmadillo (>= 0.8.100.0.0) Suggests: legion, numDeriv, testthat, knitr, rmarkdown, doMC, doParallel, foreach VignetteBuilder: knitr -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 Encoding: UTF-8 +ByteCompile: true NeedsCompilation: yes -Packaged: 2023-09-16 10:17:15 UTC; config +Packaged: 2024-04-01 12:18:10 UTC; config Author: Ivan Svetunkov [aut, cre] (Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK) -Maintainer: Ivan Svetunkov +Maintainer: Ivan Svetunkov Repository: CRAN -Date/Publication: 2023-09-17 04:30:02 UTC +Date/Publication: 2024-04-01 14:50:02 UTC diff --git a/MD5 b/MD5 index 559c2ea..d6bb4ef 100644 --- a/MD5 +++ b/MD5 @@ -1,19 +1,19 @@ -562a117839a15d4b76c9a72bd4199bc2 *DESCRIPTION -0bf119877847d6f1b4a492b15a8bec8b *NAMESPACE -8b800fde5a3d1e48d658e84d31624bca *NEWS -a62d3797825f0e330012fda89473841e *R/RcppExports.R -f8253fb8ef049a6f3f36b6a7819df39e *R/adam-es.R +2da1918be2b69e50feb86328eab97eab *DESCRIPTION +42b45b14e084f26970802fa88d5cc755 *NAMESPACE +b011be1b5b1ae15a33f9d3a24f3d4e96 *NEWS +c4f084103fc8811ad410766f8daf6a68 *R/RcppExports.R +8c77efb65c43d7da9489ce779ebeff26 *R/adam-es.R d60cfdfb03976caff6dc64f456a6010f *R/adam-msarima.R -9cb07e267c7fbb4ddf8d604d452bc62c *R/adam-sma.R -c0a099f465b7abefa7288243cbece238 *R/adam.R -f864c13a06b19c4c3e814a85e64a4fc3 *R/adamGeneral.R -cabd98681834a7c73496b4f06814a519 *R/autoadam.R +d3c9adb9d9f228937db1ec85f455091f *R/adam-sma.R +ab3645eef479521e300f3f11772e81b5 *R/adam.R +10a48fa3a3cf301df1ea369be10ce4ed *R/adamGeneral.R +58b00ac1350d2aab05de6543992342e6 *R/autoadam.R ea95a9439e93ea7f2c0a22d5f536d635 *R/autoces.R 84ec09fab0c7a9340ce8f6dcf187f67a *R/autogum.R 1584b3c2373f6233ef491536629162b1 *R/automsarima.R bdde374481cc3eb7d58511ad9577e8d4 *R/autossarima.R e6ca2c1c9ac18dd169459e06f84124e7 *R/ces.R -5e24a3f832ff393e00e43326c63c80fd *R/cma.R +838d9351ac2d631487532bc9a069e937 *R/cma.R dee51918bf1e69babdd6156d6b0dabb8 *R/depricator.R 178fb5186ef072c0d2eb6b943557708c *R/es.R 644b36637f40ce1332d934f1eec16add *R/gum.R @@ -21,7 +21,7 @@ dee51918bf1e69babdd6156d6b0dabb8 *R/depricator.R 041d01baa73ad0ec1dad2547505e8678 *R/iss.R bcd2007e8736b16f79ab5b685f74fa56 *R/methods.R 804079c011b151b646eb9d6e9a4051d6 *R/msarima.R -49540534fa230cd869497488b7357b88 *R/msdecompose.R +23d204616f0a6f3192c3f7bd83ab6de2 *R/msdecompose.R 279568641df34ca6f82a0c668e938ea2 *R/oes.R 63910e171affdc1a6988603a95025638 *R/oesg.R 5e17f699c228ba9941f1043825a7d604 *R/randomARIMA.R @@ -33,7 +33,7 @@ a26210beaf498241d089d307cb553e10 *R/simgum.R 2af7ae71a16d99aa03960cfb8b54bd28 *R/simssarima.R e6b80264b69b331f69c361c08720bca1 *R/sm.R 7fb01d6c2094b9b693aca573e60ff4f0 *R/sma.R -9a19d731523b676605bf6bbfa871af8b *R/smooth-package.R +2f369063fb706fc8473e00868d9e6f25 *R/smooth-package.R 1f6a8a5773a65ae3990a9e6178415ed3 *R/smoothCombine.R c2d9138e7acfbb6ce34906896b820b7e *R/sowhat.R 5a0d49de5494b56c7034175851c4535c *R/ssarima.R @@ -41,49 +41,49 @@ eb54820d225d371b94dfcfc3eb3c57bd *R/ssfunctions.R f24826d22eba6d731300f952efdc6db5 *R/variance-covariance.R ab5988fe4671de0c01accffc7485a8c5 *R/zzz.R e8405c90e8f2c8fe2ee501e851b26896 *README.md -7080893e02c49cd296d4424b9be55069 *build/partial.rdb -25a319928b6b47d06069d3656daa8ba0 *build/vignette.rds +20c8d3f97d93e2e318a22508275b5a02 *build/partial.rdb +9a683111b71fed68236000e2cf165a0c *build/vignette.rds 3474f00ebaee8436de9e09211e7130fa *inst/doc/adam.R 82c0124eeb40d5843fdd130c233b8cf1 *inst/doc/adam.Rmd -b2fbc1676210cbcb4c121d176d3f6802 *inst/doc/adam.html +c3833b503b83053020876613d9d1846f *inst/doc/adam.html b9e5f76276b38ed63ec627214c59c650 *inst/doc/ces.R 82ba0016e79b233c075aee7fd2f1e703 *inst/doc/ces.Rmd -8e8207849b01e7e86b17aeac2cfba1cb *inst/doc/ces.html +83c5b2bbf6dec26c180ea3c37b0753d1 *inst/doc/ces.html 68e24e5e7e0d42cdc6f6807d7f4d8c56 *inst/doc/es.R 0e97516b5daa5c8d147653db40be63ee *inst/doc/es.Rmd -4a5c5a57975e6b909496c013c1d64f8a *inst/doc/es.html +871adf6332d8e5d5b443299f1aad52c9 *inst/doc/es.html 2e237111beaf1e19e30fb3c019c888a1 *inst/doc/gum.R 1f1f86af77a9e7be83bacc96c7fb5bb4 *inst/doc/gum.Rmd -68e1837f718a5a51ba5003c52454e7b7 *inst/doc/gum.html +a2ec8bf1bc09cd7357c33a25ba5fdb55 *inst/doc/gum.html 64ec686041992e6ba10c2dd4bb169551 *inst/doc/oes.R 3560d7cf12287f7d6670067b7ad90cb1 *inst/doc/oes.Rmd -87a4b2d0749f9e337ba567017333ca12 *inst/doc/oes.html +7c41e0cfbbbc20858de9a57fa929f58a *inst/doc/oes.html da2f860c38092079751953639e2bd6b9 *inst/doc/simulate.R 3456eace8d7510180a0cdd3605c179ea *inst/doc/simulate.Rmd -26aee711c0ea1c62c9cbd2d1f2356ff8 *inst/doc/simulate.html +86c4eef65a22a908874e1c67187440f9 *inst/doc/simulate.html ca9229ba98cf31c27a1534eb73d9ef10 *inst/doc/sma.R 677443549ff9966b898196aabb087151 *inst/doc/sma.Rmd -687d35cc7697528ed3c340a179d7bdd3 *inst/doc/sma.html +ca5805b6e0fc601c2971b75186247404 *inst/doc/sma.html f2be0cff7be52faff06f2a377333a8c6 *inst/doc/smooth-Documentation.pdf 69802db80bcf5775ededb1dfc183e7a3 *inst/doc/smooth.R -e9a1a24ea9b95130cbef42a2085f378a *inst/doc/smooth.Rmd -238dfe228e9fea020a60f7d5319ff2df *inst/doc/smooth.html +c956ff6e4644fcfb50ab12fb44b3ce84 *inst/doc/smooth.Rmd +9da5fe94495348916d2edfa52e74aa31 *inst/doc/smooth.html d32873978eae67406eabac3b5c61cc24 *inst/doc/ssarima.R c5df97a93785052b7e8ac58652d7ded8 *inst/doc/ssarima.Rmd -cf06a51406be2f7c9c99f65ca66e3efa *inst/doc/ssarima.html +2400f47ab064cd2b8c3e12b9e0206c98 *inst/doc/ssarima.html 135d01678fd26cf81756b388d8033ae6 *man/accuracy.Rd -765c4d878947e743d0168ba7951f664a *man/adam.Rd +6b34ec5cb66585dec4517d3e98de993a *man/adam.Rd 2064c34918c99da32bb293f644d1049d *man/auto.ces.Rd 9d70d1633a8c0ebe1a8edb78c62f67b2 *man/auto.gum.Rd b1970940c2ac1c033449524a4361b5d2 *man/auto.ssarima.Rd 2981c36996860928a926e0d57297a905 *man/ces.Rd -7ebfc529b07de674f85d58f8594e8fe0 *man/cma.Rd -2a771d2ca635a63be4ecae3feac90520 *man/es.Rd +d2c1cbdd8e763c11a1e3450d84782090 *man/cma.Rd +7b65975013b9f309d8449ab544cfd760 *man/es.Rd 163b8603f9a45ef305b7046ef0c647f4 *man/forecast.smooth.Rd -eb6175f636bd9ba3a1440e2d28f364f9 *man/gum.Rd +29dae08a95088bfce1a5732f104abb2f *man/gum.Rd 7ce89f3f44bc255844499d9b3a7d28dc *man/isFunctions.Rd cc98946de57958be2d727e1659b6feb6 *man/msarima.Rd -09f1446e3dcabcc7514b8d9fcf3f3d31 *man/msdecompose.Rd +9791a179866a28bb975e37f83fca7e46 *man/msdecompose.Rd bc9beae49b356b2cd99143e1744a1b7a *man/multicov.Rd 303b3dead241b824860911ff0b8985e1 *man/oes.Rd da6f58f8aa4baa8edd9e146f86e0591e *man/oesg.Rd @@ -95,26 +95,26 @@ be7382f8b2d98bd730b7ab6606beece3 *man/reexports.Rd fed68f135199be9b43560dfc207ec219 *man/rmultistep.Rd ba9ecd1d8c6703453ee1175e933ab414 *man/sim.ces.Rd c55eba92adfb3630515f1d1d0e3cfbf8 *man/sim.es.Rd -050e101dde29069cd2ea1d11882741a0 *man/sim.gum.Rd +efba83e77f5f1143f3f6b08062b1725a *man/sim.gum.Rd e2e33ea158da93cf2577ab9c8a45fe69 *man/sim.oes.Rd c7c700f2093c421a10cfcef28828bae5 *man/sim.sma.Rd f9abd11f3a036809d9b84a093a3e55dd *man/sim.ssarima.Rd -656f3fb6b50b78941b6e1f2a1bd5523f *man/sma.Rd -5b254f0b1a5511310da603e55345ed2b *man/smooth.Rd +5f9e2800c60d8b22bb2cbb066da42c91 *man/sma.Rd +39019b140789537fbf7d753012562909 *man/smooth.Rd b66a2f843f81d19ddfb1476f7df3f234 *man/smoothCombine.Rd 937692b5fb25d64ce157f8815973ddd1 *man/sowhat.Rd 716e28b94fcebc0c4d0c19df77cff789 *man/ssarima.Rd 2e7b3f56789236798d3f6642813f224c *src/Makevars 7cc404a2b2d3a6aea6d80ef2d080e3c7 *src/Makevars.win -ce9de74824160030765ca2a360f46716 *src/RcppExports.cpp -5b06fafc4201f4c02590927142a06adc *src/adamGeneral.cpp -b36bcca51dd6ae49171ca576468a9a13 *src/adamGeneral.h -be656160bfd37c26ecdc5058a9967f5f *src/adamRefitter.cpp -85e0a292d0be7883590e02684057c002 *src/adamSimulator.cpp +8e457cc39ddbc90b7d23e03f830cb801 *src/RcppExports.cpp +fb257039cd14249b843880776cfe5eee *src/adamGeneral.cpp +d262d8c40bc6304e9b3ac2380951ad6c *src/adamGeneral.h +e679ff3ae4451d365bd29eb0a7bb45a6 *src/adamRefitter.cpp +1a6c8ac3da8750e1b06c5b343fb32876 *src/adamSimulator.cpp 3e778d3caa119ed42a745ca109cdecc3 *src/matrixPowerWrap.cpp 9fb0915a4e7ffe192150808e5e890c69 *src/registerDynamicSymbol.c -7e976d7ca79380b4fd2e517aa92d55fc *src/ssGeneral.cpp -80552072ea20f750d412a9f510a2d337 *src/ssGeneral.h +95e13cf27528965b1d8dc790f46ccea0 *src/ssGeneral.cpp +a2a1915024432000322f5f4f0b8561cf *src/ssGeneral.h 3b2327af0bf2cc52c6c29ba3c6c2adab *src/ssOccurrence.cpp 2251fdf370b5e6cfbafa86986dc99c80 *src/ssSimulator.cpp 4e0f43b23ba7abbb29b225614e88f276 *tests/testthat.R @@ -134,5 +134,5 @@ bb04e095c3fb0005695affa4afd2e923 *tests/testthat/test_oes.R 3456eace8d7510180a0cdd3605c179ea *vignettes/simulate.Rmd 677443549ff9966b898196aabb087151 *vignettes/sma.Rmd f2be0cff7be52faff06f2a377333a8c6 *vignettes/smooth-Documentation.pdf -e9a1a24ea9b95130cbef42a2085f378a *vignettes/smooth.Rmd +c956ff6e4644fcfb50ab12fb44b3ce84 *vignettes/smooth.Rmd c5df97a93785052b7e8ac58652d7ded8 *vignettes/ssarima.Rmd diff --git a/NAMESPACE b/NAMESPACE index 1a8c463..0fbe697 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -116,6 +116,7 @@ S3method(rstandard,smooth) S3method(rstudent,adam) S3method(rstudent,smooth) S3method(sigma,adam) +S3method(sigma,msdecompose) S3method(sigma,smooth) S3method(sigma,smooth.sim) S3method(simulate,adam) diff --git a/NEWS b/NEWS index 5b77663..4125106 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,22 @@ +smooth v4.0.1 (Release data: 2024-04-01) +======= + +Changes: +* Initialisation of ARIMA is now more precise, both for "optimal" and "backcasting", leading to a better fit of the model to the data. This was done using more meaningful initial values of states and by cutting the small step of refitting the head states in ARIMA with backcasting. +* Further tuning of the initialisation via backcasting. This time in C++ code for the preheating period. +* Rewrote the polynomialiser() function in C++. This should speed up ADAM ARIMA and msarima(). +* Renamed the internal variable "profileObseredTable" into a more meaningful "indexLookupTable". + +Bugfixes: +* Fixed a bug with multicov.adam() and simulated covariance matrix. +* adam() ARIMA and msarima() with backcasting would count number of components towards the number of estimated parameters. +* adam() would not remove components of ETS correctly in cases of small samples (e.g. obs == 2*frequency). +* Fixed the code of sma() - there was a bug in the architecture resulting in wrong fitting. +* ARIMA parameters are now initialised even when ACF/PACF returns NaNs (exotic case). +* Fixed bugs in es() not accepting xreg of a different length than y and not using the existing data in pre-estimated model via model=ourModel. Thanks to Nikos Kourentzes for spotting these! +* Another bugfix in es() and xreg: resolved the issue with the simulated data. + + smooth v4.0.0 (Release data: 2023-09-15) ======= diff --git a/R/RcppExports.R b/R/RcppExports.R index 51ed2d4..08978f1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,28 +1,32 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -adamFitterWrap <- function(matrixVt, matrixWt, matrixF, vectorG, lags, profilesObserved, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, vectorYt, vectorOt, backcast) { - .Call('_smooth_adamFitterWrap', PACKAGE = 'smooth', matrixVt, matrixWt, matrixF, vectorG, lags, profilesObserved, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, vectorYt, vectorOt, backcast) +adamFitterWrap <- function(matrixVt, matrixWt, matrixF, vectorG, lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, vectorYt, vectorOt, backcast) { + .Call('_smooth_adamFitterWrap', PACKAGE = 'smooth', matrixVt, matrixWt, matrixF, vectorG, lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, vectorYt, vectorOt, backcast) } -adamForecasterWrap <- function(matrixWt, matrixF, lags, profilesObserved, profilesRecent, E, T, S, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon) { - .Call('_smooth_adamForecasterWrap', PACKAGE = 'smooth', matrixWt, matrixF, lags, profilesObserved, profilesRecent, E, T, S, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon) +adamForecasterWrap <- function(matrixWt, matrixF, lags, indexLookupTable, profilesRecent, E, T, S, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon) { + .Call('_smooth_adamForecasterWrap', PACKAGE = 'smooth', matrixWt, matrixF, lags, indexLookupTable, profilesRecent, E, T, S, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon) } -adamErrorerWrap <- function(matrixVt, matrixWt, matrixF, lags, profilesObserved, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon, vectorYt, vectorOt) { - .Call('_smooth_adamErrorerWrap', PACKAGE = 'smooth', matrixVt, matrixWt, matrixF, lags, profilesObserved, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon, vectorYt, vectorOt) +adamErrorerWrap <- function(matrixVt, matrixWt, matrixF, lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon, vectorYt, vectorOt) { + .Call('_smooth_adamErrorerWrap', PACKAGE = 'smooth', matrixVt, matrixWt, matrixF, lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon, vectorYt, vectorOt) } -adamRefitterWrap <- function(matrixYt, matrixOt, arrayVt, arrayF, arrayWt, matrixG, E, T, S, lags, profilesObserved, arrayProfilesRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) { - .Call('_smooth_adamRefitterWrap', PACKAGE = 'smooth', matrixYt, matrixOt, arrayVt, arrayF, arrayWt, matrixG, E, T, S, lags, profilesObserved, arrayProfilesRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) +adamPolynomialiser <- function(B, arOrders, iOrders, maOrders, arEstimate, maEstimate, armaParameters, lags) { + .Call('_smooth_adamPolynomialiser', PACKAGE = 'smooth', B, arOrders, iOrders, maOrders, arEstimate, maEstimate, armaParameters, lags) } -adamReforecasterWrap <- function(arrayErrors, arrayOt, arrayF, arrayWt, matrixG, E, T, S, lags, profilesObserved, arrayProfileRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) { - .Call('_smooth_adamReforecasterWrap', PACKAGE = 'smooth', arrayErrors, arrayOt, arrayF, arrayWt, matrixG, E, T, S, lags, profilesObserved, arrayProfileRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) +adamRefitterWrap <- function(matrixYt, matrixOt, arrayVt, arrayF, arrayWt, matrixG, E, T, S, lags, indexLookupTable, arrayProfilesRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) { + .Call('_smooth_adamRefitterWrap', PACKAGE = 'smooth', matrixYt, matrixOt, arrayVt, arrayF, arrayWt, matrixG, E, T, S, lags, indexLookupTable, arrayProfilesRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) } -adamSimulatorWrap <- function(arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG, E, T, S, lags, profilesObserved, profilesRecent, nSeasonal, componentsNumber, nArima, nXreg, constant) { - .Call('_smooth_adamSimulatorWrap', PACKAGE = 'smooth', arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG, E, T, S, lags, profilesObserved, profilesRecent, nSeasonal, componentsNumber, nArima, nXreg, constant) +adamReforecasterWrap <- function(arrayErrors, arrayOt, arrayF, arrayWt, matrixG, E, T, S, lags, indexLookupTable, arrayProfileRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) { + .Call('_smooth_adamReforecasterWrap', PACKAGE = 'smooth', arrayErrors, arrayOt, arrayF, arrayWt, matrixG, E, T, S, lags, indexLookupTable, arrayProfileRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant) +} + +adamSimulatorWrap <- function(arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG, E, T, S, lags, indexLookupTable, profilesRecent, nSeasonal, componentsNumber, nArima, nXreg, constant) { + .Call('_smooth_adamSimulatorWrap', PACKAGE = 'smooth', arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG, E, T, S, lags, indexLookupTable, profilesRecent, nSeasonal, componentsNumber, nArima, nXreg, constant) } matrixPowerWrap <- function(matA, power) { diff --git a/R/adam-es.R b/R/adam-es.R index e3aa931..ad8d7d2 100644 --- a/R/adam-es.R +++ b/R/adam-es.R @@ -30,7 +30,7 @@ #' #' Also, there are posts about the functions of the package smooth on the #' website of Ivan Svetunkov: -#' \url{https://forecasting.svetunkov.ru/en/tag/smooth/} - they explain the +#' \url{https://openforecast.org/category/r-en/smooth/} - they explain the #' underlying models and how to use the functions. #' #' @@ -258,10 +258,11 @@ es <- function(y, model="ZZZ", lags=c(frequency(y)), persistence=NULL, phi=NULL, } # If a previous model provided as a model, write down the variables - if(is.smooth(model) | is.smooth.sim(model)){ + if(is.smooth(model) || is.smooth.sim(model)){ if(smoothType(model)!="ETS"){ stop("The provided model is not ETS.",call.=FALSE); } + # If this is the simulated data, extract the parameters if(is.smooth.sim(model) & !is.null(dim(model$data))){ warning("The provided model has several submodels. Choosing a random one.",call.=FALSE); @@ -273,31 +274,30 @@ es <- function(y, model="ZZZ", lags=c(frequency(y)), persistence=NULL, phi=NULL, else{ persistence <- model$persistence; initial <- model$initial; - initialSeason <- model$initialSeason; + # initialSeason <- model$initialSeason; if(any(model$probability!=1)){ occurrence <- "a"; } } phi <- model$phi; - if(is.null(xreg)){ - xreg <- model$xreg; - } - else{ - if(is.null(model$xreg)){ - xreg <- NULL; - } - else{ - if(ncol(xreg)!=ncol(model$xreg)){ - xreg <- xreg[,colnames(model$xreg)]; - } - } + if(is.null(xreg) && is.list(model$initial) && !is.null(model$initial$xreg)){ + # Exctract xreg + xreg <- model$data[,all.vars(model$formula)[-1]]; } - initialX <- model$initialX; + # initialX <- model$initialX; + + # if(is.adam(model)){ + # y <- model$data; + # } + # else{ + # y <- model$y; + # } + model <- modelType(model); - if(any(unlist(gregexpr("C",model))!=-1)){ - initial <- "o"; - } + # if(any(unlist(gregexpr("C",model))!=-1)){ + # initial <- "o"; + # } } else if(inherits(model,"ets")){ # Extract smoothing parameters @@ -353,8 +353,12 @@ es <- function(y, model="ZZZ", lags=c(frequency(y)), persistence=NULL, phi=NULL, # Merge y and xreg into one data frame if(!is.null(xreg) && is.numeric(y)){ - data <- cbind(y=as.data.frame(y),as.data.frame(xreg)); - data <- as.matrix(data) + if(is.matrix(xreg)){ + data <- as.matrix(cbind(y=as.data.frame(y),as.data.frame(xreg[1:length(y),]))); + } + else{ + data <- as.matrix(cbind(y=as.data.frame(y),as.data.frame(xreg[1:length(y)]))); + } data <- ts(data, start=start(y), frequency=frequency(y)); colnames(data)[1] <- "y"; # Give name to the explanatory variables if they do not have them @@ -371,22 +375,28 @@ es <- function(y, model="ZZZ", lags=c(frequency(y)), persistence=NULL, phi=NULL, data <- y; } - # Prepare initials if they are numeric - initialValue <- vector("list",(!is.null(initial))*1 +(!is.null(initialSeason))*1 +(!is.null(initialX))*1); - names(initialValue) <- c("level","seasonal","xreg")[c(!is.null(initial),!is.null(initialSeason),!is.null(initialX))]; - if(is.numeric(initial)){ - initialValue <- switch(length(initial), - "1"=list(level=initial[1]), - "2"=list(level=initial[1], - trend=initial[2])); - } - if(!is.null(initialSeason)){ - initialValue$seasonal <- initialSeason; - } - if(!is.null(initialX)){ - initialValue$xreg <- initialX; + # If the initial is provided in the old school style + if(!is.list(initial)){ + # Prepare initials if they are numeric + initialValue <- vector("list",(!is.null(initial))*1 +(!is.null(initialSeason))*1 +(!is.null(initialX))*1); + names(initialValue) <- c("level","seasonal","xreg")[c(!is.null(initial),!is.null(initialSeason),!is.null(initialX))]; + if(is.numeric(initial)){ + initialValue <- switch(length(initial), + "1"=list(level=initial[1]), + "2"=list(level=initial[1], + trend=initial[2])); + } + if(!is.null(initialSeason)){ + initialValue$seasonal <- initialSeason; + } + if(!is.null(initialX)){ + initialValue$xreg <- initialX; + } + if(length(initialValue)==1 && is.null(initialValue$level)){ + initialValue <- initial; + } } - if(length(initialValue)==1 && is.null(initialValue$level)){ + else{ initialValue <- initial; } diff --git a/R/adam-sma.R b/R/adam-sma.R index a62bc73..721d386 100644 --- a/R/adam-sma.R +++ b/R/adam-sma.R @@ -127,6 +127,12 @@ sma <- function(y, order=NULL, ic=c("AICc","AIC","BIC","BICc"), } else{ model <- ellipsis$model; + + if(inherits(y,"Mdata")){ + h <- y$h; + holdout <- TRUE; + y <- ts(c(y$x,y$xx),start=start(y$x),frequency=frequency(y$x)); + } } # If a previous model provided as a model, write down the variables @@ -190,39 +196,26 @@ sma <- function(y, order=NULL, ic=c("AICc","AIC","BIC","BICc"), ot[] <- 1; CreatorSMA <- function(order){ - lagsModelAll <- as.matrix(rep(1,order)); - lagsModelMax <- 1; - obsStates <- obsInSample+1; - - # Create ADAM profiles - profilesRecentTable <- matrix(0,length(lagsModelAll),lagsModelMax, - dimnames=list(lagsModelAll,NULL)); - # Create the matrix of observed profiles indices - profilesObservedTable <- matrix((1:order)-1,length(lagsModelAll),obsAll+lagsModelMax, - dimnames=list(lagsModelAll,NULL)); - - if(order>1){ - matF <- rbind(cbind(rep(1/order,order-1),diag(order-1)), - c(1/order,rep(0,order-1))); - matWt <- matrix(c(1,rep(0,order-1)),obsInSample,order,byrow=TRUE); - } - else{ - matF <- matrix(1,1,1); - matWt <- matrix(1,obsInSample,1); - } + lagsModelAll <- matrix(1:order, ncol=1); + lagsModelMax <- max(lagsModelAll); + obsStates <- obsInSample+lagsModelMax; + + # # Create ADAM profiles + adamProfiles <- adamProfileCreator(lagsModelAll, lagsModelMax, obsAll); + + indexLookupTable <- adamProfiles$lookup; + profilesRecentTable <- adamProfiles$recent; + profilesRecentTable[order,1:order] <- mean(yInSample[1:order]); + + # State space matrices + matF <- matrix(1/order,order,order); + matWt <- matrix(1,obsInSample,order,byrow=TRUE); vecG <- matrix(1/order,order); - matVt <- matrix(NA,order,obsStates); - matVt[1,1:order] <- rep(mean(yInSample[1:order]),order); - # if(order>1){ - # for(i in 2:order){ - # matVt[i,1:(order-i+1)] <- matVt[i-1,1:(order-i+1)+1] - - # matVt[1,1:(order-i+1)] * matF[i-1,1]; - # } - # } + matVt <- matrix(0,order,obsStates); #### Fitter and the losses calculation #### adamFitted <- adamFitterWrap(matVt, matWt, matF, vecG, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, order, xregNumber, constantRequired, yInSample, ot, TRUE); @@ -234,7 +227,6 @@ sma <- function(y, order=NULL, ic=c("AICc","AIC","BIC","BICc"), ICValue <- icFunction(logLik); return(ICValue); - # return(list(order=order,cfObjective=cfObjective,ICValue=ICValue,logLik=logLik)); } diff --git a/R/adam.R b/R/adam.R index c08d09e..c20947c 100644 --- a/R/adam.R +++ b/R/adam.R @@ -5,7 +5,7 @@ utils::globalVariables(c("adamFitted","algorithm","arEstimate","arOrders","arReq "horizon","iOrders","iRequired","initialArima","initialArimaEstimate", "initialArimaNumber","initialLevel","initialLevelEstimate","initialSeasonal", "initialSeasonalEstimate","initialTrend","initialTrendEstimate","lagsModelARIMA", - "lagsModelAll","lagsModelSeasonal","profilesObservedTable","profilesRecentTable", + "lagsModelAll","lagsModelSeasonal","indexLookupTable","profilesRecentTable", "other","otherParameterEstimate","lambda","lossFunction", "maEstimate","maOrders","maRequired","matVt","matWt","maxtime","modelIsTrendy", "nParamEstimated","persistenceLevel","persistenceLevelEstimate", @@ -77,6 +77,7 @@ utils::globalVariables(c("adamFitted","algorithm","arEstimate","arOrders","arReq #' @template ssAuthor #' @template ssKeywords #' +#' @template smoothRef #' @template ssADAMRef #' @template ssGeneralRef #' @template ssIntermittentRef @@ -748,14 +749,14 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) else{ profilesRecentTable <- adamProfiles$recent; } - profilesObservedTable <- adamProfiles$observed; + indexLookupTable <- adamProfiles$lookup; return(list(lagsModel=lagsModel,lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax, componentsNumberETS=componentsNumberETS, componentsNumberETSSeasonal=componentsNumberETSSeasonal, componentsNumberETSNonSeasonal=componentsNumberETS-componentsNumberETSSeasonal, componentsNamesETS=componentsNamesETS, obsStates=obsStates, modelIsTrendy=modelIsTrendy, modelIsSeasonal=modelIsSeasonal, - profilesObservedTable=profilesObservedTable, profilesRecentTable=profilesRecentTable)); + indexLookupTable=indexLookupTable, profilesRecentTable=profilesRecentTable)); } #### The function creates the necessary matrices based on the model and provided parameters #### @@ -873,8 +874,10 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # If the arma parameters were provided, fill in the persistence if(arimaModel && (!arEstimate && !maEstimate)){ # Call polynomial - arimaPolynomials <- polynomialiser(NULL, arOrders, iOrders, maOrders, - arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags); + # arimaPolynomials <- polynomialiser(NULL, arOrders, iOrders, maOrders, + # arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags); + arimaPolynomials <- lapply(adamPolynomialiser(0, arOrders, iOrders, maOrders, + arEstimate, maEstimate, armaParameters, lags), as.vector); # Fill in the transition matrix if(nrow(nonZeroARI)>0){ matF[componentsNumberETS+nonZeroARI[,2],componentsNumberETS+nonZeroARI[,2]] <- @@ -1161,6 +1164,20 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) if(initialArimaEstimate){ matVt[componentsNumberETS+1:componentsNumberARIMA, 1:initialArimaNumber] <- switch(Etype, "A"=0, "M"=1); + if(any(lags>1)){ + yDecomposition <- tail(msdecompose(yInSample, + lags[lags!=1], + type=switch(Etype, + "A"="additive", + "M"="multiplicative"))$seasonal,1)[[1]]; + } + else{ + yDecomposition <- switch(Etype, + "A"=mean(diff(yInSample[otLogical])), + "M"=exp(mean(diff(log(yInSample[otLogical]))))); + } + matVt[componentsNumberETS+componentsNumberARIMA, 1:initialArimaNumber] <- + rep(yDecomposition,ceiling(initialArimaNumber/max(lags)))[1:initialArimaNumber]; # rep(yInSample[1:initialArimaNumber],each=componentsNumberARIMA); # Failsafe mechanism in case the sample is too small @@ -1375,8 +1392,11 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # ARMA parameters. This goes before xreg in persistence if(arimaModel){ # Call the function returning ARI and MA polynomials - arimaPolynomials <- polynomialiser(B[j+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], arOrders, iOrders, maOrders, - arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags); + # arimaPolynomials <- polynomialiser(B[j+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], arOrders, iOrders, maOrders, + # arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags); + arimaPolynomials <- lapply(adamPolynomialiser(B[j+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + arOrders, iOrders, maOrders, + arEstimate, maEstimate, armaParameters, lags), as.vector); # Fill in the transition matrix if(nrow(nonZeroARI)>0){ @@ -1448,7 +1468,8 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) } # This is needed in order to propagate initials of ARIMA to all components else if(any(c(arEstimate,maEstimate))){ - if(nrow(nonZeroARI)>0 && nrow(nonZeroARI)>=nrow(nonZeroMA)){ + # if(nrow(nonZeroARI)>0 && nrow(nonZeroARI)>=nrow(nonZeroMA)){ + # if(nrow(nonZeroARI)>0){ matVt[componentsNumberETS+nonZeroARI[,2], 1:initialArimaNumber] <- switch(Etype, "A"= arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*% @@ -1457,7 +1478,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) "M"=exp(arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*% t(log(matVt[componentsNumberETS+componentsNumberARIMA, 1:initialArimaNumber])) / tail(arimaPolynomials$ariPolynomial,1))); - } + # } # else{ # matVt[componentsNumberETS+nonZeroMA[,2], # 1:initialArimaNumber] <- @@ -1663,7 +1684,12 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) } for(i in 1:length(lags)){ if(arRequired && arEstimate && arOrders[i]>0){ - B[j+c(1:arOrders[i])] <- pacfValues[c(1:arOrders[i])*lags[i]]; + if(all(!is.nan(pacfValues[c(1:arOrders[i])*lags[i]]))){ + B[j+c(1:arOrders[i])] <- pacfValues[c(1:arOrders[i])*lags[i]]; + } + else{ + B[j+c(1:arOrders[i])] <- 0.1; + } if(sum(B[j+c(1:arOrders[i])])>1){ B[j+c(1:arOrders[i])] <- B[j+c(1:arOrders[i])] / sum(B[j+c(1:arOrders[i])]) - 0.01; } @@ -1674,7 +1700,12 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) j[] <- j + arOrders[i]; } if(maRequired && maEstimate && maOrders[i]>0){ - B[j+c(1:maOrders[i])] <- acfValues[c(1:maOrders[i])*lags[i]]; + if(all(!is.nan(acfValues[c(1:maOrders[i])*lags[i]]))){ + B[j+c(1:maOrders[i])] <- acfValues[c(1:maOrders[i])*lags[i]]; + } + else{ + B[j+c(1:maOrders[i])] <- 0.1; + } if(sum(B[j+c(1:maOrders[i])])>1){ B[j+c(1:maOrders[i])] <- B[j+c(1:maOrders[i])] / sum(B[j+c(1:maOrders[i])]) - 0.01; } @@ -1764,6 +1795,8 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) Bu[j+1:initialArimaNumber] <- Inf; } else{ + # Make sure that ARIMA states are positive to avoid errors + B[j+1:initialArimaNumber] <- abs(B[j+1:initialArimaNumber]); Bl[j+1:initialArimaNumber] <- 0; Bu[j+1:initialArimaNumber] <- Inf; } @@ -1894,7 +1927,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal, componentsNumberARIMA, lags, lagsModel, lagsModelAll, lagsModelMax, - profilesObservedTable, profilesRecentTable, + indexLookupTable, profilesRecentTable, # The main matrices matVt, matWt, matF, vecG, # Persistence and phi @@ -2064,10 +2097,11 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) profilesRecentTable[] <- adamElements$matVt[,1:lagsModelMax]; # print(round(B,3)) # print(adamElements$vecG) + # print(profilesRecentTable) #### Fitter and the losses calculation #### adamFitted <- adamFitterWrap(adamElements$matVt, adamElements$matWt, adamElements$matF, adamElements$vecG, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, yInSample, ot, any(initialType==c("complete","backcasting"))); @@ -2235,7 +2269,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) else{ # Call for the Rcpp function to produce a matrix of multistep errors adamErrors <- adamErrorerWrap(adamFitted$matVt, adamElements$matWt, adamElements$matF, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, h, @@ -2274,7 +2308,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal, componentsNumberARIMA, lags, lagsModel, lagsModelAll, lagsModelMax, - profilesObservedTable, profilesRecentTable, + indexLookupTable, profilesRecentTable, matVt, matWt, matF, vecG, persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate, persistenceSeasonalEstimate, persistenceXregEstimate, @@ -2368,7 +2402,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal, componentsNumberARIMA, lags, lagsModel, lagsModelAll, lagsModelMax, - profilesObservedTable, profilesRecentTable, + indexLookupTable, profilesRecentTable, matVt, matWt, matF, vecG, persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate, persistenceSeasonalEstimate, persistenceXregEstimate, @@ -2427,7 +2461,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal, componentsNumberARIMA, lags, lagsModel, lagsModelAll, lagsModelMax, - profilesObservedTable, profilesRecentTable, + indexLookupTable, profilesRecentTable, matVt, matWt, matF, vecG, persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate, persistenceSeasonalEstimate, persistenceXregEstimate, @@ -2489,7 +2523,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # Fit the model again to extract the fitted values adamFitted <- adamFitterWrap(adamElements$matVt, adamElements$matWt, adamElements$matF, adamElements$vecG, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, yInSample, ot, any(initialType==c("complete","backcasting"))); @@ -2581,7 +2615,6 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # print(BValues$B); # Preheat the initial state of ARIMA. Do this only for optimal initials and if B is not provided - # This is also not needed for I(d) model and d>1, as the backcasting hurts in this case if(arimaModel && initialType=="optimal" && initialArimaEstimate && is.null(B)){ adamCreatedARIMA <- filler(BValues$B, etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, @@ -2607,7 +2640,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # Do initial fit to get the state values from the backcasting adamFitted <- adamFitterWrap(adamCreatedARIMA$matVt, adamCreatedARIMA$matWt, adamCreatedARIMA$matF, adamCreatedARIMA$vecG, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, yInSample, ot, TRUE); @@ -2749,7 +2782,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal, componentsNumberARIMA=componentsNumberARIMA, lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax, - profilesObservedTable=profilesObservedTable, profilesRecentTable=profilesRecentTable, + indexLookupTable=indexLookupTable, profilesRecentTable=profilesRecentTable, matVt=adamCreated$matVt, matWt=adamCreated$matWt, matF=adamCreated$matF, vecG=adamCreated$vecG, persistenceEstimate=persistenceEstimate, persistenceLevelEstimate=persistenceLevelEstimate, @@ -2799,7 +2832,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal, componentsNumberARIMA=componentsNumberARIMA, lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax, - profilesObservedTable=profilesObservedTable, profilesRecentTable=profilesRecentTable, + indexLookupTable=indexLookupTable, profilesRecentTable=profilesRecentTable, matVt=adamCreated$matVt, matWt=adamCreated$matWt, matF=adamCreated$matF, vecG=adamCreated$vecG, persistenceEstimate=persistenceEstimate, @@ -2853,7 +2886,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal, componentsNumberARIMA, lags, lagsModel, lagsModelAll, lagsModelMax, - profilesObservedTable, profilesRecentTable, + indexLookupTable, profilesRecentTable, adamCreated$matVt, adamCreated$matWt, adamCreated$matF, adamCreated$vecG, persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate, persistenceSeasonalEstimate, persistenceXregEstimate, @@ -2924,7 +2957,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # Fit the model to the data adamFitted <- adamFitterWrap(adamCreated$matVt, adamCreated$matWt, adamCreated$matF, adamCreated$vecG, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, yInSample, ot, any(initialType==c("complete","backcasting"))); @@ -3548,7 +3581,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) # Fit the model to the data adamFitted <- adamFitterWrap(matVt, matWt, matF, vecG, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, yInSample, ot, any(initialType==c("complete","backcasting"))); @@ -3607,7 +3640,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) yForecast[] <- adamForecasterWrap(tail(matWt,horizon), matF, lagsModelAll, - profilesObservedTable[,lagsModelMax+obsInSample+c(1:horizon),drop=FALSE], + indexLookupTable[,lagsModelMax+obsInSample+c(1:horizon),drop=FALSE], profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, @@ -4209,7 +4242,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) adamSelected$results[[i]]$lagsModelAll <- adamArchitect$lagsModelAll; adamSelected$results[[i]]$lagsModelMax <- adamArchitect$lagsModelMax; adamSelected$results[[i]]$profilesRecentTable <- adamArchitect$profilesRecentTable; - adamSelected$results[[i]]$profilesObservedTable <- adamArchitect$profilesObservedTable; + adamSelected$results[[i]]$indexLookupTable <- adamArchitect$indexLookupTable; adamSelected$results[[i]]$componentsNumberETS <- adamArchitect$componentsNumberETS; adamSelected$results[[i]]$componentsNumberETSSeasonal <- adamArchitect$componentsNumberETSSeasonal; adamSelected$results[[i]]$componentsNumberETSNonSeasonal <- adamArchitect$componentsNumberETSNonSeasonal; @@ -4321,7 +4354,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal, componentsNumberARIMA=componentsNumberARIMA, lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax, - profilesObservedTable=profilesObservedTable, profilesRecentTable=profilesRecentTable, + indexLookupTable=indexLookupTable, profilesRecentTable=profilesRecentTable, matVt=matVt, matWt=matWt, matF=matF, vecG=vecG, persistenceEstimate=persistenceEstimate, persistenceLevelEstimate=persistenceLevelEstimate, @@ -4356,7 +4389,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal, componentsNumberARIMA, lags, lagsModel, lagsModelAll, lagsModelMax, - profilesObservedTable, profilesRecentTable, + indexLookupTable, profilesRecentTable, matVt, matWt, matF, vecG, persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate, persistenceSeasonalEstimate, persistenceXregEstimate, @@ -4482,7 +4515,7 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal, componentsNumberARIMA=componentsNumberARIMA, lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax, - profilesObservedTable=profilesObservedTable, profilesRecentTable=profilesRecentTable, + indexLookupTable=indexLookupTable, profilesRecentTable=profilesRecentTable, matVt=matVt, matWt=matWt, matF=matF, vecG=vecG, persistenceEstimate=persistenceEstimateFI, persistenceLevelEstimate=persistenceLevelEstimateFI, persistenceTrendEstimate=persistenceTrendEstimateFI, @@ -4907,7 +4940,16 @@ adam <- function(data, model="ZXZ", lags=c(frequency(data)), orders=list(ar=c(0) } #### Small useful ADAM functions #### -# This function creates recent and observed profiles for adam +# These functions are faster than which() and tail() for vectors are. +# The main gain is in polinomialiser() +# whichFast <- function(x){ +# return(c(1:length(x))[x]); +# } +# tailFast <- function(x,...){ +# return(x[length(x)]) ; +# } + +# This function creates recent profile and the lookup table for adam #' @importFrom greybox detectdst adamProfileCreator <- function(lagsModelAll, lagsModelMax, obsAll, lags=NULL, yIndex=NULL, yClasses=NULL){ @@ -4921,34 +4963,34 @@ adamProfileCreator <- function(lagsModelAll, lagsModelMax, obsAll, # Create the matrix with profiles, based on provided lags profilesRecentTable <- matrix(0,length(lagsModelAll),lagsModelMax, dimnames=list(lagsModelAll,NULL)); - # Create the matrix of observed profiles indices - profilesObservedTable <- matrix(1,length(lagsModelAll),obsAll+lagsModelMax, + # Create the lookup table + indexLookupTable <- matrix(1,length(lagsModelAll),obsAll+lagsModelMax, dimnames=list(lagsModelAll,NULL)); - # Modify the observed ones in order to get proper indices in C++ + # Modify the lookup table in order to get proper indices in C++ profileIndices <- matrix(c(1:(lagsModelMax*length(lagsModelAll))),length(lagsModelAll)); for(i in 1:length(lagsModelAll)){ profilesRecentTable[i,1:lagsModelAll[i]] <- 1:lagsModelAll[i]; # -1 is needed to align this with C++ code - profilesObservedTable[i,lagsModelMax+c(1:obsAll)] <- rep(profileIndices[i,1:lagsModelAll[i]], + indexLookupTable[i,lagsModelMax+c(1:obsAll)] <- rep(profileIndices[i,1:lagsModelAll[i]], ceiling(obsAll/lagsModelAll[i]))[1:obsAll] -1; # Fix the head of the data, before the sample starts - profilesObservedTable[i,1:lagsModelMax] <- tail(rep(unique(profilesObservedTable[i,lagsModelMax+c(1:obsAll)]),lagsModelMax), + indexLookupTable[i,1:lagsModelMax] <- tail(rep(unique(indexLookupTable[i,lagsModelMax+c(1:obsAll)]),lagsModelMax), lagsModelMax); } # Do shifts for proper lags only: # Check lags variable for 24 / 24*7 / 24*365 / 48 / 48*7 / 48*365 / 365 / 52 # If they are there, find the DST / Leap moments - # Then amend respective observed values of profile, shifting them around + # Then amend respective lookup values of profile, shifting them around if(any(yClasses=="zoo") && !is.null(yIndex) && !is.numeric(yIndex)){ # If this is weekly data, duplicate 52, when 53 is used if(any(lags==52) && any(strftime(yIndex,format="%W")=="53")){ shiftRows <- lagsModelAll==52; # If the data does not start with 1, proceed if(all(which(strftime(yIndex,format="%W")=="53")!=1)){ - profilesObservedTable[shiftRows,which(strftime(yIndex,format="%W")=="53")] <- - profilesObservedTable[shiftRows,which(strftime(yIndex,format="%W")=="53")-1]; + indexLookupTable[shiftRows,which(strftime(yIndex,format="%W")=="53")] <- + indexLookupTable[shiftRows,which(strftime(yIndex,format="%W")=="53")-1]; } } @@ -4958,8 +5000,8 @@ adamProfileCreator <- function(lagsModelAll, lagsModelMax, obsAll, shiftRows <- lagsModelAll %in% c(365,365*48,365*24); # If the data does not start with 1/24/48, proceed (otherwise we refer to negative numbers) if(!any(which(strftime(yIndex,format="%d/%m")=="29/02") %in% shiftValue)){ - profilesObservedTable[shiftRows,which(strftime(yIndex,format="%d/%m")=="29/02")] <- - profilesObservedTable[shiftRows,which(strftime(yIndex,format="%d/%m")=="29/02")-shiftValue]; + indexLookupTable[shiftRows,which(strftime(yIndex,format="%d/%m")=="29/02")] <- + indexLookupTable[shiftRows,which(strftime(yIndex,format="%d/%m")=="29/02")-shiftValue]; } } @@ -4991,119 +5033,119 @@ adamProfileCreator <- function(lagsModelAll, lagsModelMax, obsAll, for(i in 1:nrow(dstValues$start)){ # If the end date is natural, just shift if(dstValues$end$id[i]+shiftValue<=obsAll){ - profilesObservedTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]] <- - profilesObservedTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]+shiftValue]; + indexLookupTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]] <- + indexLookupTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]+shiftValue]; } # If it isn't, we need to come up with the values for the end of sample else{ - profilesObservedTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]] <- - profilesObservedTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]-lagsModelMax+shiftValue]; + indexLookupTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]] <- + indexLookupTable[shiftRows,dstValues$start$id[i]:dstValues$end$id[i]-lagsModelMax+shiftValue]; } } } } } - return(list(recent=profilesRecentTable,observed=profilesObservedTable)); + return(list(recent=profilesRecentTable,lookup=indexLookupTable)); } #### ARI and MA polynomials function #### -polynomialiser <- function(B, arOrders, iOrders, maOrders, - arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags){ - - # Number of parameters that we have - nParamAR <- sum(arOrders); - nParamMA <- sum(maOrders); - - # Matrices with parameters - arParameters <- matrix(0, max(arOrders * lags) + 1, length(arOrders)); - iParameters <- matrix(0, max(iOrders * lags) + 1, length(iOrders)); - maParameters <- matrix(0, max(maOrders * lags) + 1, length(maOrders)); - # The first element is always 1 - arParameters[1,] <- iParameters[1,] <- maParameters[1,] <- 1; - - # nParam is used for B - nParam <- 1; - # armanParam is used for the provided arma parameters - armanParam <- 1; - # Fill in the matrices with the provided parameters - for(i in 1:length(lags)){ - if(arOrders[i]*lags[i]!=0){ - if(arEstimate){ - arParameters[1+(1:arOrders[i])*lags[i],i] <- -B[nParam+c(1:arOrders[i])-1]; - nParam <- nParam + arOrders[i]; - } - else if(!arEstimate && arRequired){ - arParameters[1+(1:arOrders[i])*lags[i],i] <- -armaParameters[armanParam+c(1:arOrders[i])-1]; - armanParam <- armanParam + arOrders[i]; - } - } - - if(iOrders[i]*lags[i] != 0){ - iParameters[1+lags[i],i] <- -1; - } - - if(maOrders[i]*lags[i]!=0){ - if(maEstimate){ - maParameters[1+(1:maOrders[i])*lags[i],i] <- B[nParam+c(1:maOrders[i])-1]; - nParam <- nParam + maOrders[i]; - } - else if(!maEstimate && maRequired){ - maParameters[1+(1:maOrders[i])*lags[i],i] <- armaParameters[armanParam+c(1:maOrders[i])-1]; - armanParam <- armanParam + maOrders[i]; - } - } - } - - # Vectors of polynomials for the ARIMA - arPolynomial <- vector("numeric", sum(arOrders * lags) + 1); - iPolynomial <- vector("numeric", sum(iOrders * lags) + 1); - maPolynomial <- vector("numeric", sum(maOrders * lags) + 1); - ariPolynomial <- vector("numeric", sum(arOrders * lags) + sum(iOrders * lags) + 1); - - # Fill in the first polynomials - arPolynomial[0:(arOrders[1]*lags[1])+1] <- arParameters[0:(arOrders[1]*lags[1])+1,1]; - iPolynomial[0:(iOrders[1]*lags[1])+1] <- iParameters[0:(iOrders[1]*lags[1])+1,1]; - maPolynomial[0:(maOrders[1]*lags[1])+1] <- maParameters[0:(maOrders[1]*lags[1])+1,1]; - - index1 <- 0 - index2 <- 0; - # Fill in all the other polynomials - for(i in 1:length(lags)){ - if(i!=1){ - if(arOrders[i]>0){ - index1[] <- tail(which(arPolynomial!=0),1); - index2[] <- tail(which(arParameters[,i]!=0),1); - arPolynomial[1:(index1+index2-1)] <- polyprod(arPolynomial[1:index1], arParameters[1:index2,i]); - } - - if(maOrders[i]>0){ - index1[] <- tail(which(maPolynomial!=0),1); - index2[] <- tail(which(maParameters[,i]!=0),1); - maPolynomial[1:(index1+index2-1)] <- polyprod(maPolynomial[1:index1], maParameters[1:index2,i]); - } - - if(iOrders[i]>0){ - index1[] <- tail(which(iPolynomial!=0),1); - index2[] <- tail(which(iParameters[,i]!=0),1); - iPolynomial[1:(index1+index2-1)] <- polyprod(iPolynomial[1:index1], iParameters[1:index2,i]); - } - } - # This part takes the power of (1-B)^D - if(iOrders[i]>1){ - for(j in 2:iOrders[i]){ - index1[] <- tail(which(iPolynomial!=0),1); - index2[] <- tail(which(iParameters[,i]!=0),1); - iPolynomial[1:(index1+index2-1)] = polyprod(iPolynomial[1:index1], iParameters[1:index2,i]); - } - } - } - # ARI polynomials - ariPolynomial[] <- polyprod(arPolynomial, iPolynomial); - - return(list(arPolynomial=arPolynomial,iPolynomial=iPolynomial, - ariPolynomial=ariPolynomial,maPolynomial=maPolynomial)); -} +# polynomialiser <- function(B, arOrders, iOrders, maOrders, +# arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags){ +# +# # Number of parameters that we have +# nParamAR <- sum(arOrders); +# nParamMA <- sum(maOrders); +# +# # Matrices with parameters +# arParameters <- matrix(0, max(arOrders * lags) + 1, length(arOrders)); +# iParameters <- matrix(0, max(iOrders * lags) + 1, length(iOrders)); +# maParameters <- matrix(0, max(maOrders * lags) + 1, length(maOrders)); +# # The first element is always 1 +# arParameters[1,] <- iParameters[1,] <- maParameters[1,] <- 1; +# +# # nParam is used for B +# nParam <- 1; +# # armanParam is used for the provided arma parameters +# armanParam <- 1; +# # Fill in the matrices with the provided parameters +# for(i in 1:length(lags)){ +# if(arOrders[i]*lags[i]!=0){ +# if(arEstimate){ +# arParameters[1+(1:arOrders[i])*lags[i],i] <- -B[nParam+c(1:arOrders[i])-1]; +# nParam[] <- nParam + arOrders[i]; +# } +# else if(!arEstimate && arRequired){ +# arParameters[1+(1:arOrders[i])*lags[i],i] <- -armaParameters[armanParam+c(1:arOrders[i])-1]; +# armanParam[] <- armanParam + arOrders[i]; +# } +# } +# +# if(iOrders[i]*lags[i] != 0){ +# iParameters[1+lags[i],i] <- -1; +# } +# +# if(maOrders[i]*lags[i]!=0){ +# if(maEstimate){ +# maParameters[1+(1:maOrders[i])*lags[i],i] <- B[nParam+c(1:maOrders[i])-1]; +# nParam[] <- nParam + maOrders[i]; +# } +# else if(!maEstimate && maRequired){ +# maParameters[1+(1:maOrders[i])*lags[i],i] <- armaParameters[armanParam+c(1:maOrders[i])-1]; +# armanParam[] <- armanParam + maOrders[i]; +# } +# } +# } +# +# # Vectors of polynomials for the ARIMA +# arPolynomial <- vector("numeric", sum(arOrders * lags) + 1); +# iPolynomial <- vector("numeric", sum(iOrders * lags) + 1); +# maPolynomial <- vector("numeric", sum(maOrders * lags) + 1); +# ariPolynomial <- vector("numeric", sum(arOrders * lags) + sum(iOrders * lags) + 1); +# +# # Fill in the first polynomials +# arPolynomial[0:(arOrders[1]*lags[1])+1] <- arParameters[0:(arOrders[1]*lags[1])+1,1]; +# iPolynomial[0:(iOrders[1]*lags[1])+1] <- iParameters[0:(iOrders[1]*lags[1])+1,1]; +# maPolynomial[0:(maOrders[1]*lags[1])+1] <- maParameters[0:(maOrders[1]*lags[1])+1,1]; +# +# index1 <- 0; +# index2 <- 0; +# # Fill in all the other polynomials +# for(i in 1:length(lags)){ +# if(i!=1){ +# if(arOrders[i]>0){ +# index1[] <- tailFast(whichFast(arPolynomial!=0)); +# index2[] <- tailFast(whichFast(arParameters[,i]!=0)); +# arPolynomial[1:(index1+index2-1)] <- polyprod(arPolynomial[1:index1], arParameters[1:index2,i]); +# } +# +# if(maOrders[i]>0){ +# index1[] <- tailFast(whichFast(maPolynomial!=0)); +# index2[] <- tailFast(whichFast(maParameters[,i]!=0)); +# maPolynomial[1:(index1+index2-1)] <- polyprod(maPolynomial[1:index1], maParameters[1:index2,i]); +# } +# +# if(iOrders[i]>0){ +# index1[] <- tailFast(whichFast(iPolynomial!=0)); +# index2[] <- tailFast(whichFast(iParameters[,i]!=0)); +# iPolynomial[1:(index1+index2-1)] <- polyprod(iPolynomial[1:index1], iParameters[1:index2,i]); +# } +# } +# # This part takes the power of (1-B)^D +# if(iOrders[i]>1){ +# for(j in 2:iOrders[i]){ +# index1[] <- tailFast(whichFast(iPolynomial!=0)); +# index2[] <- tailFast(whichFast(iParameters[,i]!=0)); +# iPolynomial[1:(index1+index2-1)] = polyprod(iPolynomial[1:index1], iParameters[1:index2,i]); +# } +# } +# } +# # ARI polynomials +# ariPolynomial[] <- polyprod(arPolynomial, iPolynomial); +# +# return(list(arPolynomial=arPolynomial,iPolynomial=iPolynomial, +# ariPolynomial=ariPolynomial,maPolynomial=maPolynomial)); +# } #### Technical methods #### #' @export @@ -5755,7 +5797,7 @@ plot.adam <- function(x, which=c(1,2,4,6), level=0.95, legend=FALSE, x$states <- cbind(actuals(x),x$states,residuals(x)); colnames(x$states) <- statesNames; if(ncol(x$states)>10){ - message("Too many states. Plotting them one by one on several graphs."); + message("Too many states. Plotting them one by one on several plots."); if(is.null(ellipsis$main)){ ellipsisMain <- NULL; } @@ -6613,7 +6655,8 @@ xtable.summary.adam <- function(x, caption = NULL, label = NULL, align = NULL, d #' @importFrom greybox coefbootstrap #' @export -coefbootstrap.adam <- function(object, nsim=100, size=floor(0.5*nobs(object)), +coefbootstrap.adam <- function(object, nsim=100, + size=floor(0.5*nobs(object)), replace=FALSE, prob=NULL, parallel=FALSE, ...){ startTime <- Sys.time(); @@ -7017,7 +7060,7 @@ rmultistep.adam <- function(object, h=10, ...){ adamProfiles <- adamProfileCreator(lagsModelAll, lagsModelMax, obsInSample, lagsOriginal, time(actuals(object)), yClasses); profilesRecentTable <- adamProfiles$recent; - profilesObservedTable <- adamProfiles$observed; + indexLookupTable <- adamProfiles$lookup; # Fill in the profile. This is done in Errorer as well, but this is just in case profilesRecentTable[] <- t(object$states[1:lagsModelMax,,drop=FALSE]); @@ -7025,7 +7068,7 @@ rmultistep.adam <- function(object, h=10, ...){ # Return multi-step errors matrix if(any(yClasses=="ts")){ return(ts(adamErrorerWrap(t(object$states), object$measurement, object$transition, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, h, @@ -7034,7 +7077,7 @@ rmultistep.adam <- function(object, h=10, ...){ } else{ return(zoo(adamErrorerWrap(t(object$states), object$measurement, object$transition, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, h, @@ -7657,10 +7700,10 @@ forecast.adam <- function(object, h=10, newdata=NULL, occurrence=NULL, yForecastIndex <- yIndex[obsInSample]+diff(tail(yIndex,2))*c(1:h); } - # Get the observed profiles - profilesObservedTable <- adamProfileCreator(lagsModelAll, lagsModelMax, obsInSample+h, + # Get the lookup table + indexLookupTable <- adamProfileCreator(lagsModelAll, lagsModelMax, obsInSample+h, lags(object), c(yIndex,yForecastIndex), - yClasses)$observed[,-c(1:(obsInSample+lagsModelMax)),drop=FALSE]; + yClasses)$lookup[,-c(1:(obsInSample+lagsModelMax)),drop=FALSE]; # All the important matrices matVt <- t(object$states[obsStates-(lagsModelMax:1)+1,,drop=FALSE]); @@ -7813,7 +7856,7 @@ forecast.adam <- function(object, h=10, newdata=NULL, occurrence=NULL, if(Ttype!="M" && (Stype!="M" | (Stype=="M" & h<=lagsModelMin)) || any(interval==c("nonparametric","semiparametric","empirical","approximate"))){ adamForecast <- adamForecasterWrap(matWt, matF, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal, componentsNumberARIMA, xregNumber, constantRequired, @@ -8016,7 +8059,7 @@ forecast.adam <- function(object, h=10, newdata=NULL, occurrence=NULL, array(matF,c(dim(matF),nsim)), matWt, matrix(vecG, componentsNumberETS+componentsNumberARIMA+xregNumber+constantRequired, nsim), EtypeModified, Ttype, Stype, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, componentsNumberETSSeasonal, componentsNumberETS, componentsNumberARIMA, xregNumber, constantRequired)$matrixYt; @@ -8916,7 +8959,7 @@ reapply.adam <- function(object, nsim=1000, bootstrap=FALSE, heuristics=NULL, .. xregParametersEstimated <- 0; xregParametersPersistence <- 0; } - profilesObservedTable <- adamProfileCreator(lagsModelAll, lagsModelMax, obsInSample)$observed; + indexLookupTable <- adamProfileCreator(lagsModelAll, lagsModelMax, obsInSample)$lookup; # Generate the data from the multivariate normal randomParameters <- mvrnorm(nsim, coef(object), vcovAdam); @@ -9180,9 +9223,12 @@ reapply.adam <- function(object, nsim=1000, bootstrap=FALSE, heuristics=NULL, .. for(i in 1:nsim){ # Call the function returning ARI and MA polynomials - arimaPolynomials <- polynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], - arOrders, iOrders, maOrders, arRequired, maRequired, arEstimate, maEstimate, - armaParameters, lags); + # arimaPolynomials <- polynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + # arOrders, iOrders, maOrders, arRequired, maRequired, arEstimate, maEstimate, + # armaParameters, lags); + arimaPolynomials <- lapply(adamPolynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + arOrders, iOrders, maOrders, + arEstimate, maEstimate, armaParameters, lags), as.vector) # Fill in the transition and persistence matrices if(nrow(nonZeroARI)>0){ @@ -9250,9 +9296,12 @@ reapply.adam <- function(object, nsim=1000, bootstrap=FALSE, heuristics=NULL, .. # Call the function returning ARI and MA polynomials ### This is not optimal, as the polynomialiser() is called twice (for parameters and here), ### but this is simpler - arimaPolynomials <- polynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], - arOrders, iOrders, maOrders, arRequired, maRequired, arEstimate, maEstimate, - armaParameters, lags); + # arimaPolynomials <- polynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + # arOrders, iOrders, maOrders, arRequired, maRequired, arEstimate, maEstimate, + # armaParameters, lags); + arimaPolynomials <- lapply(adamPolynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + arOrders, iOrders, maOrders, + arEstimate, maEstimate, armaParameters, lags), as.vector) profilesRecentArray[j+componentsNumberARIMA, 1:initialArimaNumber, i] <- randomParameters[i, k+1:initialArimaNumber]; profilesRecentArray[j+nonZeroARI[,2], 1:initialArimaNumber, i] <- @@ -9270,9 +9319,12 @@ reapply.adam <- function(object, nsim=1000, bootstrap=FALSE, heuristics=NULL, .. else{ for(i in 1:nsim){ # Call the function returning ARI and MA polynomials - arimaPolynomials <- polynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], - arOrders, iOrders, maOrders, arRequired, maRequired, arEstimate, maEstimate, - armaParameters, lags); + # arimaPolynomials <- polynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + # arOrders, iOrders, maOrders, arRequired, maRequired, arEstimate, maEstimate, + # armaParameters, lags); + arimaPolynomials <- lapply(adamPolynomialiser(randomParameters[i,polyIndex+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], + arOrders, iOrders, maOrders, + arEstimate, maEstimate, armaParameters, lags), as.vector) profilesRecentArray[componentsNumberETS+componentsNumberARIMA, 1:initialArimaNumber, i] <- randomParameters[i, k+1:initialArimaNumber]; profilesRecentArray[j+nonZeroMA[,2], 1:initialArimaNumber, i] <- @@ -9320,7 +9372,7 @@ reapply.adam <- function(object, nsim=1000, bootstrap=FALSE, heuristics=NULL, .. # Refit the model with the new parameter adamRefitted <- adamRefitterWrap(yt, ot, arrVt, arrF, arrWt, matG, Etype, Ttype, Stype, - lagsModelAll, profilesObservedTable, profilesRecentArray, + lagsModelAll, indexLookupTable, profilesRecentArray, componentsNumberETSSeasonal, componentsNumberETS, componentsNumberARIMA, xregNumber, constantRequired); arrVt[] <- adamRefitted$states; @@ -9546,8 +9598,8 @@ reforecast.adam <- function(object, h=10, newdata=NULL, occurrence=NULL, obsStates <- nrow(object$states); obsInSample <- nobs(object); - profilesObservedTable <- adamProfileCreator(lagsModelAll, lagsModelMax, - obsInSample+h)$observed[,-c(1:(obsInSample+lagsModelMax)),drop=FALSE]; + indexLookupTable <- adamProfileCreator(lagsModelAll, lagsModelMax, + obsInSample+h)$lookup[,-c(1:(obsInSample+lagsModelMax)),drop=FALSE]; yClasses <- class(actuals(object)); @@ -9847,7 +9899,7 @@ reforecast.adam <- function(object, h=10, newdata=NULL, occurrence=NULL, arrWt, objectRefitted$persistence, EtypeModified, Ttype, Stype, - lagsModelAll, profilesObservedTable, profilesRecentArray, + lagsModelAll, indexLookupTable, profilesRecentArray, componentsNumberETSSeasonal, componentsNumberETS, componentsNumberARIMA, xregNumber, constantRequired)$matrixYt; @@ -10007,9 +10059,9 @@ multicov.adam <- function(object, type=c("analytical","empirical","simulated"), Etype <- errorType(object); Stype <- substr(modelType(object),nchar(modelType(object)),nchar(modelType(object))); - # Get the observed profiles - profilesObservedTable <- adamProfileCreator(lagsModelAll, lagsModelMax, - obsInSample+h)$observed[,-c(1:(obsInSample+lagsModelMax)),drop=FALSE]; + # Get the lookup table + indexLookupTable <- adamProfileCreator(lagsModelAll, lagsModelMax, + obsInSample+h)$lookup[,-c(1:(obsInSample+lagsModelMax)),drop=FALSE]; profilesRecentTable <- object$profile; lagsModelMin <- lagsModelAll[lagsModelAll!=1]; @@ -10107,7 +10159,7 @@ multicov.adam <- function(object, type=c("analytical","empirical","simulated"), array(matF,c(dim(matF),nsim)), matWt, matrix(vecG, componentsNumberETS+componentsNumberARIMA+xregNumber+constantRequired, nsim), EtypeModified, Ttype, Stype, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, componentsNumberETSSeasonal, componentsNumberETS, componentsNumberARIMA, xregNumber, constantRequired)$matrixYt; @@ -10121,6 +10173,10 @@ multicov.adam <- function(object, type=c("analytical","empirical","simulated"), yForecast[i] <- mean(ySimulated[i,],na.rm=TRUE); } ySimulated[i,] <- ySimulated[i,]-yForecast[i]; + # If it is the multiplicative error, return epsilon_t + if(Etype=="M"){ + ySimulated[i,] <- ySimulated[i,]/yForecast[i]; + } } covarMat <- (ySimulated %*% t(ySimulated))/nsim; @@ -10412,7 +10468,7 @@ simulate.adam <- function(object, nsim=1, seed=NULL, obs=nobs(object), ...){ xregParametersPersistence <- 0; } profiles <- adamProfileCreator(lagsModelAll, lagsModelMax, obsInSample); - profilesObservedTable <- profiles$observed; + indexLookupTable <- profiles$lookup; profilesRecentTable <- profiles$recent; #### Prepare the necessary matrices #### @@ -10494,7 +10550,7 @@ simulate.adam <- function(object, nsim=1, seed=NULL, obs=nobs(object), ...){ matrix(rbinom(obsInSample*nsim, 1, pt), obsInSample, nsim), arrF, matWt, matG, EtypeModified, Ttype, Stype, - lagsModelAll, profilesObservedTable, profilesRecentTable, + lagsModelAll, indexLookupTable, profilesRecentTable, componentsNumberETSSeasonal, componentsNumberETS, componentsNumberARIMA, xregNumber, constantRequired); diff --git a/R/adamGeneral.R b/R/adamGeneral.R index 8d63f2e..5e38ec2 100644 --- a/R/adamGeneral.R +++ b/R/adamGeneral.R @@ -2497,7 +2497,7 @@ parametersChecker <- function(data, model, lags, formulaToUse, orders, constant= } } - if(arimaModel && obsNonzero < initialArimaNumber && !select){ + if(arimaModel && obsNonzero < (initialType=="optimal")*initialArimaNumber && !select){ warning(paste0("In-sample size is ",obsNonzero,", while number of ARIMA components is ",initialArimaNumber, ". Cannot fit the model."),call.=FALSE) stop("Not enough observations for such a complicated model.",call.=FALSE); @@ -2644,27 +2644,36 @@ parametersChecker <- function(data, model, lags, formulaToUse, orders, constant= else if(obsNonzero > (3 + nParamExo) && any(modelDo==c("estimate","use"))){ # We don't have enough observations for seasonal models with damped trend if((obsNonzero <= (6 + lagsModelMax + 1 + nParamExo))){ - model <- model[!(nchar(model)==4 & - substr(model,nchar(model),nchar(model))=="A")]; - model <- model[!(nchar(model)==4 & - substr(model,nchar(model),nchar(model))=="M")]; + if(nchar(model)==4){ + model <- paste0(substr(model,1,2),substr(model,4,4)); + } + # model <- model[!(nchar(model)==4 & + # substr(model,nchar(model),nchar(model))=="A")]; + # model <- model[!(nchar(model)==4 & + # substr(model,nchar(model),nchar(model))=="M")]; } # We don't have enough observations for seasonal models with trend if((obsNonzero <= (5 + lagsModelMax + 1 + nParamExo))){ - model <- model[!(substr(model,2,2)!="N" & - substr(model,nchar(model),nchar(model))!="N")]; + model <- paste0(substr(model,1,1),"N",substr(model,3,3)); + # model <- model[!(substr(model,2,2)!="N" & + # substr(model,nchar(model),nchar(model))!="N")]; } # We don't have enough observations for seasonal models if(obsNonzero <= 2*lagsModelMax){ - model <- model[substr(model,nchar(model),nchar(model))=="N"]; + model <- paste0(substr(model,1,2),"N"); + # model <- model[substr(model,nchar(model),nchar(model))=="N"]; } # We don't have enough observations for damped trend if(obsNonzero <= (6 + nParamExo)){ - model <- model[nchar(model)!=4]; + if(nchar(model)==4){ + model <- paste0(substr(model,1,2),substr(model,4,4)); + } + # model <- model[nchar(model)!=4]; } # We don't have enough observations for any trend if(obsNonzero <= (5 + nParamExo)){ - model <- model[substr(model,2,2)=="N"]; + model <- paste0(substr(model,1,1),"N",substr(model,3,3)); + # model <- model[substr(model,2,2)=="N"]; } } # Extreme cases of small samples @@ -2790,7 +2799,7 @@ parametersChecker <- function(data, model, lags, formulaToUse, orders, constant= "take more time to converge to the optimum. Consider either setting maxeval parameter ", "to a higher value (e.g. maxeval=10000, which will take ~25 times more than this) ", "or using initial='backcasting'."), - call.=FALSE); + call.=FALSE, immediate.=TRUE); } } else{ diff --git a/R/autoadam.R b/R/autoadam.R index 280fdbb..2110e64 100644 --- a/R/autoadam.R +++ b/R/autoadam.R @@ -503,6 +503,11 @@ auto.adam <- function(data, model="ZXZ", lags=c(frequency(data)), ##### Loop for differences ##### # Prepare table with differences # expand.grid() can be used instead... + # iOrders <- as.matrix(cbind(expand.grid(lapply(iMax, seq, 0)), constant=1)); + # Remove the row with all zeroes, duplicate orders for the constant + # iOrders <- rbind(iOrders[-nrow(iOrders),],iOrders); + # First rows should not have constant + # iOrders[1:(floor(nrow(iOrders)/2)),ncol(iOrders)] <- 0; if(any(iMax!=0)){ iOrders[,1] <- rep(c(0:iMax[1]),times=prod(iMax[-1]+1)); if(ordersLength>1){ diff --git a/R/cma.R b/R/cma.R index 46056f2..cc265b1 100644 --- a/R/cma.R +++ b/R/cma.R @@ -6,7 +6,7 @@ #' shifts it back in time. Otherwise an AR(order+1) model is constructed #' with the preset parameters: #' -#' phi_i = {0.5,1,1,...,0.5} / order +#' \deqn{phi_i = {0.5,1,1,...,0.5} / order} #' #' This then corresponds to the centered MA with 0.5 weight for the #' first observation and 0.5 weight for an additional one. e.g. if this is diff --git a/R/msdecompose.R b/R/msdecompose.R index 9dabb93..757755f 100644 --- a/R/msdecompose.R +++ b/R/msdecompose.R @@ -7,6 +7,7 @@ #' function and order specified in \code{lags} variable in order to smooth the #' original series and obtain level, trend and seasonal components of the series. #' +#' @template smoothRef #' @template ssAuthor #' @template ssKeywords #' @@ -336,6 +337,7 @@ residuals.msdecompose <- function(object, ...){ } } +#' @export sigma.msdecompose <- function(object, ...){ if(errorType(object)=="A"){ return(sqrt(mean(residuals(object)^2,na.rm=TRUE))); diff --git a/R/smooth-package.R b/R/smooth-package.R index fe02abd..42c50c5 100644 --- a/R/smooth-package.R +++ b/R/smooth-package.R @@ -50,7 +50,6 @@ #' } #' #' @name smooth -#' @docType package #' @aliases smooth-package #' @author Ivan Svetunkov #' diff --git a/build/partial.rdb b/build/partial.rdb index 0456224c093b22d0564c9d2d598547fc919122e3..98af4dbb5518edb3e764d662088ae1f89cf6b186 100644 GIT binary patch delta 18 ZcmcDuogl)PHc?c8!(`u!@8%2)3;;6y1zZ3C delta 18 ZcmcDuogl)PGEr24L#VTw!<>PE0RSl{1E&B0 diff --git a/build/vignette.rds b/build/vignette.rds index 31f2b5c2de0d856bae3e86e63f64c4c843c6c3c7..db4de876a17d4737c4a403c812cf3f1197af5aa5 100644 GIT binary patch literal 492 zcmV;MP)?(wDJ?8>qc=av@$Kv zZreBfH`{yQW^dW%F?N5LcR zF4Z{h{Oi%kxz~eXbRIe92B4$J{MEG~>^%D_uXAjeN$*z#Cr(Hwh`ht^JLek?a61TH z;THtP2ghMV*P;w|8r=)qgY4(>+8+->>ysr^po(Cn=?{2*)HUL?(5%*vi%{$s*Wd-W i{eJ)Fts-7B6BLEZdBH4OtZRIa=Z z(Ur7FS|_>XU-wpbRx53!9D2$@vYzH?=FRx+uO%U5NXDZ98IJHcI>&8_>lLmA8RH^k z&Ppa%H?j7>lbaK3nGu1MOQS`zlb{@w zbp4>16r2_%tDQIiz1AgEf%hIIfA~l*=oc;P3hwFUeXSKf6s)2fUr?xBAioaW|8d@7 zo1XL)6c|>gqril7*HH&&5WaZOU32whRTHeb`vHp~ z;THrp1gCLCSE7n`n%oN4gVB%mRX7e}>$5e~pvq{a`FD7@(*|)~Xkql-Di%A#HGIVF ics%}ju83F6Mn&;*Su)2~TZ8Wr{rv}~*BU0o1ONbaH0|{O diff --git a/inst/doc/adam.html b/inst/doc/adam.html index 0bad95d..d4203af 100644 --- a/inst/doc/adam.html +++ b/inst/doc/adam.html @@ -12,7 +12,7 @@ - + Augmented Dynamic Adaptive Model @@ -362,7 +362,7 @@

Augmented Dynamic Adaptive Model

Ivan Svetunkov

-

2023-09-16

+

2024-04-01

@@ -410,35 +410,35 @@

ADAM ETS

#> Model estimated using adam() function: ETS(MMM) #> Response variable: AirPassengers #> Distribution used in the estimation: Normal -#> Loss function type: likelihood; Loss function value: 468.9192 +#> Loss function type: likelihood; Loss function value: 468.8889 #> Coefficients: #> Estimate Std. Error Lower 2.5% Upper 97.5% -#> alpha 0.8451 0.0844 0.6780 1.0000 * -#> beta 0.0205 0.0265 0.0000 0.0727 -#> gamma 0.0000 0.0373 0.0000 0.0736 -#> level 120.2667 14.4276 91.6885 148.7758 * -#> trend 1.0017 0.0100 0.9819 1.0216 * -#> seasonal_1 0.9132 0.0078 0.8987 0.9365 * -#> seasonal_2 0.8996 0.0082 0.8851 0.9229 * -#> seasonal_3 1.0301 0.0096 1.0156 1.0533 * -#> seasonal_4 0.9866 0.0084 0.9721 1.0098 * -#> seasonal_5 0.9852 0.0073 0.9707 1.0085 * -#> seasonal_6 1.1164 0.0095 1.1019 1.1397 * -#> seasonal_7 1.2328 0.0118 1.2183 1.2561 * -#> seasonal_8 1.2262 0.0109 1.2117 1.2495 * -#> seasonal_9 1.0671 0.0093 1.0527 1.0904 * -#> seasonal_10 0.9279 0.0090 0.9134 0.9512 * -#> seasonal_11 0.8058 0.0078 0.7914 0.8291 * +#> alpha 0.8451 0.0841 0.6784 1.0000 * +#> beta 0.0205 0.0264 0.0000 0.0725 +#> gamma 0.0000 0.0381 0.0000 0.0753 +#> level 120.3078 13.1685 94.2236 146.3289 * +#> trend 1.0017 0.0100 0.9818 1.0216 * +#> seasonal_1 0.9150 0.0081 0.9007 0.9381 * +#> seasonal_2 0.9001 0.0081 0.8859 0.9233 * +#> seasonal_3 1.0299 0.0097 1.0156 1.0530 * +#> seasonal_4 0.9867 0.0084 0.9725 1.0099 * +#> seasonal_5 0.9856 0.0072 0.9714 1.0088 * +#> seasonal_6 1.1172 0.0095 1.1030 1.1404 * +#> seasonal_7 1.2344 0.0117 1.2202 1.2576 * +#> seasonal_8 1.2251 0.0106 1.2108 1.2482 * +#> seasonal_9 1.0663 0.0093 1.0521 1.0895 * +#> seasonal_10 0.9261 0.0089 0.9119 0.9493 * +#> seasonal_11 0.8046 0.0077 0.7904 0.8278 * #> -#> Error standard deviation: 0.0376 +#> Error standard deviation: 0.0375 #> Sample size: 132 #> Number of estimated parameters: 17 #> Number of degrees of freedom: 115 #> Information criteria: #> AIC AICc BIC BICc -#> 971.8385 977.2069 1020.8461 1033.9526 +#> 971.7778 977.1463 1020.7855 1033.8919 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 @@ -451,7 +451,7 @@

ADAM ETS

#> Time elapsed: 0.12 seconds #> Model estimated using adam() function: ETS(MMM) #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 468.9192 +#> Loss function type: likelihood; Loss function value: 468.8889 #> Persistence vector g: #> alpha beta gamma #> 0.8451 0.0205 0.0000 @@ -461,17 +461,17 @@

ADAM ETS

#> Number of degrees of freedom: 115 #> Information criteria: #> AIC AICc BIC BICc -#> 971.8385 977.2069 1020.8461 1033.9526 +#> 971.7778 977.1463 1020.7855 1033.8919 #> #> Forecast errors: -#> ME: -4.532; MAE: 15.668; RMSE: 21.652 -#> sCE: -20.72%; Asymmetry: -12.1%; sMAE: 5.969%; sMSE: 0.68% -#> MASE: 0.651; RMSSE: 0.691; rMAE: 0.206; rRMSE: 0.21 +#> ME: -4.589; MAE: 15.524; RMSE: 21.552 +#> sCE: -20.977%; Asymmetry: -12.8%; sMAE: 5.914%; sMSE: 0.674% +#> MASE: 0.645; RMSSE: 0.688; rMAE: 0.204; rRMSE: 0.209

Also, 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):

plot(forecast(testModel,h=18,interval="simulated"))
-

+

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):

@@ -479,7 +479,7 @@

ADAM ETS

plot(testModel,which=c(1:11)) par(mfcol=c(1,1)) plot(testModel,which=12) -

+

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 @@ -495,10 +495,10 @@

ADAM ETS

#> Time elapsed: 0.02 seconds #> Model estimated using adam() function: ETS(AAN) #> Distribution assumed in the model: Normal -#> Loss function type: custom; Loss function value: 599.2829 +#> Loss function type: custom; Loss function value: 599.2251 #> Persistence vector g: #> alpha beta -#> 1.0000 0.2253 +#> 1.0000 0.2266 #> #> Sample size: 138 #> Number of estimated parameters: 4 @@ -506,10 +506,10 @@

ADAM ETS

#> Information criteria are unavailable for the chosen loss & distribution. #> #> Forecast errors: -#> ME: 3.02; MAE: 3.133; RMSE: 3.871 -#> sCE: 15.942%; Asymmetry: 91.8%; sMAE: 1.378%; sMSE: 0.029% -#> MASE: 2.63; RMSSE: 2.523; rMAE: 1.011; rRMSE: 1.01 -

+#> ME: 3.016; MAE: 3.129; RMSE: 3.867 +#> sCE: 15.922%; Asymmetry: 91.7%; sMAE: 1.377%; sMSE: 0.029% +#> MASE: 2.627; RMSSE: 2.521; rMAE: 1.01; 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 @@ -527,7 +527,7 @@

ADAM ETS

nu:

testModel <- adam(BJsales, "MMN", silent=FALSE, distribution="dgnorm", shape=3,
                   h=12, holdout=TRUE)
-

+

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 @@ -540,26 +540,26 @@

ADAM ETS

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.55 seconds +#> Time elapsed: 0.56 seconds #> Model estimated using adam() function: ETS(MAM) #> Distribution assumed in the model: Gamma -#> Loss function type: likelihood; Loss function value: 467.2981 +#> Loss function type: likelihood; Loss function value: 466.5937 #> Persistence vector g: #> alpha beta gamma -#> 0.7691 0.0053 0.0000 +#> 0.7764 0.0000 0.0000 #> #> Sample size: 132 #> Number of estimated parameters: 17 #> Number of degrees of freedom: 115 #> Information criteria: #> AIC AICc BIC BICc -#> 968.5961 973.9646 1017.6038 1030.7102 +#> 967.1875 972.5559 1016.1951 1029.3016 #> #> Forecast errors: -#> ME: 9.537; MAE: 20.784; RMSE: 26.106 -#> sCE: 43.598%; Asymmetry: 64.8%; sMAE: 7.918%; sMSE: 0.989% -#> MASE: 0.863; RMSSE: 0.833; rMAE: 0.273; rRMSE: 0.254 -

+#> ME: 10.678; MAE: 21.661; RMSE: 26.648 +#> sCE: 48.815%; Asymmetry: 66.5%; sMAE: 8.252%; sMSE: 1.031% +#> MASE: 0.899; RMSSE: 0.85; rMAE: 0.285; rRMSE: 0.259 +

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 @@ -572,68 +572,68 @@

ADAM ETS

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 412.6246 389.5669 385.2802 436.2369 -#> Feb 1960 407.7283 378.7151 373.3670 437.6366 -#> Mar 1960 470.0235 429.8654 422.5230 511.6786 -#> Apr 1960 453.1200 411.1176 403.4699 496.8268 -#> May 1960 453.9442 410.7956 402.9500 498.8901 -#> Jun 1960 515.9682 464.4714 455.1338 569.7223 -#> Jul 1960 572.6622 512.9232 502.1196 635.1440 -#> Aug 1960 571.3754 510.0484 498.9773 635.6037 -#> Sep 1960 499.1215 444.9185 435.1408 555.9206 -#> Oct 1960 436.8293 388.4180 379.6965 487.6094 -#> Nov 1960 381.0985 338.1248 330.3919 426.2136 -#> Dec 1960 429.7057 379.5381 370.5320 482.4666 -#> Jan 1961 437.2582 383.1664 373.4963 494.3226 -#> Feb 1961 431.9490 374.9482 364.8088 492.3048 -#> Mar 1961 497.8071 427.2176 414.7363 572.8824 -#> Apr 1961 479.7731 407.9884 395.3571 556.3898 -#> May 1961 480.5163 407.1042 394.2120 558.9822 -#> Jun 1961 546.0243 460.6961 445.7441 637.3716 +#> Jan 1960 412.7963 389.9166 385.6619 436.2219 +#> Feb 1960 407.2683 378.6459 373.3673 436.7623 +#> Mar 1960 468.6856 429.4083 422.2202 509.3979 +#> Apr 1960 451.6047 410.7938 403.3531 494.0283 +#> May 1960 452.4222 410.7801 403.1952 495.7416 +#> Jun 1960 514.3279 465.1691 456.2333 565.5463 +#> Jul 1960 571.2106 514.8504 504.6241 630.0119 +#> Aug 1960 569.9350 512.9405 502.6071 629.4331 +#> Sep 1960 497.4044 447.9023 438.9249 549.0698 +#> Oct 1960 433.9085 390.4578 382.5806 479.2706 +#> Nov 1960 378.7129 340.7720 333.8938 418.3236 +#> Dec 1960 427.5252 383.8846 375.9820 473.1243 +#> Jan 1961 435.1441 388.7155 380.3307 483.7554 +#> Feb 1961 429.2176 381.0829 372.4180 479.7370 +#> Mar 1961 493.8320 434.4863 423.8542 556.3402 +#> Apr 1961 475.7266 416.0207 405.3583 538.7638 +#> May 1961 476.4809 416.3052 405.5640 540.0364 +#> Jun 1961 541.5579 473.1998 460.9976 613.7535 #> Upper bound (97.5%) -#> Jan 1960 440.8930 -#> Feb 1960 443.5806 -#> Mar 1960 520.0179 -#> Apr 1960 505.6094 -#> May 1960 507.9326 -#> Jun 1960 580.5632 -#> Jul 1960 647.7740 -#> Aug 1960 648.6068 -#> Sep 1960 567.4270 -#> Oct 1960 497.9081 -#> Nov 1960 435.3726 -#> Dec 1960 493.1995 -#> Jan 1961 505.9721 -#> Feb 1961 504.6780 -#> Mar 1961 588.3503 -#> Apr 1961 572.2380 -#> May 1961 575.2389 -#> Jun 1961 656.3307 +#> Jan 1960 440.8401 +#> Feb 1960 442.6214 +#> Mar 1960 517.5417 +#> Apr 1960 502.5429 +#> May 1960 504.4435 +#> Jun 1960 575.8535 +#> Jul 1960 641.8638 +#> Aug 1960 641.4335 +#> Sep 1960 559.4879 +#> Oct 1960 488.4205 +#> Nov 1960 426.3137 +#> Dec 1960 482.3312 +#> Jan 1961 493.5936 +#> Feb 1961 489.9897 +#> Mar 1961 569.0779 +#> Apr 1961 551.6443 +#> May 1961 553.0280 +#> Jun 1961 628.5107 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       412.6246          430.9093          436.2369          446.3480
-#> Feb 1960       407.7283          430.8493          437.6366          450.5587
-#> Mar 1960       470.0235          502.1743          511.6786          529.8266
-#> Apr 1960       453.1200          486.8270          496.8268          515.9494
-#> May 1960       453.9442          488.5976          498.8901          518.5818
-#> Jun 1960       515.9682          557.3906          569.7223          593.3382
-#> Jul 1960       572.6622          620.7856          635.1440          662.6663
-#> Aug 1960       571.3754          620.8272          635.6037          663.9448
-#> Sep 1960       499.1215          542.8470          555.9206          581.0019
-#> Oct 1960       436.8293          475.9114          487.6094          510.0618
-#> Nov 1960       381.0985          415.8130          426.2136          446.1838
-#> Dec 1960       429.7057          470.2851          482.4666          505.8752
-#> Jan 1961       437.2582          481.1129          494.3226          519.7428
-#> Feb 1961       431.9490          478.2896          492.3048          519.3197
-#> Mar 1961       497.8071          555.3845          572.8824          606.6772
-#> Apr 1961       479.7731          538.4800          556.3898          591.0344
-#> May 1961       480.5163          540.6183          558.9822          594.5276
-#> Jun 1961       546.0243          615.9648          637.3716          678.8359
+#> Jan 1960 412.7963 430.9374 436.2219 446.2504 +#> Feb 1960 407.2683 430.0713 436.7623 449.4989 +#> Mar 1960 468.6856 500.1145 509.3979 527.1183 +#> Apr 1960 451.6047 484.3306 494.0283 512.5643 +#> May 1960 452.4222 485.8329 495.7416 514.6875 +#> Jun 1960 514.3279 553.8151 565.5463 587.9929 +#> Jul 1960 571.2106 616.5282 630.0119 655.8281 +#> Aug 1960 569.9350 615.7827 629.4331 655.5753 +#> Sep 1960 497.4044 537.2186 549.0698 571.7643 +#> Oct 1960 433.9085 468.8628 479.2706 499.2034 +#> Nov 1960 378.7129 409.2353 418.3236 435.7297 +#> Dec 1960 427.5252 462.6545 473.1243 493.1839 +#> Jan 1961 435.1441 472.5745 483.7554 505.1975 +#> Feb 1961 429.2176 468.0933 479.7370 502.0912 +#> Mar 1961 493.8320 541.8896 556.3402 584.1282 +#> Apr 1961 475.7266 524.1617 538.7638 566.8736 +#> May 1961 476.4809 525.3098 540.0364 568.3902 +#> Jun 1961 541.5579 597.0253 613.7535 645.9606

A 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, @@ -697,26 +697,26 @@

ADAM ETS

testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
                   silent=FALSE, h=336, holdout=TRUE)
 testModel
-#> Time elapsed: 0.79 seconds
+#> Time elapsed: 0.8 seconds
 #> Model estimated using adam() function: ETS(MMdM)[48, 336]
 #> Distribution assumed in the model: Gamma
-#> Loss function type: likelihood; Loss function value: 22013.03
+#> Loss function type: likelihood; Loss function value: 22193.97
 #> Persistence vector g:
 #>  alpha   beta gamma1 gamma2 
-#> 0.9789 0.0000 0.0211 0.0211 
-#> Damping parameter: 0.9988
+#> 0.9825 0.0000 0.0175 0.0175 
+#> Damping parameter: 0.9388
 #> Sample size: 3024
 #> Number of estimated parameters: 6
 #> Number of degrees of freedom: 3018
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 44038.07 44038.10 44074.15 44074.27 
+#> 44399.94 44399.96 44436.02 44436.13 
 #> 
 #> Forecast errors:
-#> ME: 568.186; MAE: 921.821; RMSE: 1133.784
-#> sCE: 628.35%; Asymmetry: 64.2%; sMAE: 3.034%; sMSE: 0.139%
-#> MASE: 1.27; RMSSE: 1.112; rMAE: 0.128; rRMSE: 0.129
-

+#> ME: 41.303; MAE: 723.168; RMSE: 954.172 +#> sCE: 46.478%; Asymmetry: 8.6%; sMAE: 2.422%; sMSE: 0.102% +#> MASE: 0.984; RMSSE: 0.924; rMAE: 0.107; rRMSE: 0.115 +

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 @@ -730,26 +730,26 @@

ADAM ETS

testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
                   silent=FALSE, h=336, holdout=TRUE, maxeval=10000)
 testModel
-#> Time elapsed: 4.92 seconds
+#> Time elapsed: 2.02 seconds
 #> Model estimated using adam() function: ETS(MMdM)[48, 336]
 #> Distribution assumed in the model: Gamma
-#> Loss function type: likelihood; Loss function value: 19493.85
+#> Loss function type: likelihood; Loss function value: 22193.97
 #> Persistence vector g:
 #>  alpha   beta gamma1 gamma2 
-#> 0.0214 0.0000 0.1850 0.2381 
-#> Damping parameter: 0.9008
+#> 0.9825 0.0000 0.0175 0.0175 
+#> Damping parameter: 0.9464
 #> Sample size: 3024
 #> Number of estimated parameters: 6
 #> Number of degrees of freedom: 3018
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 38999.69 38999.72 39035.78 39035.89 
+#> 44399.94 44399.96 44436.02 44436.13 
 #> 
 #> Forecast errors:
-#> ME: 102.323; MAE: 160.013; RMSE: 193.081
-#> sCE: 113.158%; Asymmetry: 66.9%; sMAE: 0.527%; sMSE: 0.004%
-#> MASE: 0.22; RMSSE: 0.189; rMAE: 0.022; rRMSE: 0.022
-

+#> ME: 41.303; MAE: 723.168; RMSE: 954.172 +#> sCE: 46.478%; Asymmetry: 8.6%; sMAE: 2.422%; sMSE: 0.102% +#> MASE: 0.984; RMSSE: 0.924; rMAE: 0.107; rRMSE: 0.115 +

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, @@ -764,26 +764,26 @@

ADAM ETS

testModel <- adam(y, "MMdM", lags=c(1,48,336), initial="backcasting",
                   silent=FALSE, h=336, holdout=TRUE, B=testModel$B)
 testModel
-#> Time elapsed: 0.96 seconds
+#> Time elapsed: 0.67 seconds
 #> Model estimated using adam() function: ETS(MMdM)[48, 336]
 #> Distribution assumed in the model: Gamma
-#> Loss function type: likelihood; Loss function value: 19493.85
+#> Loss function type: likelihood; Loss function value: 22193.97
 #> Persistence vector g:
 #>  alpha   beta gamma1 gamma2 
-#> 0.0214 0.0000 0.1850 0.2381 
-#> Damping parameter: 0.906
+#> 0.9825 0.0000 0.0175 0.0175 
+#> Damping parameter: 0.9537
 #> Sample size: 3024
 #> Number of estimated parameters: 6
 #> Number of degrees of freedom: 3018
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 38999.69 38999.72 39035.78 39035.89 
+#> 44399.94 44399.96 44436.02 44436.13 
 #> 
 #> Forecast errors:
-#> ME: 102.323; MAE: 160.013; RMSE: 193.081
-#> sCE: 113.158%; Asymmetry: 66.9%; sMAE: 0.527%; sMSE: 0.004%
-#> MASE: 0.22; RMSSE: 0.189; rMAE: 0.022; rRMSE: 0.022
-

+#> ME: 41.303; MAE: 723.168; RMSE: 954.172 +#> sCE: 46.478%; Asymmetry: 8.6%; sMAE: 2.422%; sMSE: 0.102% +#> MASE: 0.984; RMSSE: 0.924; rMAE: 0.107; rRMSE: 0.115 +

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 @@ -794,25 +794,25 @@

ADAM ETS

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: 0.62 seconds
+#> Time elapsed: 0.57 seconds
 #> Model estimated using adam() function: ETS(MMdM)[48, 336]
 #> Distribution assumed in the model: Gamma
-#> Loss function type: likelihood; Loss function value: 21830.26
+#> Loss function type: likelihood; Loss function value: 21765.29
 #> Persistence vector g:
 #>  alpha   beta gamma1 gamma2 
-#> 0.9671 0.1000 0.0318 0.0329 
-#> Damping parameter: 0.6265
+#> 0.9344 0.1000 0.0638 0.0656 
+#> Damping parameter: 0.4112
 #> Sample size: 3024
 #> Number of estimated parameters: 5
 #> Number of degrees of freedom: 3019
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 43670.52 43670.54 43700.60 43700.68 
+#> 43540.58 43540.60 43570.65 43570.73 
 #> 
 #> Forecast errors:
-#> ME: 612.352; MAE: 945.443; RMSE: 1160.15
-#> sCE: 677.192%; Asymmetry: 66.9%; sMAE: 3.112%; sMSE: 0.146%
-#> MASE: 1.303; RMSSE: 1.138; rMAE: 0.131; rRMSE: 0.132
+#> ME: -18.126; MAE: 692.908; RMSE: 915.669 +#> sCE: -20.397%; Asymmetry: 0.1%; sMAE: 2.321%; sMSE: 0.094% +#> MASE: 0.943; RMSSE: 0.887; rMAE: 0.103; rRMSE: 0.11

The 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 @@ -820,25 +820,25 @@

ADAM ETS

testModel <- adam(rpois(120,0.5), "MNN", silent=FALSE, h=12, holdout=TRUE,
                   occurrence="odds-ratio")
 testModel
-#> Time elapsed: 0.05 seconds
+#> Time elapsed: 0.04 seconds
 #> Model estimated using adam() function: iETS(MNN)[O]
 #> Occurrence model type: Odds ratio
 #> Distribution assumed in the model: Mixture of Bernoulli and Gamma
-#> Loss function type: likelihood; Loss function value: 68.7142
+#> Loss function type: likelihood; Loss function value: 39.789
 #> Persistence vector g:
-#> alpha 
-#>     0 
+#>  alpha 
+#> 0.7545 
 #> 
 #> Sample size: 108
 #> Number of estimated parameters: 5
 #> Number of degrees of freedom: 103
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 294.7652 294.9960 308.1758 299.3518 
+#> 232.9765 233.2073 246.3872 237.5632 
 #> 
 #> Forecast errors:
-#> Asymmetry: -18.254%; sMSE: 23.078%; rRMSE: 0.793; sPIS: 368.43%; sCE: -45.267%
-

+#> Asymmetry: -31.604%; sMSE: 38.824%; rRMSE: 0.794; sPIS: -79.727%; sCE: 31.189% +

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. @@ -850,35 +850,35 @@

ADAM ETS

"adam:" #> [1] "adam:" adamModel -#> Time elapsed: 2.17 seconds +#> Time elapsed: 2.25 seconds #> Model estimated: ETS(CCC) #> Loss function type: likelihood #> #> Number of models combined: 30 #> Sample size: 132 -#> Average number of estimated parameters: 21.8083 -#> Average number of degrees of freedom: 110.1917 +#> Average number of estimated parameters: 20.6558 +#> Average number of degrees of freedom: 111.3442 #> #> Forecast errors: -#> ME: -1.759; MAE: 14.402; RMSE: 20.569 -#> sCE: -8.04%; sMAE: 5.487%; sMSE: 0.614% -#> MASE: 0.598; RMSSE: 0.656; rMAE: 0.19; rRMSE: 0.2 +#> ME: 2.537; MAE: 15.688; RMSE: 21.994 +#> sCE: 11.599%; sMAE: 5.976%; sMSE: 0.702% +#> MASE: 0.651; RMSSE: 0.702; rMAE: 0.206; rRMSE: 0.214 "es():" #> [1] "es():" esModel -#> Time elapsed: 2.07 seconds +#> Time elapsed: 2.14 seconds #> Model estimated: ETS(CCC) #> Loss function type: likelihood #> #> Number of models combined: 30 #> Sample size: 132 -#> Average number of estimated parameters: 21.3232 -#> Average number of degrees of freedom: 110.6768 +#> Average number of estimated parameters: 21.1644 +#> Average number of degrees of freedom: 110.8356 #> #> Forecast errors: -#> ME: 2.906; MAE: 15.858; RMSE: 22.111 -#> sCE: 13.287%; sMAE: 6.041%; sMSE: 0.71% -#> MASE: 0.658; RMSSE: 0.706; rMAE: 0.209; rRMSE: 0.215 +#> ME: 4.754; MAE: 17.114; RMSE: 22.996 +#> sCE: 21.733%; sMAE: 6.52%; sMSE: 0.767% +#> MASE: 0.711; RMSSE: 0.734; rMAE: 0.225; rRMSE: 0.223

ADAM ARIMA

@@ -889,7 +889,7 @@

ADAM ARIMA

testModel <- adam(BJsales, "NNN", silent=FALSE, orders=c(0,2,2),
                   h=12, holdout=TRUE)
 testModel
-#> Time elapsed: 0.04 seconds
+#> Time elapsed: 0.03 seconds
 #> Model estimated using adam() function: ARIMA(0,2,2)
 #> Distribution assumed in the model: Normal
 #> Loss function type: likelihood; Loss function value: 240.5944
@@ -923,60 +923,60 @@ 
ADAM ARIMA

                   orders=list(ar=c(1,1),i=c(1,1),ma=c(2,2)), distribution="dlnorm",
                   h=12, holdout=TRUE)
 testModel
-#> Time elapsed: 0.58 seconds
+#> Time elapsed: 0.37 seconds
 #> Model estimated using adam() function: SARIMA(1,1,2)[1](1,1,2)[12]
 #> Distribution assumed in the model: Log-Normal
-#> Loss function type: likelihood; Loss function value: 506.1549
+#> Loss function type: likelihood; Loss function value: 617.5117
 #> ARMA parameters of the model:
 #> AR:
 #>  phi1[1] phi1[12] 
-#>   0.0925   0.1077 
+#>   0.8286   0.4681 
 #> MA:
 #>  theta1[1]  theta2[1] theta1[12] theta2[12] 
-#>     0.0703     0.3247    -0.4464    -0.0607 
+#>    -0.5840     0.0896    -0.1938    -0.1226 
 #> 
 #> Sample size: 132
 #> Number of estimated parameters: 33
 #> Number of degrees of freedom: 99
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 1078.310 1101.208 1173.442 1229.345 
+#> 1301.024 1323.921 1396.156 1452.059 
 #> 
 #> Forecast errors:
-#> ME: -33.764; MAE: 33.764; RMSE: 37.234
-#> sCE: -154.352%; Asymmetry: -100%; sMAE: 12.863%; sMSE: 2.012%
-#> MASE: 1.402; RMSSE: 1.188; rMAE: 0.444; rRMSE: 0.362
-

+#> ME: -94.42; MAE: 94.42; RMSE: 100.143 +#> sCE: -431.648%; Asymmetry: -100%; sMAE: 35.971%; sMSE: 14.555% +#> MASE: 3.92; RMSSE: 3.196; rMAE: 1.242; rRMSE: 0.972
+

Second, if you want the model with intercept / drift, you can do it using constant parameter:

testModel <- adam(AirPassengers, "NNN", silent=FALSE, lags=c(1,12), constant=TRUE,
                   orders=list(ar=c(1,1),i=c(1,1),ma=c(2,2)), distribution="dnorm",
                   h=12, holdout=TRUE)
 testModel
-#> Time elapsed: 0.49 seconds
+#> Time elapsed: 0.34 seconds
 #> Model estimated using adam() function: SARIMA(1,1,2)[1](1,1,2)[12] with drift
 #> Distribution assumed in the model: Normal
-#> Loss function type: likelihood; Loss function value: 489.421
+#> Loss function type: likelihood; Loss function value: 493.1328
 #> ARMA parameters of the model:
 #> AR:
 #>  phi1[1] phi1[12] 
-#>  -0.6043  -0.2731 
+#>  -0.8840  -0.2784 
 #> MA:
 #>  theta1[1]  theta2[1] theta1[12] theta2[12] 
-#>     0.3958     0.0039     0.1525     0.1035 
+#>     0.8752    -0.1018    -0.1000     0.0965 
 #> 
 #> Sample size: 132
 #> Number of estimated parameters: 34
 #> Number of degrees of freedom: 98
 #> Information criteria:
 #>      AIC     AICc      BIC     BICc 
-#> 1046.842 1071.378 1144.857 1204.760 
+#> 1054.265 1078.802 1152.281 1212.183 
 #> 
 #> Forecast errors:
-#> ME: -16.487; MAE: 18.657; RMSE: 23.574
-#> sCE: -75.371%; Asymmetry: -85.2%; sMAE: 7.108%; sMSE: 0.807%
-#> MASE: 0.775; RMSSE: 0.752; rMAE: 0.245; rRMSE: 0.229
-

+#> ME: -6.83; MAE: 13.958; RMSE: 17.923 +#> sCE: -31.226%; Asymmetry: -49.2%; sMAE: 5.317%; sMSE: 0.466% +#> MASE: 0.58; RMSSE: 0.572; rMAE: 0.184; rRMSE: 0.174 +

If the model contains non-zero differences, then the constant acts as a drift. Third, you can specify parameters of ARIMA via the arma parameter in the following manner:

@@ -985,10 +985,10 @@

ADAM ARIMA

arma=list(ar=c(0.1,0.1), ma=c(-0.96, 0.03, -0.12, 0.03)), h=12, holdout=TRUE) testModel -#> Time elapsed: 0.26 seconds +#> Time elapsed: 0.18 seconds #> Model estimated using adam() function: SARIMA(1,1,2)[1](1,1,2)[12] #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 534.9627 +#> Loss function type: likelihood; Loss function value: 534.8804 #> ARMA parameters of the model: #> AR: #> phi1[1] phi1[12] @@ -1002,13 +1002,13 @@

ADAM ARIMA

#> Number of degrees of freedom: 105 #> Information criteria: #> AIC AICc BIC BICc -#> 1123.925 1138.464 1201.761 1237.255 +#> 1123.761 1138.299 1201.596 1237.091 #> #> Forecast errors: #> ME: 9.575; MAE: 17.082; RMSE: 19.148 #> sCE: 43.773%; Asymmetry: 61.9%; sMAE: 6.508%; sMSE: 0.532% #> MASE: 0.709; RMSSE: 0.611; rMAE: 0.225; rRMSE: 0.186 -

+

Finally, the initials for the states can also be provided, although getting the correct ones might be a challenging task (you also need to know how many of them to provide; checking @@ -1018,30 +1018,30 @@

ADAM ARIMA

initial=list(arima=AirPassengers[1:24]), h=12, holdout=TRUE) testModel -#> Time elapsed: 0.31 seconds +#> Time elapsed: 0.28 seconds #> Model estimated using adam() function: SARIMA(1,1,2)[1](1,1,0)[12] #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 493.7137 +#> Loss function type: likelihood; Loss function value: 492.2311 #> ARMA parameters of the model: #> AR: #> phi1[1] phi1[12] -#> 0.100 0.153 +#> 0.1018 0.0950 #> MA: #> theta1[1] theta2[1] -#> -0.3336 0.0581 +#> -0.3302 0.0588 #> #> Sample size: 132 #> Number of estimated parameters: 31 #> Number of degrees of freedom: 101 #> Information criteria: #> AIC AICc BIC BICc -#> 1049.427 1069.267 1138.794 1187.232 +#> 1046.462 1066.302 1135.829 1184.266 #> #> Forecast errors: -#> ME: -23.835; MAE: 23.967; RMSE: 29.528 -#> sCE: -108.963%; Asymmetry: -97.9%; sMAE: 9.13%; sMSE: 1.265% -#> MASE: 0.995; RMSSE: 0.942; rMAE: 0.315; rRMSE: 0.287 -

+#> ME: -21.3; MAE: 21.948; RMSE: 27.37 +#> sCE: -97.376%; Asymmetry: -95%; sMAE: 8.361%; sMSE: 1.087% +#> MASE: 0.911; RMSSE: 0.874; rMAE: 0.289; rRMSE: 0.266 +

If you work with ADAM ARIMA model, then there is no such thing as “usual” bounds for the parameters, so the function will use the bounds="admissible", checking the AR / MA polynomials in @@ -1068,7 +1068,7 @@

ADAM ETSX / ARIMAX / ETSX+ARIMA

brief example:

BJData <- cbind(BJsales,BJsales.lead)
 testModel <- adam(BJData, "AAN", h=18, silent=FALSE)
-

+

If you work with data.frame or similar structures, then you can use them directly, ADAM will extract the response variable either assuming that it is in the first column or from the provided formula (if you @@ -1082,7 +1082,7 @@

ADAM ETSX / ARIMAX / ETSX+ARIMA

#> Time elapsed: 0.13 seconds #> Model estimated using adam() function: ETSX(ANN) #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 210.9197 +#> Loss function type: likelihood; Loss function value: 211.0843 #> Persistence vector g (excluding xreg): #> alpha #> 1 @@ -1092,13 +1092,13 @@

ADAM ETSX / ARIMAX / ETSX+ARIMA

#> Number of degrees of freedom: 126 #> Information criteria: #> AIC AICc BIC BICc -#> 433.8393 434.5113 451.1361 452.7768 +#> 434.1687 434.8407 451.4655 453.1061 #> #> Forecast errors: -#> ME: 0.744; MAE: 1.299; RMSE: 1.782 -#> sCE: 5.924%; Asymmetry: 44.2%; sMAE: 0.575%; sMSE: 0.006% -#> MASE: 1.065; RMSSE: 1.141; rMAE: 0.58; rRMSE: 0.71 -

+#> ME: 0.758; MAE: 1.304; RMSE: 1.79 +#> sCE: 6.042%; Asymmetry: 45.1%; sMAE: 0.577%; sMSE: 0.006% +#> MASE: 1.069; RMSSE: 1.145; rMAE: 0.582; rRMSE: 0.713 +

Similarly to es(), there is a support for variables selection, but via the regressors parameter instead of xregDo, which will then use stepwise() @@ -1189,29 +1189,29 @@

ADAM ETSX / ARIMAX / ETSX+ARIMA

#> Model estimated using adam() function: ETSX(AAN)+ARIMA(1,0,0) #> Response variable: y #> Distribution used in the estimation: Normal -#> Loss function type: likelihood; Loss function value: 74.487 +#> Loss function type: likelihood; Loss function value: 73.9255 #> Coefficients: #> Estimate Std. Error Lower 2.5% Upper 97.5% -#> alpha 0.5579 0.1044 0.3512 0.7642 * -#> beta 0.0000 0.0147 0.0000 0.0291 -#> phi1[1] 0.7581 0.1384 0.4842 1.0000 * -#> level 42.1940 4.6637 32.9602 51.4122 * -#> trend 0.0752 0.0253 0.0251 0.1252 * -#> ARIMAState1 1.9223 1.3764 -0.8029 4.6429 -#> xLag3 5.0456 0.1245 4.7991 5.2917 * -#> xLag7 0.9462 0.1325 0.6838 1.2082 * -#> xLag4 4.1336 0.1557 3.8254 4.4413 * -#> xLag6 2.3273 0.1568 2.0169 2.6371 * -#> xLag5 3.1050 0.1598 2.7886 3.4210 * +#> alpha 0.3774 0.1257 0.1286 0.6258 * +#> beta 0.0000 0.0128 0.0000 0.0253 +#> phi1[1] 0.9062 0.0847 0.7385 1.0000 * +#> level 40.0446 9.1078 22.0118 58.0469 * +#> trend 0.0554 0.0275 0.0010 0.1097 * +#> ARIMAState1 0.5975 2.3814 -4.1175 5.3045 +#> xLag3 5.1053 0.1384 4.8313 5.3788 * +#> xLag7 0.9556 0.1494 0.6598 1.2508 * +#> xLag4 4.2029 0.2121 3.7829 4.6221 * +#> xLag6 2.3560 0.2469 1.8671 2.8440 * +#> xLag5 3.2282 0.2567 2.7198 3.7356 * #> -#> Error standard deviation: 0.4462 +#> Error standard deviation: 0.4443 #> Sample size: 132 #> Number of estimated parameters: 12 #> Number of degrees of freedom: 120 #> Information criteria: #> AIC AICc BIC BICc -#> 172.9740 175.5959 207.5677 213.9686 -

+#> 171.8510 174.4728 206.4446 212.8456 +

This might be handy, when you explore a high frequency data, want to add calendar events, apply ETS and add AR/MA errors to it.

Finally, if you estimate ETSX or ARIMAX model and want to speed @@ -1224,24 +1224,24 @@

ADAM ETSX / ARIMAX / ETSX+ARIMA

#> Model estimated using adam() function: ETSX(AAN) #> Response variable: y #> Distribution used in the estimation: Normal -#> Loss function type: likelihood; Loss function value: 62.3464 +#> Loss function type: likelihood; Loss function value: 62.3043 #> Coefficients: #> Estimate Std. Error Lower 2.5% Upper 97.5% #> alpha 1.0000 0.0871 0.8277 1.0000 * #> beta 0.0047 0.0101 0.0000 0.0246 -#> xLag3 4.8015 3.3300 -1.7895 11.3852 -#> xLag7 0.7732 3.3408 -5.8392 7.3781 -#> xLag4 3.6139 3.0342 -2.3917 9.6128 -#> xLag6 1.6868 3.0303 -4.3111 7.6780 -#> xLag5 2.4306 2.8286 -3.1679 8.0229 +#> xLag3 4.8014 3.3115 -1.7530 11.3484 +#> xLag7 0.7730 3.3222 -5.8027 7.3412 +#> xLag4 3.6140 3.0176 -2.3587 9.5800 +#> xLag6 1.6868 3.0136 -4.2780 7.6448 +#> xLag5 2.4304 2.8129 -3.1371 7.9918 #> -#> Error standard deviation: 0.4004 +#> Error standard deviation: 0.4002 #> Sample size: 132 #> Number of estimated parameters: 8 #> Number of degrees of freedom: 124 #> Information criteria: #> AIC AICc BIC BICc -#> 140.6929 141.8636 163.7553 166.6135 +#> 140.6086 141.7793 163.6710 166.5292

Auto ADAM

@@ -1262,29 +1262,29 @@

Auto ADAM

#> Selecting ARMA... | #> The best ARIMA is selected. ds , Selecting ARIMA orders... #> Selecting differences... -#> Selecting ARMA... | +#> Selecting ARMA... |- #> The best ARIMA is selected. Done! testModel -#> Time elapsed: 0.44 seconds -#> Model estimated using auto.adam() function: ETS(AAdN) with drift +#> Time elapsed: 0.49 seconds +#> Model estimated using auto.adam() function: ETS(AAdN) #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 237.1294 +#> Loss function type: likelihood; Loss function value: 238.0241 #> Persistence vector g: -#> alpha beta -#> 0.9541 0.2851 -#> Damping parameter: 0.8363 +#> alpha beta +#> 0.955 0.296 +#> Damping parameter: 0.8456 #> Sample size: 138 -#> Number of estimated parameters: 7 -#> Number of degrees of freedom: 131 +#> Number of estimated parameters: 6 +#> Number of degrees of freedom: 132 #> Information criteria: #> AIC AICc BIC BICc -#> 488.2588 489.1204 508.7496 510.8721 +#> 488.0481 488.6893 505.6116 507.1914 #> #> Forecast errors: -#> ME: 0.622; MAE: 1.2; RMSE: 1.5 -#> sCE: 3.286%; Asymmetry: 48.9%; sMAE: 0.528%; sMSE: 0.004% -#> MASE: 1.007; RMSSE: 0.978; rMAE: 0.387; rRMSE: 0.391
-

+#> ME: 2.807; MAE: 2.966; RMSE: 3.65 +#> sCE: 14.815%; Asymmetry: 87.5%; sMAE: 1.305%; sMSE: 0.026% +#> MASE: 2.489; RMSSE: 2.379; rMAE: 0.957; rRMSE: 0.952 +

This process can also be done in parallel on either the automatically selected number of cores (e.g. parallel=TRUE) or on the specified by user (e.g. parallel=4):

@@ -1297,33 +1297,33 @@

Auto ADAM

distribution=c("dnorm","dlaplace","ds","dgnorm"), h=12, holdout=TRUE) testModel -#> Time elapsed: 0.47 seconds +#> Time elapsed: 0.38 seconds #> Model estimated using auto.adam() function: ETS(AAN)+ARIMA(2,2,2) #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 243.9582 +#> Loss function type: likelihood; Loss function value: 242.2371 #> Persistence vector g: #> alpha beta -#> 0.0143 0.0000 +#> 0.0194 0.0000 #> #> ARMA parameters of the model: #> AR: #> phi1[1] phi2[1] -#> -0.7434 -0.2804 +#> -0.5962 -0.1945 #> MA: #> theta1[1] theta2[1] -#> -0.0786 -0.2418 +#> -0.2037 -0.2357 #> #> Sample size: 138 #> Number of estimated parameters: 13 #> Number of degrees of freedom: 125 #> Information criteria: #> AIC AICc BIC BICc -#> 513.9164 516.8519 551.9707 559.2026 +#> 510.4741 513.4096 548.5284 555.7603 #> #> Forecast errors: -#> ME: 2.745; MAE: 2.901; RMSE: 3.554 -#> sCE: 14.492%; Asymmetry: 87.9%; sMAE: 1.276%; sMSE: 0.024% -#> MASE: 2.435; RMSSE: 2.316; rMAE: 0.936; rRMSE: 0.927 +#> ME: 2.803; MAE: 2.95; RMSE: 3.622 +#> sCE: 14.795%; Asymmetry: 88.4%; sMAE: 1.297%; sMSE: 0.025% +#> MASE: 2.476; RMSSE: 2.361; rMAE: 0.951; rRMSE: 0.945

However, this way the function will just use ARIMA(2,2,2) and fit it together with ETS. If you want it to select the most appropriate ARIMA orders from the provided (e.g. up to AR(2), I(1) and MA(2)), you need to @@ -1337,26 +1337,26 @@

Auto ADAM

#> Selecting ARMA... | #> The best ARIMA is selected. Done! testModel -#> Time elapsed: 0.17 seconds -#> Model estimated using auto.adam() function: ETS(AAdN) with drift +#> Time elapsed: 0.16 seconds +#> Model estimated using auto.adam() function: ETS(AAdN) #> Distribution assumed in the model: Normal -#> Loss function type: likelihood; Loss function value: 237.1294 +#> Loss function type: likelihood; Loss function value: 238.0241 #> Persistence vector g: -#> alpha beta -#> 0.9541 0.2851 -#> Damping parameter: 0.8363 +#> alpha beta +#> 0.955 0.296 +#> Damping parameter: 0.8456 #> Sample size: 138 -#> Number of estimated parameters: 7 -#> Number of degrees of freedom: 131 +#> Number of estimated parameters: 6 +#> Number of degrees of freedom: 132 #> Information criteria: #> AIC AICc BIC BICc -#> 488.2588 489.1204 508.7496 510.8721 +#> 488.0481 488.6893 505.6116 507.1914 #> #> Forecast errors: -#> ME: 0.622; MAE: 1.2; RMSE: 1.5 -#> sCE: 3.286%; Asymmetry: 48.9%; sMAE: 0.528%; sMSE: 0.004% -#> MASE: 1.007; RMSSE: 0.978; rMAE: 0.387; rRMSE: 0.391 -

+#> ME: 2.807; MAE: 2.966; RMSE: 3.65 +#> sCE: 14.815%; Asymmetry: 87.5%; sMAE: 1.305%; sMSE: 0.026% +#> MASE: 2.489; RMSSE: 2.379; rMAE: 0.957; rRMSE: 0.952 +

Knowing how to work with adam(), you can use similar principles, when dealing with auto.adam(). Just keep in mind that the provided persistence, phi, @@ -1376,26 +1376,26 @@

Auto ADAM

#> The best ARIMA is selected. #> Dealing with outliers... testModel -#> Time elapsed: 4.94 seconds +#> Time elapsed: 4.71 seconds #> Model estimated using auto.adam() function: ETSX(MMM) #> Distribution assumed in the model: Gamma -#> Loss function type: likelihood; Loss function value: 464.5879 +#> Loss function type: likelihood; Loss function value: 464.5904 #> Persistence vector g (excluding xreg): #> alpha beta gamma -#> 0.7633 0.0000 0.0000 +#> 0.7634 0.0000 0.0000 #> #> Sample size: 132 #> Number of estimated parameters: 18 #> Number of degrees of freedom: 114 #> Information criteria: #> AIC AICc BIC BICc -#> 965.1758 971.2289 1017.0663 1031.8443 +#> 965.1809 971.2340 1017.0713 1031.8494 #> #> Forecast errors: -#> ME: -8.461; MAE: 15.346; RMSE: 22.431 -#> sCE: -38.678%; Asymmetry: -49%; sMAE: 5.846%; sMSE: 0.73% -#> MASE: 0.637; RMSSE: 0.716; rMAE: 0.202; rRMSE: 0.218 -

+#> ME: -8.629; MAE: 15.42; RMSE: 22.535 +#> sCE: -39.448%; Asymmetry: -50%; sMAE: 5.875%; sMSE: 0.737% +#> MASE: 0.64; RMSSE: 0.719; rMAE: 0.203; rRMSE: 0.219 +

If you specify outliers="select", the function will create leads and lags 1 of the outliers and then select the most appropriate ones via the regressors parameter of adam.

diff --git a/inst/doc/ces.html b/inst/doc/ces.html index 9bab67a..c24b859 100644 --- a/inst/doc/ces.html +++ b/inst/doc/ces.html @@ -12,7 +12,7 @@ - + ces() - Complex Exponential Smoothing @@ -340,7 +340,7 @@

ces() - Complex Exponential Smoothing

Ivan Svetunkov

-

2023-09-16

+

2024-04-01

@@ -356,7 +356,7 @@

2023-09-16

For the same series from M3 dataset ces() can be constructed using:

ces(BJsales, h=12, holdout=TRUE, silent=FALSE)
-
## Time elapsed: 0.12 seconds
+
## Time elapsed: 0.13 seconds
 ## Model estimated: CES(n)
 ## a0 + ia1: 1.9981+1.0034i
 ## Initial values were produced using backcasting.
@@ -381,7 +381,7 @@ 
2023-09-16

 
If we want automatic model selection, then we use
 auto.ces() function:

 
auto.ces(BJsales, h=12, holdout=TRUE, interval="p", silent=FALSE)

-
## Time elapsed: 0.04 seconds
+
## Time elapsed: 0.05 seconds
 ## Model estimated: CES(n)
 ## a0 + ia1: 1.9981+1.0034i
 ## Initial values were produced using backcasting.
diff --git a/inst/doc/es.html b/inst/doc/es.html
index 9732140..0f814c1 100644
--- a/inst/doc/es.html
+++ b/inst/doc/es.html
@@ -12,7 +12,7 @@
 
 
 
-
+
 
 es() - Exponential Smoothing
 
@@ -340,7 +340,7 @@
 
 
es() - Exponential Smoothing

 
Ivan Svetunkov

-
2023-09-16

+
2024-04-01

 
 
 
@@ -358,26 +358,26 @@ 
2023-09-16

 
ourModel <- es(BJsales, h=12, holdout=TRUE, silent=FALSE)

 
## Forming the pool of models based on... ANN , AAN , Estimation progress:    33 %44 %56 %67 %78 %89 %100 %... Done!

 
ourModel

-
## Time elapsed: 0.2 seconds
+
## Time elapsed: 0.21 seconds
 ## Model estimated using es() function: ETS(AMdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 237.7324
+## Loss function type: likelihood; Loss function value: 237.8564
 ## Persistence vector g:
-## alpha  beta 
-## 1.000 0.249 
-## Damping parameter: 0.9034
+##  alpha   beta 
+## 0.9984 0.2184 
+## Damping parameter: 0.9012
 ## Sample size: 138
 ## Number of estimated parameters: 6
 ## Number of degrees of freedom: 132
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 487.4647 488.1060 505.0283 506.6080 
+## 487.7127 488.3540 505.2763 506.8560 
 ## 
 ## Forecast errors:
-## ME: 2.83; MAE: 2.986; RMSE: 3.676
-## sCE: 14.939%; Asymmetry: 87.8%; sMAE: 1.314%; sMSE: 0.026%
-## MASE: 2.507; RMSSE: 2.396; rMAE: 0.963; rRMSE: 0.959

-

+## ME: 2.888; MAE: 3.029; RMSE: 3.734
+## sCE: 15.244%; Asymmetry: 88.8%; sMAE: 1.332%; sMSE: 0.027%
+## MASE: 2.543; RMSSE: 2.434; rMAE: 0.977; rRMSE: 0.975

+

 
In this case function uses branch and bound algorithm to form a pool
 of models to check and after that constructs a model with the lowest
 information criterion. As we can see, it also produces an output with
@@ -399,24 +399,24 @@ 
2023-09-16

 
If we need prediction interval, then we can use the
 forecast() method:

 
plot(forecast(ourModel, h=12, interval="prediction"))

-

+

 
The same model can be reused for different purposes, for example to
 produce forecasts based on newly available data:

 
es(BJsales, model=ourModel, h=12, holdout=FALSE)

 
## Time elapsed: 0 seconds
 ## Model estimated using es() function: ETS(AMdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 255.4039
+## Loss function type: likelihood; Loss function value: 255.555
 ## Persistence vector g:
-## alpha  beta 
-## 1.000 0.249 
-## Damping parameter: 0.9034
+##  alpha   beta 
+## 0.9984 0.2184 
+## Damping parameter: 0.9012
 ## Sample size: 150
 ## Number of estimated parameters: 1
 ## Number of degrees of freedom: 149
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 512.8078 512.8348 515.8184 515.8862

+## 513.1099 513.1370 516.1206 516.1883

 
We can also extract the type of model in order to reuse it later:

 
modelType(ourModel)

 
## [1] "AMdN"

@@ -450,17 +450,17 @@ 
2023-09-16

 
## Time elapsed: 0.02 seconds
 ## Model estimated using es() function: ETS(AMdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 255.2728
+## Loss function type: likelihood; Loss function value: 255.2509
 ## Persistence vector g:
 ##  alpha   beta 
-## 0.9716 0.2812 
-## Damping parameter: 0.8747
+## 0.9719 0.2835 
+## Damping parameter: 0.8698
 ## Sample size: 150
 ## Number of estimated parameters: 4
 ## Number of degrees of freedom: 146
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 518.5456 518.8215 530.5881 531.2793

+## 518.5019 518.7777 530.5444 531.2355

 
# Provided persistence
 es(BJsales, model=modelType(ourModel),
    h=12, holdout=FALSE,
@@ -468,17 +468,17 @@ 
2023-09-16

 
## Time elapsed: 0.02 seconds
 ## Model estimated using es() function: ETS(AMdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 255.3497
+## Loss function type: likelihood; Loss function value: 256.1161
 ## Persistence vector g:
-## alpha  beta 
-## 1.000 0.249 
-## Damping parameter: 0.888
+##  alpha   beta 
+## 0.9984 0.2184 
+## Damping parameter: 0.944
 ## Sample size: 150
 ## Number of estimated parameters: 4
 ## Number of degrees of freedom: 146
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 518.6993 518.9752 530.7419 531.4330

+## 520.2321 520.5080 532.2747 532.9658

 
or provide some arbitrary values:

 
es(BJsales, model=modelType(ourModel),
    h=12, holdout=FALSE,
@@ -486,81 +486,81 @@ 
2023-09-16

 
## Time elapsed: 0.03 seconds
 ## Model estimated using es() function: ETS(AMdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 255.3529
+## Loss function type: likelihood; Loss function value: 255.3409
 ## Persistence vector g:
 ##  alpha   beta 
-## 0.9893 0.2677 
-## Damping parameter: 0.8928
+## 0.9882 0.2671 
+## Damping parameter: 0.8925
 ## Sample size: 150
 ## Number of estimated parameters: 5
 ## Number of degrees of freedom: 145
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 520.7058 521.1224 535.7589 536.8028

+## 520.6818 521.0984 535.7349 536.7788

 
Using some other parameters may lead to completely different model
 and forecasts (see discussion of the additional parameters in the online textbook about
 ADAM):

 
es(BJsales, h=12, holdout=TRUE, loss="MSEh", bounds="a", ic="BIC")

-
## Time elapsed: 0.79 seconds
+
## Time elapsed: 0.84 seconds
 ## Model estimated using es() function: ETS(MAN)
 ## Distribution assumed in the model: Normal
 ## Loss function type: MSEh; Loss function value: 0.0018
 ## Persistence vector g:
 ##  alpha   beta 
-## 1.5008 0.0000 
+## 1.5365 0.0002 
 ## 
 ## Sample size: 138
 ## Number of estimated parameters: 4
 ## Number of degrees of freedom: 134
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 1022.671 1022.971 1034.380 1035.121 
+## 1021.898 1022.199 1033.607 1034.348 
 ## 
 ## Forecast errors:
-## ME: -0.446; MAE: 1.206; RMSE: 1.345
-## sCE: -2.356%; Asymmetry: -43.6%; sMAE: 0.531%; sMSE: 0.004%
-## MASE: 1.012; RMSSE: 0.877; rMAE: 0.389; rRMSE: 0.351

+## ME: -1.396; MAE: 1.604; RMSE: 1.809
+## sCE: -7.369%; Asymmetry: -84.7%; sMAE: 0.705%; sMSE: 0.006%
+## MASE: 1.346; RMSSE: 1.179; rMAE: 0.517; rRMSE: 0.472

 
You can play around with all the available parameters to see what’s
 their effect on the final model.

 
In order to combine forecasts we need to use “C” letter:

 
es(BJsales, model="CCN", h=12, holdout=TRUE)

-
## Time elapsed: 0.22 seconds
+
## Time elapsed: 0.23 seconds
 ## Model estimated: ETS(CCN)
 ## Loss function type: likelihood
 ## 
 ## Number of models combined: 10
 ## Sample size: 138
-## Average number of estimated parameters: 6.3367
-## Average number of degrees of freedom: 131.6633
+## Average number of estimated parameters: 6.2884
+## Average number of degrees of freedom: 131.7116
 ## 
 ## Forecast errors:
-## ME: 2.842; MAE: 2.993; RMSE: 3.686
-## sCE: 15.002%; sMAE: 1.317%; sMSE: 0.026%
-## MASE: 2.513; RMSSE: 2.403; rMAE: 0.966; rRMSE: 0.962

+## ME: 2.861; MAE: 3.007; RMSE: 3.706
+## sCE: 15.102%; sMAE: 1.323%; sMSE: 0.027%
+## MASE: 2.524; RMSSE: 2.416; rMAE: 0.97; rRMSE: 0.967

 
Model selection from a specified pool and forecasts combination are
 called using respectively:

 
# Select the best model in the pool
 es(BJsales, model=c("ANN","AAN","AAdN","MNN","MAN","MAdN"),
    h=12, holdout=TRUE)

-
## Time elapsed: 0.1 seconds
+
## Time elapsed: 0.11 seconds
 ## Model estimated using es() function: ETS(AAdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 238.2715
+## Loss function type: likelihood; Loss function value: 238.0241
 ## Persistence vector g:
-##  alpha   beta 
-## 0.9534 0.2925 
-## Damping parameter: 0.8622
+## alpha  beta 
+## 0.955 0.296 
+## Damping parameter: 0.8456
 ## Sample size: 138
 ## Number of estimated parameters: 6
 ## Number of degrees of freedom: 132
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 488.5431 489.1843 506.1066 507.6863 
+## 488.0481 488.6893 505.6116 507.1914 
 ## 
 ## Forecast errors:
-## ME: 2.814; MAE: 2.969; RMSE: 3.655
-## sCE: 14.854%; Asymmetry: 87.8%; sMAE: 1.306%; sMSE: 0.026%
-## MASE: 2.492; RMSSE: 2.382; rMAE: 0.958; rRMSE: 0.954

+## ME: 2.807; MAE: 2.966; RMSE: 3.65
+## sCE: 14.815%; Asymmetry: 87.5%; sMAE: 1.305%; sMSE: 0.026%
+## MASE: 2.489; RMSSE: 2.379; rMAE: 0.957; rRMSE: 0.952

 
# Combine the pool of models
 es(BJsales, model=c("CCC","ANN","AAN","AAdN","MNN","MAN","MAdN"),
    h=12, holdout=TRUE)

@@ -570,46 +570,46 @@ 
2023-09-16

 ## 
 ## Number of models combined: 6
 ## Sample size: 138
-## Average number of estimated parameters: 6.4484
-## Average number of degrees of freedom: 131.5516
+## Average number of estimated parameters: 6.6291
+## Average number of degrees of freedom: 131.3709
 ## 
 ## Forecast errors:
-## ME: 2.851; MAE: 2.998; RMSE: 3.694
-## sCE: 15.051%; sMAE: 1.319%; sMSE: 0.026%
-## MASE: 2.517; RMSSE: 2.408; rMAE: 0.967; rRMSE: 0.964

+## ME: 2.83; MAE: 2.983; RMSE: 3.673
+## sCE: 14.936%; sMAE: 1.312%; sMSE: 0.026%
+## MASE: 2.504; RMSSE: 2.394; rMAE: 0.962; rRMSE: 0.959

Now we introduce explanatory variable in ETS:

x <- BJsales.lead

and fit an ETSX model with the exogenous variable first:

es(BJsales, model="ZZZ", h=12, holdout=TRUE,
    xreg=x)
-
## Time elapsed: 1.02 seconds
+
## Time elapsed: 1.01 seconds
 ## Model estimated using es() function: ETSX(AAdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 237.5567
+## Loss function type: likelihood; Loss function value: 237.5664
 ## Persistence vector g (excluding xreg):
 ##  alpha   beta 
-## 0.9498 0.2933 
-## Damping parameter: 0.8812
+## 0.9321 0.3047 
+## Damping parameter: 0.8768
 ## Sample size: 138
 ## Number of estimated parameters: 7
 ## Number of degrees of freedom: 131
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 489.1133 489.9748 509.6041 511.7266 
+## 489.1327 489.9943 509.6235 511.7460 
 ## 
 ## Forecast errors:
-## ME: 2.875; MAE: 2.998; RMSE: 3.7
-## sCE: 15.178%; Asymmetry: 90%; sMAE: 1.319%; sMSE: 0.026%
-## MASE: 2.517; RMSSE: 2.412; rMAE: 0.967; rRMSE: 0.966

+## ME: 2.872; MAE: 2.994; RMSE: 3.696
+## sCE: 15.161%; Asymmetry: 89.9%; sMAE: 1.317%; sMSE: 0.026%
+## MASE: 2.514; RMSSE: 2.409; rMAE: 0.966; rRMSE: 0.965

If we want to check if lagged x can be used for forecasting purposes, we can use xregExpander() function from greybox package:

es(BJsales, model="ZZZ", h=12, holdout=TRUE,
    xreg=xregExpander(x), regressors="use")
-
## Time elapsed: 0.49 seconds
+
## Time elapsed: 0.5 seconds
 ## Model estimated using es() function: ETSX(ANN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 251.4178
+## Loss function type: likelihood; Loss function value: 252.7445
 ## Persistence vector g (excluding xreg):
 ## alpha 
 ##     1 
@@ -619,35 +619,35 @@ 
2023-09-16

 ## Number of degrees of freedom: 132
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 514.8355 515.4767 532.3990 533.9788 
+## 517.4889 518.1301 535.0524 536.6322 
 ## 
 ## Forecast errors:
-## ME: 2.447; MAE: 2.925; RMSE: 3.462
-## sCE: 12.917%; Asymmetry: 73.6%; sMAE: 1.287%; sMSE: 0.023%
-## MASE: 2.455; RMSSE: 2.257; rMAE: 0.943; rRMSE: 0.904

+## ME: 2.738; MAE: 2.983; RMSE: 3.648
+## sCE: 14.453%; Asymmetry: 82.7%; sMAE: 1.312%; sMSE: 0.026%
+## MASE: 2.504; RMSSE: 2.378; rMAE: 0.962; rRMSE: 0.952

We can also construct a model with selected exogenous (based on IC):

es(BJsales, model="ZZZ", h=12, holdout=TRUE,
    xreg=xregExpander(x), regressors="select")
-
## Time elapsed: 0.96 seconds
+
## Time elapsed: 0.99 seconds
 ## Model estimated using es() function: ETS(AMdN)
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 237.7324
+## Loss function type: likelihood; Loss function value: 237.8564
 ## Persistence vector g:
-## alpha  beta 
-## 1.000 0.249 
-## Damping parameter: 0.9034
+##  alpha   beta 
+## 0.9984 0.2184 
+## Damping parameter: 0.9012
 ## Sample size: 138
 ## Number of estimated parameters: 6
 ## Number of degrees of freedom: 132
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 487.4647 488.1060 505.0283 506.6080 
+## 487.7127 488.3540 505.2763 506.8560 
 ## 
 ## Forecast errors:
-## ME: 2.83; MAE: 2.986; RMSE: 3.676
-## sCE: 14.939%; Asymmetry: 87.8%; sMAE: 1.314%; sMSE: 0.026%
-## MASE: 2.507; RMSSE: 2.396; rMAE: 0.963; rRMSE: 0.959

+## ME: 2.888; MAE: 3.029; RMSE: 3.734
+## sCE: 15.244%; Asymmetry: 88.8%; sMAE: 1.332%; sMSE: 0.027%
+## MASE: 2.543; RMSSE: 2.434; rMAE: 0.977; rRMSE: 0.975

Finally, if you work with M or M3 data, and need to test a function on a specific time series, you can use the following simplified call:

diff --git a/inst/doc/gum.html b/inst/doc/gum.html index 2896903..27d54c8 100644 --- a/inst/doc/gum.html +++ b/inst/doc/gum.html @@ -12,7 +12,7 @@ - + gum() - Generalised Univariate Model @@ -340,7 +340,7 @@

gum() - Generalised Univariate Model

Ivan Svetunkov

-

2023-09-16

+

2024-04-01

@@ -354,7 +354,7 @@

2023-09-16

A simple call by default constructs GUM\((1^1,1^m)\), where \(m\) is frequency of the data. So for our example with AirPassengers data, we will have GUM\((1^1,1^{12})\):

gum(AirPassengers, h=18, holdout=TRUE)
-
## Time elapsed: 0.6 seconds
+
## Time elapsed: 0.63 seconds
 ## Model estimated: GUM(1[1],1[12])
 ## Persistence vector g:
 ##        [,1]   [,2]
@@ -381,7 +381,7 @@ 
2023-09-16

 ## MASE: 0.6; sMAE: 5.4%; sMSE: 0.5%; rMAE: 0.216; rRMSE: 0.233

 
But some different orders and lags can be specified. For example:

 
gum(AirPassengers, h=18, holdout=TRUE, orders=c(2,1), lags=c(1,12))

-
## Time elapsed: 0.8 seconds
+
## Time elapsed: 0.83 seconds
 ## Model estimated: GUM(2[1],1[12])
 ## Persistence vector g:
 ##        [,1]   [,2]  [,3]
@@ -415,32 +415,32 @@ 
2023-09-16

 ## Searching for appropriate lags:  —\|/—\|/—\|/We found them!
 ## Searching for appropriate orders:  —\|/—\|/—Orders found.
 ## Reestimating the model. Done!

-
## Time elapsed: 13 seconds
-## Model estimated: GUM(1[1],3[12])
+
## Time elapsed: 13.83 seconds
+## Model estimated: GUM(2[1],2[12])
 ## Persistence vector g:
-##        [,1]   [,2]    [,3]   [,4]
-## [1,] 0.6227 0.2699 -0.0083 0.5413
+##         [,1]   [,2]    [,3]   [,4]
+## [1,] -0.1415 0.6703 -0.0281 0.8222
 ## Transition matrix F: 
 ##        [,1]   [,2]   [,3]   [,4]
-## [1,] 0.9086 0.0001 0.0058 0.0000
-## [2,] 0.4310 0.9214 0.1803 0.0000
-## [3,] 0.2741 0.0084 0.0004 0.9978
-## [4,] 0.0002 0.9940 0.0003 0.0906
+## [1,] 0.0905 0.9642 0.0001 0.0021
+## [2,] 0.4753 0.3107 0.0004 0.0001
+## [3,] 0.1894 0.4591 0.0008 0.9844
+## [4,] 0.9971 0.2650 0.7327 0.4419
 ## Measurement vector w: 1, 1, 1, 1
 ## Initial values were produced using backcasting.
 ## 
-## Loss function type: likelihood; Loss function value: 524.9815
-## Error standard deviation: 10.03
+## Loss function type: likelihood; Loss function value: 532.0226
+## Error standard deviation: 10.5326
 ## Sample size: 144
 ## Number of estimated parameters: 21
 ## Number of provided parameters: 4
 ## Number of degrees of freedom: 123
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 1091.963 1099.537 1154.329 1173.149 
+## 1106.045 1113.619 1168.411 1187.231 
 ## 
 ## 95% parametric prediction interval was constructed

-

+

 
In addition to standard values that other functions accept, GUM
 accepts predefined values for transition matrix, measurement and
 persistence vectors. For example, something more common can be passed to
@@ -448,7 +448,7 @@ 
2023-09-16

 
    transition <- matrix(c(1,0,0,1,1,0,0,0,1),3,3)
     measurement <- c(1,1,1)
     gum(AirPassengers, h=18, holdout=TRUE, orders=c(2,1), lags=c(1,12), transition=transition, measurement=measurement)

-
## Time elapsed: 0.14 seconds
+
## Time elapsed: 0.15 seconds
 ## Model estimated: GUM(2[1],1[12])
 ## Persistence vector g:
 ##        [,1]   [,2]   [,3]
diff --git a/inst/doc/oes.html b/inst/doc/oes.html
index 0dbc9f8..a59992a 100644
--- a/inst/doc/oes.html
+++ b/inst/doc/oes.html
@@ -12,7 +12,7 @@
 
 
 
-
+
 
 oes() - occurrence part of iETS model
 
@@ -340,7 +340,7 @@
 
 
oes() - occurrence part of iETS model

 
Ivan Svetunkov

-
2023-09-16

+
2024-04-01

 
 
 
@@ -526,17 +526,17 @@ 
iETS\(_F\)

 ## Underlying ETS model: oETS[F](MNN)
 ## Vector of initials:
 ##  level 
-## 0.6182 
+## 0.6091 
 ## 
-## Error standard deviation: 1.0937
+## Error standard deviation: 1.0945
 ## Sample size: 110
 ## Number of estimated parameters: 1
 ## Number of degrees of freedom: 109
 ## Information criteria: 
 ##      AIC     AICc      BIC     BICc 
-## 148.2884 148.3254 150.9889 151.0759

+## 149.2137 149.2507 151.9141 152.0012

 
plot(oETSFModel1)

-

+

 
The occurrence part of the model is supported by adam()
 function. For example, here’s how the iETS(M,M,N)\(_F\) can be constructed:

 
adam(y, "MMN", occurrence="fixed", h=10, holdout=TRUE, silent=FALSE)

@@ -544,21 +544,21 @@ 
iETS\(_F\)

 ## Model estimated using adam() function: iETS(MMN)[F]
 ## Occurrence model type: Fixed probability
 ## Distribution assumed in the model: Mixture of Bernoulli and Gamma
-## Loss function type: likelihood; Loss function value: 134.2167
+## Loss function type: likelihood; Loss function value: 131.0314
 ## Persistence vector g:
-## alpha  beta 
-##     0     0 
+##  alpha   beta 
+## 0.0929 0.0000 
 ## 
 ## Sample size: 110
 ## Number of estimated parameters: 6
 ## Number of degrees of freedom: 104
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 426.7219 427.2988 442.9247 439.5802 
+## 421.2765 421.8534 437.4794 434.1348 
 ## 
 ## Forecast errors:
-## Asymmetry: 46.336%; sMSE: 80.876%; rRMSE: 0.618; sPIS: -2007.855%; sCE: 285.432%

-

+## Asymmetry: 65.276%; sMSE: 141.086%; rRMSE: 0.652; sPIS: -3639.762%; sCE: 626.321%

+

 
 

 
iETS\(_O\)

@@ -612,41 +612,41 @@ 
iETS\(_O\)

 ## Underlying ETS model: oETS[O](MMN)
 ## Smoothing parameters:
 ## level trend 
-## 0.001 0.001 
+##     0     0 
 ## Vector of initials:
 ##  level  trend 
-## 0.1215 1.0190 
+## 0.0827 1.0605 
 ## 
-## Error standard deviation: 2.2848
+## Error standard deviation: 2.1412
 ## Sample size: 110
 ## Number of estimated parameters: 4
 ## Number of degrees of freedom: 106
 ## Information criteria: 
 ##      AIC     AICc      BIC     BICc 
-##  93.7540  94.1349 104.5559 105.4512
+## 103.7227 104.1037 114.5247 115.4200 
plot(oETSOModel)
-

+

And here’s the full iETS(M,M,N)\(_O\) model:

adam(y, "MMN", occurrence="o", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.07 seconds
 ## Model estimated using adam() function: iETS(MMN)[O]
 ## Occurrence model type: Odds ratio
 ## Distribution assumed in the model: Mixture of Bernoulli and Gamma
-## Loss function type: likelihood; Loss function value: 134.2167
+## Loss function type: likelihood; Loss function value: 131.0314
 ## Persistence vector g:
-## alpha  beta 
-##     0     0 
+##  alpha   beta 
+## 0.0929 0.0000 
 ## 
 ## Sample size: 110
 ## Number of estimated parameters: 9
 ## Number of degrees of freedom: 101
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 372.1874 372.7644 396.4918 379.0457 
+## 375.7856 376.3625 400.0899 382.6439 
 ## 
 ## Forecast errors:
-## Asymmetry: -65.758%; sMSE: 101.715%; rRMSE: 0.693; sPIS: 2287.89%; sCE: -521.942%
-

+## Asymmetry: -40.488%; sMSE: 107.673%; rRMSE: 0.569; sPIS: 824.807%; sCE: -216.586% +

This should give the same results as before, meaning that we ask explicitly for the adam() function to use the earlier estimated model:

@@ -707,42 +707,42 @@

iETS\(_I\)

## Occurrence state space model estimated: Inverse odds ratio
 ## Underlying ETS model: oETS[I](MMN)
 ## Smoothing parameters:
-##  level  trend 
-## 0.0147 0.0000 
+## level trend 
+## 0.033 0.000 
 ## Vector of initials:
 ##   level   trend 
-## 34.9548  0.9130 
+## 19.9393  0.9158 
 ## 
-## Error standard deviation: 4.0968
+## Error standard deviation: 7.863
 ## Sample size: 110
 ## Number of estimated parameters: 4
 ## Number of degrees of freedom: 106
 ## Information criteria: 
 ##      AIC     AICc      BIC     BICc 
-##  95.8162  96.1972 106.6181 107.5135
+## 106.3568 106.7378 117.1587 118.0541 
plot(oETSIModel)
-

+

And here’s the full iETS(M,M,N)\(_O\) model:

adam(y, "MMN", occurrence="i", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.07 seconds
 ## Model estimated using adam() function: iETS(MMN)
 ## Occurrence model type: Inverse odds ratio
 ## Distribution assumed in the model: Mixture of Bernoulli and Gamma
-## Loss function type: likelihood; Loss function value: 134.2167
+## Loss function type: likelihood; Loss function value: 131.0314
 ## Persistence vector g:
-## alpha  beta 
-##     0     0 
+##  alpha   beta 
+## 0.0929 0.0000 
 ## 
 ## Sample size: 110
 ## Number of estimated parameters: 9
 ## Number of degrees of freedom: 101
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 374.2497 374.8266 398.5540 381.1080 
+## 378.4196 378.9966 402.7240 385.2779 
 ## 
 ## Forecast errors:
-## Asymmetry: -65.77%; sMSE: 101.738%; rRMSE: 0.693; sPIS: 2288.308%; sCE: -522.097%
-

+## Asymmetry: -38.282%; sMSE: 106.822%; rRMSE: 0.567; sPIS: 662.135%; sCE: -190.683% +

Once again, an earlier estimated model can be used in the univariate forecasting functions:

adam(y, "MMN", occurrence=oETSIModel, h=10, holdout=TRUE, silent=FALSE)
@@ -801,20 +801,20 @@

iETS\(_D\)

## Underlying ETS model: oETS[D](MMN) ## Smoothing parameters: ## level trend -## 0.0382 0.0000 +## 0.0983 0.0000 ## Vector of initials: ## level trend -## 0.1960 1.0166 +## 0.1806 1.0002 ## -## Error standard deviation: 1.264 +## Error standard deviation: 2.0949 ## Sample size: 110 ## Number of estimated parameters: 4 ## Number of degrees of freedom: 106 ## Information criteria: ## AIC AICc BIC BICc -## 99.6179 99.9988 110.4198 111.3151 +## 116.1501 116.5310 126.9520 127.8473 
plot(oETSDModel)
-

+

The usage of the model in case of univariate forecasting functions is the same as in the cases of other occurrence models, discussed above:

@@ -823,21 +823,21 @@

iETS\(_D\)

## Model estimated using adam() function: iETS(MMN)[D] ## Occurrence model type: Direct ## Distribution assumed in the model: Mixture of Bernoulli and Gamma -## Loss function type: likelihood; Loss function value: 134.2167 +## Loss function type: likelihood; Loss function value: 131.0314 ## Persistence vector g: -## alpha beta -## 0 0 +## alpha beta +## 0.0929 0.0000 ## ## Sample size: 110 ## Number of estimated parameters: 5 ## Number of degrees of freedom: 105 ## Information criteria: ## AIC AICc BIC BICc -## 370.0513 370.6283 383.5537 384.9097 +## 380.2129 380.7898 393.7153 395.0712 ## ## Forecast errors: -## Asymmetry: -66.356%; sMSE: 102.408%; rRMSE: 0.696; sPIS: 2330.461%; sCE: -529.041% -

+## Asymmetry: -33.861%; sMSE: 104.774%; rRMSE: 0.562; sPIS: 438.341%; sCE: -141.717% +

iETS\(_G\)

@@ -876,9 +876,9 @@

iETS\(_G\)

## Number of degrees of freedom: 104 ## Information criteria: ## AIC AICc BIC BICc -## 98.8576 99.6731 115.0604 116.9771 +## 113.2555 114.0710 129.4583 131.3750 
plot(oETSGModel1)
-

+

The oes() function accepts occurrence="g" and in this case calls for oesg() with the same types of ETS models for both parts:

@@ -892,9 +892,9 @@

iETS\(_G\)

## Number of degrees of freedom: 106 ## Information criteria: ## AIC AICc BIC BICc -## 105.8325 106.2135 116.6344 117.5298 +## 110.2579 110.6389 121.0599 121.9552 
plot(oETSGModel2)
-

+

Finally, the more flexible way to construct iETS model would be to do it in two steps: either using oesg() or oes() and then using the es() with the provided model in @@ -905,21 +905,21 @@

iETS\(_G\)

## Model estimated using es() function: iETS(MMN)[G] ## Occurrence model type: General ## Distribution assumed in the model: Mixture of Bernoulli and Normal -## Loss function type: likelihood; Loss function value: 146.7158 +## Loss function type: likelihood; Loss function value: 129.3672 ## Persistence vector g: ## alpha beta -## 0.0098 0.0000 +## 0.1890 0.0001 ## ## Sample size: 110 ## Number of estimated parameters: 13 ## Number of degrees of freedom: 97 ## Information criteria: ## AIC AICc BIC BICc -## 451.3876 451.9645 486.4938 450.2459 +## 450.9306 451.5075 486.0368 449.7889 ## ## Forecast errors: -## Asymmetry: 99.686%; sMSE: 327.94%; rRMSE: 1.245; sPIS: -9024.098%; sCE: 1602.907% -

+## Asymmetry: 100%; sMSE: 496.019%; rRMSE: 1.222; sPIS: -10867.109%; sCE: 1986.235% +

iETS\(_A\)

@@ -933,21 +933,21 @@

iETS\(_A\)

## Occurrence state space model estimated: Odds ratio
 ## Underlying ETS model: oETS[O](MNN)
 ## Smoothing parameters:
-## level 
-## 0.103 
+##  level 
+## 0.1041 
 ## Vector of initials:
 ##  level 
-## 0.1179 
+## 0.0925 
 ## 
-## Error standard deviation: 2.2021
+## Error standard deviation: 2.6804
 ## Sample size: 110
 ## Number of estimated parameters: 2
 ## Number of degrees of freedom: 108
 ## Information criteria: 
 ##      AIC     AICc      BIC     BICc 
-## 101.8325 101.9447 107.2335 107.4970
+## 106.2579 106.3701 111.6589 111.9225 
plot(oETSAModel)
-

+

The main restriction of the iETS models at the moment (smooth v.2.5.0) is that there is no model selection between the ETS models for the occurrence part. This needs to be done diff --git a/inst/doc/simulate.html b/inst/doc/simulate.html index 9463578..8cf54ab 100644 --- a/inst/doc/simulate.html +++ b/inst/doc/simulate.html @@ -12,7 +12,7 @@ - + Simulate functions of the package @@ -340,7 +340,7 @@

Simulate functions of the package

Ivan Svetunkov

-

2023-09-16

+

2024-04-01

@@ -369,16 +369,16 @@

Exponential Smoothing

We can plot produced data, states or residuals in order to see what was generated. This is done using:

plot(ourSimulation$data)
-

+

If only one time series has been generated, we can use a simpler command in order to plot it:

plot(ourSimulation)
-

+

Now let’s use more complicated model and be more specific, providing persistence vector:

ourSimulation <- sim.es("MAdM", frequency=12, obs=120, phi=0.95, persistence=c(0.1,0.05,0.01))
 plot(ourSimulation)
-

+

High values of smoothing parameters are not advised for models with multiplicative components, because they may lead to explosive data. As for randomizer the default values seem to work fine in the majority of @@ -386,7 +386,7 @@

Exponential Smoothing

(for example, some values taken from some estimated model):

ourSimulation <- sim.es("MAdM", frequency=12, obs=120, phi=0.95, persistence=c(0.1,0.05,0.01), randomizer="rlnorm", meanlog=0, sdlog=0.015)
 plot(ourSimulation)
-

+

It is advised to use lower values for sdlog and sd for models with multiplicative components. Once again, using higher values may lead to data with explosive behaviour.

@@ -395,7 +395,7 @@

Exponential Smoothing

and specify initials in this case:

ourSimulation <- sim.es("MNN", frequency=12, obs=120, probability=0.2, initial=10, persistence=0.1)
 plot(ourSimulation)
-

+

If we want to have several time series generated using the same parameters then we can use nsim parameter:

ourSimulation <- sim.es("MNN", frequency=12, obs=120, probability=0.2, initial=10, persistence=0.1, nsim=50)
@@ -413,7 +413,7 @@

Exponential Smoothing

plot(x) plot(ourData$data[,1]) par(mfcol=c(1,1))
-

+

As we see the level is the same and variance is similar for both series. Achievement unlocked!

@@ -443,12 +443,12 @@

SARIMA

or residuals for each time series in order to see what was generated. Here’s an example:

plot(ourSimulation$data[,5])
-

+

Now let’s use more complicated model. For example, data from SARIMA(0,1,1)(1,0,2)_12 with drift can be generated using:

ourSimulation <- sim.ssarima(orders=list(ar=c(0,1),i=c(1,0),ma=c(1,2)), lags=c(1,12), constant=TRUE, obs=120)
 plot(ourSimulation)
-

+

If we want to provide some specific parameters, then we should follow the structure: from lower lag to lag from lower order to higher order. For example, the same model with predefined MA terms will be:

@@ -457,14 +457,14 @@

SARIMA

## Data generated from: SARIMA(0,1,1)[1](1,0,2)[12] with drift
 ## Number of generated series: 1
 ## AR parameters: 
-##        Lag 12
-## AR(1) -0.0362
+##       Lag 12
+## AR(1) -0.836
 ## MA parameters: 
 ##       Lag 1 Lag 12
 ## MA(1)   0.5    0.2
 ## MA(2)   0.0    0.3
-## Constant value: -14.411
-## True likelihood: -580.5031
+## Constant value: -17.909 +## True likelihood: -514.4043

We can create time series with several frequencies For example, some sort of daily series from SARIMA(1,0,2)(0,1,1)_7(1,0,1)_30 can be generated with a command:

@@ -473,16 +473,16 @@

SARIMA

## Data generated from: SARIMA(1,0,2)[1](0,1,1)[7](1,0,1)[30]
 ## Number of generated series: 1
 ## AR parameters: 
-##        Lag 1  Lag 30
-## AR(1) 0.9188 -0.0975
+##        Lag 1 Lag 30
+## AR(1) 0.1365 -3e-04
 ## MA parameters: 
-##         Lag 1  Lag 7 Lag 30
-## MA(1) -0.0769 -0.093 -0.095
-## MA(2)  0.4108  0.000  0.000
+##        Lag 1   Lag 7 Lag 30
+## MA(1) 0.2225 -0.5478 0.1979
+## MA(2) 0.2635  0.0000 0.0000
 ## Constant value: 0
-## True likelihood: -1820.9788
+## True likelihood: -1631.6521 
plot(ourSimulation)
-

+

sim.ssarima also supports intermittent data, which is defined via probability parameter, similar to sim.es():

@@ -491,16 +491,16 @@

SARIMA

## Data generated from: iSARIMA(1,0,2)[1](0,1,1)[7]
 ## Number of generated series: 1
 ## AR parameters: 
-##         Lag 1
-## AR(1) -0.9167
+##        Lag 1
+## AR(1) 0.0026
 ## MA parameters: 
-##        Lag 1  Lag 7
-## MA(1) 0.8052 0.2275
-## MA(2) 0.5686 0.0000
+##         Lag 1  Lag 7
+## MA(1)  0.8294 0.6553
+## MA(2) -0.1607 0.0000
 ## Constant value: 0
-## True likelihood: -529.1351
+## True likelihood: -588.7486 
plot(ourSimulation)
-

+

Finally we can use simulate() function in a similar manner as with sim.es(). For example:

x <- ts(100 + c(1:100) + rnorm(100,0,15),frequency=12)
@@ -516,7 +516,7 @@ 
SARIMA

 plot(x)
 plot(ourData)
 par(mfcol=c(1,1))
-

+

As we see series demonstrate similarities in dynamics and have similar variances.

@@ -546,26 +546,26 @@

Complex Exponential Smoothing

produced data, states or residuals for each time series in order to see what was generated. Here’s an example:

plot(ourSimulation)
-

+

We can also see a brief summary of our simulated data:

ourSimulation
## Data generated from: CES(n)
 ## Number of generated series: 1
-## Smoothing parameter a: 2.1874+1.082i
-## True likelihood: -563.9533
+## Smoothing parameter a: 1.7167+1.0799i +## True likelihood: -552.2716

We can produce one out of three possible seasonal CES models. For example, Let’s generate data from “Simple CES”, which does not have a level:

ourSimulation <- sim.ces("s",frequency=24, obs=240, nsim=1)
 plot(ourSimulation)
-

+

Now let’s be more creative and mess around with the generated initial values of the previous model. We will make some of them equal to zero and regenerate the data:

ourSimulation$initial[c(1:5,20:24),] <- 0
 ourSimulation <- sim.ces("s",frequency=24, obs=120, nsim=1, initial=ourSimulation$initial, randomizer="rt", df=4)
 plot(ourSimulation)
-

+

The resulting generated series has properties close to the ones that solar irradiation data has: changing amplitude of seasonality without changes in level. We have also chosen a random number generator, based @@ -576,13 +576,13 @@

Complex Exponential Smoothing

10 of such time series:

ourSimulation <- sim.ces("p",b=0.2,frequency=12, obs=240, nsim=10)
 plot(ourSimulation)
-

+

Finally, the most complicated CES model is the one with “full” seasonality, implying that there are two complex exponential smoothing models inside: one for seasonal and the other for non-seasonal part:

ourSimulation <- sim.ces("f",frequency=12, obs=240, nsim=10)
 plot(ourSimulation)
-

+

The generated smoothing parameters may sometimes lead to explosive behaviour and produce meaningless time series. That is why it is advised to use parameters of a model fitted to time series of interest. For @@ -597,7 +597,7 @@

Complex Exponential Smoothing

plot(x) plot(ourData) par(mfcol=c(1,1)) -

+

Generalised Exponential Smoothing

@@ -617,7 +617,7 @@

Generalised Exponential Smoothing

plot(x) plot(ourData) par(mfcol=c(1,1))
-

+

Note that GUM is still an ongoing research and its properties are currently understudied.

@@ -627,13 +627,13 @@

Simple Moving Average

simulate the data for it. Here how it can be done:

ourSimulation <- sim.sma(10,frequency=12, obs=240, nsim=1)
 plot(ourSimulation)
-

+

As usual, you can use simulate function as well:

x <- ts(rnorm(100,100,5), frequency=12)
 ourModel <- sma(x)
 ourData <- simulate(ourModel, nsim=50, obs=1000)
 plot(ourData)
-

+

diff --git a/inst/doc/sma.html b/inst/doc/sma.html index 088443a..89c4d21 100644 --- a/inst/doc/sma.html +++ b/inst/doc/sma.html @@ -12,7 +12,7 @@ - + sma() - Simple Moving Average @@ -340,7 +340,7 @@

sma() - Simple Moving Average

Ivan Svetunkov

-

2023-09-16

+

2024-04-01

@@ -368,14 +368,14 @@

2023-09-16

5472, 8634.7, 5032, 6236, 6356, 9857.8, 6322.2, 7907, 13842.4, 13665.1, 3272), .Tsp = c(1983, 1992.5, 12), class = "ts") sma(y, h=18, silent=FALSE) -
## Order 1 - 2119.2666; Order 58 - 2137.8509; Order 115 - 2166.153
-## Order 1 - 2119.2666; Order 29 - 2110.5061; Order 58 - 2137.8509
-## Order 1 - 2119.2666; Order 15 - 2088.1727; Order 29 - 2110.5061
-## Order 15 - 2088.1727; Order 22 - 2101.2155; Order 29 - 2110.5061
-## Order 15 - 2088.1727; Order 18 - 2093.1735; Order 22 - 2101.2155
-## Order 15 - 2088.1727; Order 16 - 2088.1798; Order 18 - 2093.1735
-## Order 12 - 2087.6726
-
## Time elapsed: 0.01 seconds
+
## Order 1 - 2119.2666; Order 58 - 2135.2066; Order 115 - 2157.6146
+## Order 1 - 2119.2666; Order 29 - 2109.953; Order 58 - 2135.2066
+## Order 1 - 2119.2666; Order 15 - 2088.1914; Order 29 - 2109.953
+## Order 15 - 2088.1914; Order 22 - 2101.1263; Order 29 - 2109.953
+## Order 15 - 2088.1914; Order 18 - 2093.311; Order 22 - 2101.1263
+## Order 15 - 2088.1914; Order 16 - 2088.1965; Order 18 - 2093.311
+## Order 12 - 2087.6591

+
## Time elapsed: 0.03 seconds
 ## Model estimated using sma() function: SMA(12)
 ## Distribution assumed in the model: Normal
 ## Loss function type: MSE; Loss function value: 4403812
diff --git a/inst/doc/smooth.Rmd b/inst/doc/smooth.Rmd
index 8f5e69f..bd2d8a6 100644
--- a/inst/doc/smooth.Rmd
+++ b/inst/doc/smooth.Rmd
@@ -14,7 +14,7 @@ knitr::opts_chunk$set(fig.width=6, fig.height=4, fig.path='Figs/', fig.show='hol
                       warning=FALSE, message=FALSE)
 ```
 
-This vignette explains how to use functions in `smooth` package, what they produce, what each field in outputs and what returned values mean. Underlying statistical models are not discussed here, but if you want to know more about them, then there is a document "[Statistical models underlying functions of 'smooth' package for R](https://github.com/config-i1/smooth/blob/master/vignettes/smooth-Documentation.pdf)". Some of the features of the package are also explained in my [blog](https://forecasting.svetunkov.ru/en/tag/smooth/).
+This vignette explains how to use functions in `smooth` package, what they produce, what each field in outputs and what returned values mean. Underlying statistical models are not discussed here, but if you want to know more about them, then there is a working paper "[Smooth forecasting with the smooth package in R](https://doi.org/10.48550/arXiv.2301.01790)". Some of the features of the package are also explained in my [blog](https://openforecast.org/category/r-en/smooth/). Finally, I have a monograph [Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM)](https://openforecast.org/adam/), which explains the model underlying the majority of forecasting functions in the package.
 
 The package includes the following functions:
 
diff --git a/inst/doc/smooth.html b/inst/doc/smooth.html
index a78d187..4c585a1 100644
--- a/inst/doc/smooth.html
+++ b/inst/doc/smooth.html
@@ -12,7 +12,7 @@
 
 
 
-
+
 
 smooth: forecasting using state-space models
 
@@ -240,16 +240,20 @@
 
smooth: forecasting using state-space
 models

 
Ivan Svetunkov

-
2023-09-16

+
2024-04-01

 
 
 
 
This vignette explains how to use functions in smooth
 package, what they produce, what each field in outputs and what returned
 values mean. Underlying statistical models are not discussed here, but
-if you want to know more about them, then there is a document “Statistical
-models underlying functions of ‘smooth’ package for R”. Some of the
-features of the package are also explained in my blog.

+if you want to know more about them, then there is a working paper “Smooth forecasting with
+the smooth package in R”. Some of the features of the package are
+also explained in my blog. Finally,
+I have a monograph Forecasting
+and Analytics with the Augmented Dynamic Adaptive Model (ADAM),
+which explains the model underlying the majority of forecasting
+functions in the package.

 
The package includes the following functions:

 
    
 
  1. adam() - ADAM - Advanced Dynamic Adaptive
diff --git a/inst/doc/ssarima.html b/inst/doc/ssarima.html
index 5f070bd..88720bc 100644
--- a/inst/doc/ssarima.html
+++ b/inst/doc/ssarima.html
@@ -12,7 +12,7 @@
 
 
 
-
+
 
 ssarima() - State-Space ARIMA
 
@@ -362,7 +362,7 @@
 
 
    ssarima() - State-Space ARIMA
    
 
    Ivan Svetunkov
    
-
    2023-09-16
    
+
    2024-04-01
    
 
 
 
@@ -398,29 +398,29 @@ 
    2023-09-16
    
 
    Some more complicated model can be defined using parameter
 orders the following way:
    
 
    ssarima(AirPassengers, orders=list(ar=c(0,1),i=c(1,0),ma=c(1,1)), lags=c(1,12), h=12)
    
-
    ## Time elapsed: 0.05 seconds
+
    ## Time elapsed: 0.02 seconds
 ## Model estimated: SARIMA(0,1,1)[1](1,0,1)[12]
 ## Matrix of AR terms:
 ##       Lag 12
 ## AR(1)      1
 ## Matrix of MA terms:
-##         Lag 1 Lag 12
-## MA(1) -0.3269 0.1017
+##       Lag 1 Lag 12
+## MA(1) -0.31 0.1851
 ## Initial values were produced using backcasting.
 ## 
-## Loss function type: likelihood; Loss function value: 554.9947
-## Error standard deviation: 11.58
+## Loss function type: likelihood; Loss function value: 558.2012
+## Error standard deviation: 11.8407
 ## Sample size: 144
 ## Number of estimated parameters: 4
 ## Number of degrees of freedom: 140
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 1117.989 1118.277 1129.869 1130.584
    
+## 1124.402 1124.690 1136.282 1136.997
    
 
    This would construct seasonal ARIMA(0,1,1)(1,0,1)\(_{12}\).
    
 
    We could try selecting orders manually, but this can also be done
 automatically via auto.ssarima() function:
    
 
    auto.ssarima(AirPassengers, h=12)
    
-
    ## Time elapsed: 2.4 seconds
+
    ## Time elapsed: 1.55 seconds
 ## Model estimated: SARIMA(0,1,3)[1](0,1,0)[12]
 ## Matrix of MA terms:
 ##         Lag 1
@@ -443,7 +443,7 @@ 
    2023-09-16
    
 because of potential overfitting of first observations when non-zero
 order of AR is selected:
    
 
    auto.ssarima(AirPassengers, h=12, initial="backcasting")
    
-
    ## Time elapsed: 2.4 seconds
+
    ## Time elapsed: 1.55 seconds
 ## Model estimated: SARIMA(0,1,3)[1](0,1,0)[12]
 ## Matrix of MA terms:
 ##         Lag 1
@@ -461,7 +461,7 @@ 
    2023-09-16
    
 ##      AIC     AICc      BIC     BICc 
 ## 1107.841 1108.129 1119.720 1120.435
    
 
    auto.ssarima(AirPassengers, h=12, initial="optimal")
    
-
    ## Time elapsed: 22.85 seconds
+
    ## Time elapsed: 14.43 seconds
 ## Model estimated: SARIMA(0,1,3)[1](0,1,0)[12] with drift
 ## Matrix of MA terms:
 ##         Lag 1
@@ -505,22 +505,22 @@ 
    2023-09-16
    
 ## Model estimated: SARIMAX(0,1,3)[1](0,1,0)[12] with drift
 ## Matrix of MA terms:
 ##         Lag 1
-## MA(1) -0.2541
-## MA(2)  0.1523
-## MA(3) -0.2556
-## Constant value is: 0.4061
+## MA(1) -0.2672
+## MA(2)  0.1403
+## MA(3) -0.2293
+## Constant value is: 0.3824
 ## Initial values were provided by user.
 ## Xreg coefficients were estimated in a normal style
 ## 
-## Loss function type: likelihood; Loss function value: 550.2307
-## Error standard deviation: 11.085
+## Loss function type: likelihood; Loss function value: 550.0097
+## Error standard deviation: 11.068
 ## Sample size: 144
 ## Number of estimated parameters: 1
 ## Number of provided parameters: 19
 ## Number of degrees of freedom: 143
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 1102.461 1102.490 1105.431 1105.501 
+## 1102.019 1102.047 1104.989 1105.059 
 ## 
 ## 95% parametric prediction interval was constructed
    
 
    Finally, we can combine several SARIMA models:
    
@@ -557,25 +557,25 @@ 
    MSARIMA
    
 auto.ssarima() and ssarima(). Here’s just one
 example of what can be done with it:
    
 
    msarima(AirPassengers, orders=list(ar=c(0,0,1),i=c(1,0,0),ma=c(1,1,1)),lags=c(1,6,12),h=12, silent=FALSE)
    
-
    ## Time elapsed: 0.29 seconds
+
    ## Time elapsed: 0.23 seconds
 ## Model estimated using msarima() function: SARIMA(0,1,1)[1](0,0,1)[6](1,0,1)[12]
 ## Distribution assumed in the model: Normal
-## Loss function type: likelihood; Loss function value: 579.5247
+## Loss function type: likelihood; Loss function value: 553.982
 ## ARMA parameters of the model:
 ## AR:
 ## phi1[12] 
-##   0.9841 
+##   0.9943 
 ## MA:
 ##  theta1[1]  theta1[6] theta1[12] 
-##     0.2906    -0.2172    -0.2202 
+##    -0.2988     0.0011    -0.0219 
 ## 
 ## Sample size: 144
 ## Number of estimated parameters: 24
 ## Number of degrees of freedom: 120
 ## Information criteria:
 ##      AIC     AICc      BIC     BICc 
-## 1207.049 1217.133 1278.325 1303.383
    
-
    
+## 1155.964 1166.048 1227.239 1252.297
    
+
    
 
    The forecasts of the two models might differ due to the different
 state space form. The detailed explanation of MSARIMA is given in
 Chapter 9 of ADAM
diff --git a/man/adam.Rd b/man/adam.Rd
index 2b98ff5..d3b8969 100644
--- a/man/adam.Rd
+++ b/man/adam.Rd
@@ -421,8 +421,16 @@ x <- simulate(ourModel)
 }
 \references{
 \itemize{
-\item Svetunkov, I. (2021) Time Series Analysis and Forecasting with ADAM:
-Lancaster, UK. \url{https://openforecast.org/adam/}.
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
+\item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
+}
+
+\itemize{
+\item Svetunkov, I. (2023). Forecasting and Analytics with the Augmented
+Dynamic Adaptive Model (ADAM) (1st ed.). Chapman and Hall/CRC.
+\doi{10.1201/9781003452652}, online version: \url{https://openforecast.org/adam/}.
 }
 
 \itemize{
diff --git a/man/cma.Rd b/man/cma.Rd
index ddbcb0e..391a49f 100644
--- a/man/cma.Rd
+++ b/man/cma.Rd
@@ -50,7 +50,7 @@ If the order is odd, then the function constructs SMA(order) and
 shifts it back in time. Otherwise an AR(order+1) model is constructed
 with the preset parameters:
 
-phi_i = {0.5,1,1,...,0.5} / order
+\deqn{phi_i = {0.5,1,1,...,0.5} / order}
 
 This then corresponds to the centered MA with 0.5 weight for the
 first observation and 0.5 weight for an additional one. e.g. if this is
@@ -74,11 +74,10 @@ summary(ourModel)
 }
 \references{
 \itemize{
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
 \item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
-models and how to use them: \url{https://forecasting.svetunkov.ru/en/tag/smooth/}.
-\item Svetunkov I. (2017). Statistical models underlying functions of 'smooth'
-package for R. Working Paper of Department of Management Science, Lancaster
-University 2017:1, 1-52.
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
 }
 }
 \seealso{
diff --git a/man/es.Rd b/man/es.Rd
index 88b8dae..beec549 100644
--- a/man/es.Rd
+++ b/man/es.Rd
@@ -300,7 +300,7 @@ vignette: \code{vignette("es","smooth")}.
 
 Also, there are posts about the functions of the package smooth on the
 website of Ivan Svetunkov:
-\url{https://forecasting.svetunkov.ru/en/tag/smooth/} - they explain the
+\url{https://openforecast.org/category/r-en/smooth/} - they explain the
 underlying models and how to use the functions.
 }
 \examples{
diff --git a/man/gum.Rd b/man/gum.Rd
index f4eeda0..4d81881 100644
--- a/man/gum.Rd
+++ b/man/gum.Rd
@@ -258,11 +258,10 @@ plot(forecast(ourModel))}
 }
 \references{
 \itemize{
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
 \item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
-models and how to use them: \url{https://forecasting.svetunkov.ru/en/tag/smooth/}.
-\item Svetunkov I. (2017). Statistical models underlying functions of 'smooth'
-package for R. Working Paper of Department of Management Science, Lancaster
-University 2017:1, 1-52.
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
 }
 
 \itemize{
diff --git a/man/msdecompose.Rd b/man/msdecompose.Rd
index 8a51b3f..68a6d56 100644
--- a/man/msdecompose.Rd
+++ b/man/msdecompose.Rd
@@ -44,6 +44,14 @@ ourModel <- msdecompose(AirPassengers, lags=c(12), type="m")
 plot(ourModel)
 plot(forecast(ourModel, model="AAN", h=12))
 
+}
+\references{
+\itemize{
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
+\item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
+}
 }
 \seealso{
 \code{\link[stats]{filter}}
diff --git a/man/sim.gum.Rd b/man/sim.gum.Rd
index a7bd6c0..a50fab9 100644
--- a/man/sim.gum.Rd
+++ b/man/sim.gum.Rd
@@ -102,11 +102,10 @@ simulate(ourModel,nsim=10)
 }
 \references{
 \itemize{
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
 \item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
-models and how to use them: \url{https://forecasting.svetunkov.ru/en/tag/smooth/}.
-\item Svetunkov I. (2017). Statistical models underlying functions of 'smooth'
-package for R. Working Paper of Department of Management Science, Lancaster
-University 2017:1, 1-52.
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
 }
 }
 \seealso{
diff --git a/man/sma.Rd b/man/sma.Rd
index 008d75b..ffe012b 100644
--- a/man/sma.Rd
+++ b/man/sma.Rd
@@ -156,11 +156,10 @@ plot(forecast(ourModel, h=18, interval="empirical"))
 }
 \references{
 \itemize{
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
 \item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
-models and how to use them: \url{https://forecasting.svetunkov.ru/en/tag/smooth/}.
-\item Svetunkov I. (2017). Statistical models underlying functions of 'smooth'
-package for R. Working Paper of Department of Management Science, Lancaster
-University 2017:1, 1-52.
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
 }
 
 \itemize{
diff --git a/man/smooth.Rd b/man/smooth.Rd
index c7f6059..16d2ee2 100644
--- a/man/smooth.Rd
+++ b/man/smooth.Rd
@@ -1,6 +1,5 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/smooth-package.R
-\docType{package}
 \name{smooth}
 \alias{smooth}
 \alias{smooth-package}
@@ -98,11 +97,10 @@ for Time Series Forecasting. Not yet published.
 }
 
 \itemize{
+\item Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790.
+\doi{10.48550/arXiv.2301.01790}.
 \item Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying
-models and how to use them: \url{https://forecasting.svetunkov.ru/en/tag/smooth/}.
-\item Svetunkov I. (2017). Statistical models underlying functions of 'smooth'
-package for R. Working Paper of Department of Management Science, Lancaster
-University 2017:1, 1-52.
+models and how to use them: \url{https://openforecast.org/category/r-en/smooth/}.
 }
 
 \itemize{
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 39140c3..70f9e4a 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -12,8 +12,8 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 // adamFitterWrap
-RcppExport SEXP adamFitterWrap(arma::mat matrixVt, arma::mat& matrixWt, arma::mat& matrixF, arma::vec& vectorG, arma::uvec& lags, arma::umat& profilesObserved, arma::mat& profilesRecent, char const& Etype, char const& Ttype, char const& Stype, unsigned int const& componentsNumberETS, unsigned int const& nSeasonal, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant, arma::vec& vectorYt, arma::vec& vectorOt, bool const& backcast);
-RcppExport SEXP _smooth_adamFitterWrap(SEXP matrixVtSEXP, SEXP matrixWtSEXP, SEXP matrixFSEXP, SEXP vectorGSEXP, SEXP lagsSEXP, SEXP profilesObservedSEXP, SEXP profilesRecentSEXP, SEXP EtypeSEXP, SEXP TtypeSEXP, SEXP StypeSEXP, SEXP componentsNumberETSSEXP, SEXP nSeasonalSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP, SEXP vectorYtSEXP, SEXP vectorOtSEXP, SEXP backcastSEXP) {
+RcppExport SEXP adamFitterWrap(arma::mat matrixVt, arma::mat& matrixWt, arma::mat& matrixF, arma::vec& vectorG, arma::uvec& lags, arma::umat& indexLookupTable, arma::mat& profilesRecent, char const& Etype, char const& Ttype, char const& Stype, unsigned int const& componentsNumberETS, unsigned int const& nSeasonal, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant, arma::vec& vectorYt, arma::vec& vectorOt, bool const& backcast);
+RcppExport SEXP _smooth_adamFitterWrap(SEXP matrixVtSEXP, SEXP matrixWtSEXP, SEXP matrixFSEXP, SEXP vectorGSEXP, SEXP lagsSEXP, SEXP indexLookupTableSEXP, SEXP profilesRecentSEXP, SEXP EtypeSEXP, SEXP TtypeSEXP, SEXP StypeSEXP, SEXP componentsNumberETSSEXP, SEXP nSeasonalSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP, SEXP vectorYtSEXP, SEXP vectorOtSEXP, SEXP backcastSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -22,7 +22,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< arma::mat& >::type matrixF(matrixFSEXP);
     Rcpp::traits::input_parameter< arma::vec& >::type vectorG(vectorGSEXP);
     Rcpp::traits::input_parameter< arma::uvec& >::type lags(lagsSEXP);
-    Rcpp::traits::input_parameter< arma::umat& >::type profilesObserved(profilesObservedSEXP);
+    Rcpp::traits::input_parameter< arma::umat& >::type indexLookupTable(indexLookupTableSEXP);
     Rcpp::traits::input_parameter< arma::mat& >::type profilesRecent(profilesRecentSEXP);
     Rcpp::traits::input_parameter< char const& >::type Etype(EtypeSEXP);
     Rcpp::traits::input_parameter< char const& >::type Ttype(TtypeSEXP);
@@ -35,20 +35,20 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< arma::vec& >::type vectorYt(vectorYtSEXP);
     Rcpp::traits::input_parameter< arma::vec& >::type vectorOt(vectorOtSEXP);
     Rcpp::traits::input_parameter< bool const& >::type backcast(backcastSEXP);
-    rcpp_result_gen = Rcpp::wrap(adamFitterWrap(matrixVt, matrixWt, matrixF, vectorG, lags, profilesObserved, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, vectorYt, vectorOt, backcast));
+    rcpp_result_gen = Rcpp::wrap(adamFitterWrap(matrixVt, matrixWt, matrixF, vectorG, lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, vectorYt, vectorOt, backcast));
     return rcpp_result_gen;
 END_RCPP
 }
 // adamForecasterWrap
-RcppExport SEXP adamForecasterWrap(arma::mat& matrixWt, arma::mat& matrixF, arma::uvec& lags, arma::umat& profilesObserved, arma::mat& profilesRecent, char const& E, char const& T, char const& S, unsigned int const& componentsNumberETS, unsigned int const& nSeasonal, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant, unsigned int const& horizon);
-RcppExport SEXP _smooth_adamForecasterWrap(SEXP matrixWtSEXP, SEXP matrixFSEXP, SEXP lagsSEXP, SEXP profilesObservedSEXP, SEXP profilesRecentSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP componentsNumberETSSEXP, SEXP nSeasonalSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP, SEXP horizonSEXP) {
+RcppExport SEXP adamForecasterWrap(arma::mat& matrixWt, arma::mat& matrixF, arma::uvec& lags, arma::umat& indexLookupTable, arma::mat& profilesRecent, char const& E, char const& T, char const& S, unsigned int const& componentsNumberETS, unsigned int const& nSeasonal, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant, unsigned int const& horizon);
+RcppExport SEXP _smooth_adamForecasterWrap(SEXP matrixWtSEXP, SEXP matrixFSEXP, SEXP lagsSEXP, SEXP indexLookupTableSEXP, SEXP profilesRecentSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP componentsNumberETSSEXP, SEXP nSeasonalSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP, SEXP horizonSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::mat& >::type matrixWt(matrixWtSEXP);
     Rcpp::traits::input_parameter< arma::mat& >::type matrixF(matrixFSEXP);
     Rcpp::traits::input_parameter< arma::uvec& >::type lags(lagsSEXP);
-    Rcpp::traits::input_parameter< arma::umat& >::type profilesObserved(profilesObservedSEXP);
+    Rcpp::traits::input_parameter< arma::umat& >::type indexLookupTable(indexLookupTableSEXP);
     Rcpp::traits::input_parameter< arma::mat& >::type profilesRecent(profilesRecentSEXP);
     Rcpp::traits::input_parameter< char const& >::type E(ESEXP);
     Rcpp::traits::input_parameter< char const& >::type T(TSEXP);
@@ -59,13 +59,13 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< unsigned int const& >::type nXreg(nXregSEXP);
     Rcpp::traits::input_parameter< bool const& >::type constant(constantSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type horizon(horizonSEXP);
-    rcpp_result_gen = Rcpp::wrap(adamForecasterWrap(matrixWt, matrixF, lags, profilesObserved, profilesRecent, E, T, S, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon));
+    rcpp_result_gen = Rcpp::wrap(adamForecasterWrap(matrixWt, matrixF, lags, indexLookupTable, profilesRecent, E, T, S, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon));
     return rcpp_result_gen;
 END_RCPP
 }
 // adamErrorerWrap
-RcppExport SEXP adamErrorerWrap(arma::mat matrixVt, arma::mat matrixWt, arma::mat matrixF, arma::uvec lags, arma::umat profilesObserved, arma::mat profilesRecent, char Etype, char Ttype, char Stype, unsigned int& componentsNumberETS, unsigned int& nSeasonal, unsigned int nArima, unsigned int nXreg, bool constant, unsigned int horizon, arma::vec vectorYt, arma::vec vectorOt);
-RcppExport SEXP _smooth_adamErrorerWrap(SEXP matrixVtSEXP, SEXP matrixWtSEXP, SEXP matrixFSEXP, SEXP lagsSEXP, SEXP profilesObservedSEXP, SEXP profilesRecentSEXP, SEXP EtypeSEXP, SEXP TtypeSEXP, SEXP StypeSEXP, SEXP componentsNumberETSSEXP, SEXP nSeasonalSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP, SEXP horizonSEXP, SEXP vectorYtSEXP, SEXP vectorOtSEXP) {
+RcppExport SEXP adamErrorerWrap(arma::mat matrixVt, arma::mat matrixWt, arma::mat matrixF, arma::uvec lags, arma::umat indexLookupTable, arma::mat profilesRecent, char Etype, char Ttype, char Stype, unsigned int& componentsNumberETS, unsigned int& nSeasonal, unsigned int nArima, unsigned int nXreg, bool constant, unsigned int horizon, arma::vec vectorYt, arma::vec vectorOt);
+RcppExport SEXP _smooth_adamErrorerWrap(SEXP matrixVtSEXP, SEXP matrixWtSEXP, SEXP matrixFSEXP, SEXP lagsSEXP, SEXP indexLookupTableSEXP, SEXP profilesRecentSEXP, SEXP EtypeSEXP, SEXP TtypeSEXP, SEXP StypeSEXP, SEXP componentsNumberETSSEXP, SEXP nSeasonalSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP, SEXP horizonSEXP, SEXP vectorYtSEXP, SEXP vectorOtSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -73,7 +73,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< arma::mat >::type matrixWt(matrixWtSEXP);
     Rcpp::traits::input_parameter< arma::mat >::type matrixF(matrixFSEXP);
     Rcpp::traits::input_parameter< arma::uvec >::type lags(lagsSEXP);
-    Rcpp::traits::input_parameter< arma::umat >::type profilesObserved(profilesObservedSEXP);
+    Rcpp::traits::input_parameter< arma::umat >::type indexLookupTable(indexLookupTableSEXP);
     Rcpp::traits::input_parameter< arma::mat >::type profilesRecent(profilesRecentSEXP);
     Rcpp::traits::input_parameter< char >::type Etype(EtypeSEXP);
     Rcpp::traits::input_parameter< char >::type Ttype(TtypeSEXP);
@@ -86,13 +86,31 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< unsigned int >::type horizon(horizonSEXP);
     Rcpp::traits::input_parameter< arma::vec >::type vectorYt(vectorYtSEXP);
     Rcpp::traits::input_parameter< arma::vec >::type vectorOt(vectorOtSEXP);
-    rcpp_result_gen = Rcpp::wrap(adamErrorerWrap(matrixVt, matrixWt, matrixF, lags, profilesObserved, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon, vectorYt, vectorOt));
+    rcpp_result_gen = Rcpp::wrap(adamErrorerWrap(matrixVt, matrixWt, matrixF, lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype, componentsNumberETS, nSeasonal, nArima, nXreg, constant, horizon, vectorYt, vectorOt));
+    return rcpp_result_gen;
+END_RCPP
+}
+// adamPolynomialiser
+RcppExport SEXP adamPolynomialiser(arma::vec const& B, arma::uvec const& arOrders, arma::uvec const& iOrders, arma::uvec const& maOrders, bool const& arEstimate, bool const& maEstimate, SEXP armaParameters, arma::uvec const& lags);
+RcppExport SEXP _smooth_adamPolynomialiser(SEXP BSEXP, SEXP arOrdersSEXP, SEXP iOrdersSEXP, SEXP maOrdersSEXP, SEXP arEstimateSEXP, SEXP maEstimateSEXP, SEXP armaParametersSEXP, SEXP lagsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::vec const& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< arma::uvec const& >::type arOrders(arOrdersSEXP);
+    Rcpp::traits::input_parameter< arma::uvec const& >::type iOrders(iOrdersSEXP);
+    Rcpp::traits::input_parameter< arma::uvec const& >::type maOrders(maOrdersSEXP);
+    Rcpp::traits::input_parameter< bool const& >::type arEstimate(arEstimateSEXP);
+    Rcpp::traits::input_parameter< bool const& >::type maEstimate(maEstimateSEXP);
+    Rcpp::traits::input_parameter< SEXP >::type armaParameters(armaParametersSEXP);
+    Rcpp::traits::input_parameter< arma::uvec const& >::type lags(lagsSEXP);
+    rcpp_result_gen = Rcpp::wrap(adamPolynomialiser(B, arOrders, iOrders, maOrders, arEstimate, maEstimate, armaParameters, lags));
     return rcpp_result_gen;
 END_RCPP
 }
 // adamRefitterWrap
-RcppExport SEXP adamRefitterWrap(arma::mat matrixYt, arma::mat matrixOt, arma::cube arrayVt, arma::cube arrayF, arma::cube arrayWt, arma::mat matrixG, char const& E, char const& T, char const& S, arma::uvec lags, arma::umat profilesObserved, arma::cube arrayProfilesRecent, unsigned int const& nSeasonal, unsigned int const& componentsNumberETS, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant);
-RcppExport SEXP _smooth_adamRefitterWrap(SEXP matrixYtSEXP, SEXP matrixOtSEXP, SEXP arrayVtSEXP, SEXP arrayFSEXP, SEXP arrayWtSEXP, SEXP matrixGSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP lagsSEXP, SEXP profilesObservedSEXP, SEXP arrayProfilesRecentSEXP, SEXP nSeasonalSEXP, SEXP componentsNumberETSSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP) {
+RcppExport SEXP adamRefitterWrap(arma::mat matrixYt, arma::mat matrixOt, arma::cube arrayVt, arma::cube arrayF, arma::cube arrayWt, arma::mat matrixG, char const& E, char const& T, char const& S, arma::uvec lags, arma::umat indexLookupTable, arma::cube arrayProfilesRecent, unsigned int const& nSeasonal, unsigned int const& componentsNumberETS, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant);
+RcppExport SEXP _smooth_adamRefitterWrap(SEXP matrixYtSEXP, SEXP matrixOtSEXP, SEXP arrayVtSEXP, SEXP arrayFSEXP, SEXP arrayWtSEXP, SEXP matrixGSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP lagsSEXP, SEXP indexLookupTableSEXP, SEXP arrayProfilesRecentSEXP, SEXP nSeasonalSEXP, SEXP componentsNumberETSSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -106,20 +124,20 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< char const& >::type T(TSEXP);
     Rcpp::traits::input_parameter< char const& >::type S(SSEXP);
     Rcpp::traits::input_parameter< arma::uvec >::type lags(lagsSEXP);
-    Rcpp::traits::input_parameter< arma::umat >::type profilesObserved(profilesObservedSEXP);
+    Rcpp::traits::input_parameter< arma::umat >::type indexLookupTable(indexLookupTableSEXP);
     Rcpp::traits::input_parameter< arma::cube >::type arrayProfilesRecent(arrayProfilesRecentSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nSeasonal(nSeasonalSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type componentsNumberETS(componentsNumberETSSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nArima(nArimaSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nXreg(nXregSEXP);
     Rcpp::traits::input_parameter< bool const& >::type constant(constantSEXP);
-    rcpp_result_gen = Rcpp::wrap(adamRefitterWrap(matrixYt, matrixOt, arrayVt, arrayF, arrayWt, matrixG, E, T, S, lags, profilesObserved, arrayProfilesRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant));
+    rcpp_result_gen = Rcpp::wrap(adamRefitterWrap(matrixYt, matrixOt, arrayVt, arrayF, arrayWt, matrixG, E, T, S, lags, indexLookupTable, arrayProfilesRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant));
     return rcpp_result_gen;
 END_RCPP
 }
 // adamReforecasterWrap
-RcppExport SEXP adamReforecasterWrap(arma::cube arrayErrors, arma::cube arrayOt, arma::cube arrayF, arma::cube arrayWt, arma::mat matrixG, char const& E, char const& T, char const& S, arma::uvec& lags, arma::umat const& profilesObserved, arma::cube arrayProfileRecent, unsigned int const& nSeasonal, unsigned int const& componentsNumberETS, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant);
-RcppExport SEXP _smooth_adamReforecasterWrap(SEXP arrayErrorsSEXP, SEXP arrayOtSEXP, SEXP arrayFSEXP, SEXP arrayWtSEXP, SEXP matrixGSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP lagsSEXP, SEXP profilesObservedSEXP, SEXP arrayProfileRecentSEXP, SEXP nSeasonalSEXP, SEXP componentsNumberETSSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP) {
+RcppExport SEXP adamReforecasterWrap(arma::cube arrayErrors, arma::cube arrayOt, arma::cube arrayF, arma::cube arrayWt, arma::mat matrixG, char const& E, char const& T, char const& S, arma::uvec& lags, arma::umat const& indexLookupTable, arma::cube arrayProfileRecent, unsigned int const& nSeasonal, unsigned int const& componentsNumberETS, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant);
+RcppExport SEXP _smooth_adamReforecasterWrap(SEXP arrayErrorsSEXP, SEXP arrayOtSEXP, SEXP arrayFSEXP, SEXP arrayWtSEXP, SEXP matrixGSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP lagsSEXP, SEXP indexLookupTableSEXP, SEXP arrayProfileRecentSEXP, SEXP nSeasonalSEXP, SEXP componentsNumberETSSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -132,20 +150,20 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< char const& >::type T(TSEXP);
     Rcpp::traits::input_parameter< char const& >::type S(SSEXP);
     Rcpp::traits::input_parameter< arma::uvec& >::type lags(lagsSEXP);
-    Rcpp::traits::input_parameter< arma::umat const& >::type profilesObserved(profilesObservedSEXP);
+    Rcpp::traits::input_parameter< arma::umat const& >::type indexLookupTable(indexLookupTableSEXP);
     Rcpp::traits::input_parameter< arma::cube >::type arrayProfileRecent(arrayProfileRecentSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nSeasonal(nSeasonalSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type componentsNumberETS(componentsNumberETSSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nArima(nArimaSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nXreg(nXregSEXP);
     Rcpp::traits::input_parameter< bool const& >::type constant(constantSEXP);
-    rcpp_result_gen = Rcpp::wrap(adamReforecasterWrap(arrayErrors, arrayOt, arrayF, arrayWt, matrixG, E, T, S, lags, profilesObserved, arrayProfileRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant));
+    rcpp_result_gen = Rcpp::wrap(adamReforecasterWrap(arrayErrors, arrayOt, arrayF, arrayWt, matrixG, E, T, S, lags, indexLookupTable, arrayProfileRecent, nSeasonal, componentsNumberETS, nArima, nXreg, constant));
     return rcpp_result_gen;
 END_RCPP
 }
 // adamSimulatorWrap
-RcppExport SEXP adamSimulatorWrap(arma::cube arrayVt, arma::mat matrixErrors, arma::mat matrixOt, arma::cube arrayF, arma::mat matrixWt, arma::mat matrixG, char const& E, char const& T, char const& S, arma::uvec lags, arma::umat profilesObserved, arma::mat profilesRecent, unsigned int const& nSeasonal, unsigned int const& componentsNumber, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant);
-RcppExport SEXP _smooth_adamSimulatorWrap(SEXP arrayVtSEXP, SEXP matrixErrorsSEXP, SEXP matrixOtSEXP, SEXP arrayFSEXP, SEXP matrixWtSEXP, SEXP matrixGSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP lagsSEXP, SEXP profilesObservedSEXP, SEXP profilesRecentSEXP, SEXP nSeasonalSEXP, SEXP componentsNumberSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP) {
+RcppExport SEXP adamSimulatorWrap(arma::cube arrayVt, arma::mat matrixErrors, arma::mat matrixOt, arma::cube arrayF, arma::mat matrixWt, arma::mat matrixG, char const& E, char const& T, char const& S, arma::uvec lags, arma::umat indexLookupTable, arma::mat profilesRecent, unsigned int const& nSeasonal, unsigned int const& componentsNumber, unsigned int const& nArima, unsigned int const& nXreg, bool const& constant);
+RcppExport SEXP _smooth_adamSimulatorWrap(SEXP arrayVtSEXP, SEXP matrixErrorsSEXP, SEXP matrixOtSEXP, SEXP arrayFSEXP, SEXP matrixWtSEXP, SEXP matrixGSEXP, SEXP ESEXP, SEXP TSEXP, SEXP SSEXP, SEXP lagsSEXP, SEXP indexLookupTableSEXP, SEXP profilesRecentSEXP, SEXP nSeasonalSEXP, SEXP componentsNumberSEXP, SEXP nArimaSEXP, SEXP nXregSEXP, SEXP constantSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -159,14 +177,14 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< char const& >::type T(TSEXP);
     Rcpp::traits::input_parameter< char const& >::type S(SSEXP);
     Rcpp::traits::input_parameter< arma::uvec >::type lags(lagsSEXP);
-    Rcpp::traits::input_parameter< arma::umat >::type profilesObserved(profilesObservedSEXP);
+    Rcpp::traits::input_parameter< arma::umat >::type indexLookupTable(indexLookupTableSEXP);
     Rcpp::traits::input_parameter< arma::mat >::type profilesRecent(profilesRecentSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nSeasonal(nSeasonalSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type componentsNumber(componentsNumberSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nArima(nArimaSEXP);
     Rcpp::traits::input_parameter< unsigned int const& >::type nXreg(nXregSEXP);
     Rcpp::traits::input_parameter< bool const& >::type constant(constantSEXP);
-    rcpp_result_gen = Rcpp::wrap(adamSimulatorWrap(arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG, E, T, S, lags, profilesObserved, profilesRecent, nSeasonal, componentsNumber, nArima, nXreg, constant));
+    rcpp_result_gen = Rcpp::wrap(adamSimulatorWrap(arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG, E, T, S, lags, indexLookupTable, profilesRecent, nSeasonal, componentsNumber, nArima, nXreg, constant));
     return rcpp_result_gen;
 END_RCPP
 }
diff --git a/src/adamGeneral.cpp b/src/adamGeneral.cpp
index de4d477..8a1a792 100644
--- a/src/adamGeneral.cpp
+++ b/src/adamGeneral.cpp
@@ -9,7 +9,7 @@ using namespace Rcpp;
 
 // # Fitter for univariate models
 List adamFitter(arma::mat &matrixVt, arma::mat const &matrixWt, arma::mat &matrixF, arma::vec const &vectorG,
-                arma::uvec &lags, arma::umat const &profilesObserved, arma::mat profilesRecent,
+                arma::uvec &lags, arma::umat const &indexLookupTable, arma::mat profilesRecent,
                 char const &E, char const &T, char const &S,
                 unsigned int const &nNonSeasonal, unsigned int const &nSeasonal,
                 unsigned int const &nArima, unsigned int const &nXreg, bool const &constant,
@@ -41,9 +41,10 @@ List adamFitter(arma::mat &matrixVt, arma::mat const &matrixWt, arma::mat &matri
         // Refine the head (in order for it to make sense)
         // This is only needed for ETS(*,Z,*) models, with trend.
         // if(!backcast){
-        for (int i=0; i=lagsModelMax; i=i-1) {
                 /* # Measurement equation and the error term */
-                vecYfit(i-lagsModelMax) = adamWvalue(profilesRecent(profilesObserved.col(i)),
+                vecYfit(i-lagsModelMax) = adamWvalue(profilesRecent(indexLookupTable.col(i)),
                         matrixWt.row(i-lagsModelMax), E, T, S,
                         nETS, nNonSeasonal, nSeasonal, nArima, nXreg, nComponents, constant);
 
@@ -98,28 +99,34 @@ List adamFitter(arma::mat &matrixVt, arma::mat const &matrixWt, arma::mat &matri
                 }
 
                 /* # Transition equation */
-                profilesRecent(profilesObserved.col(i)) = adamFvalue(profilesRecent(profilesObserved.col(i)),
+                profilesRecent(indexLookupTable.col(i)) = adamFvalue(profilesRecent(indexLookupTable.col(i)),
                                matrixF, E, T, S, nETS, nNonSeasonal, nSeasonal, nArima, nComponents, constant) +
-                                   adamGvalue(profilesRecent(profilesObserved.col(i)), matrixF,
+                                   adamGvalue(profilesRecent(indexLookupTable.col(i)), matrixF,
                                               matrixWt.row(i-lagsModelMax), E, T, S,
                                               nETS, nNonSeasonal, nSeasonal, nArima, nXreg, nComponents, constant,
                                               vectorG, vecErrors(i-lagsModelMax));
             }
 
-            // Fill in the head of the series
-            for (int i=lagsModelMax-1; i>=0; i=i-1) {
-                profilesRecent(profilesObserved.col(i)) = adamFvalue(profilesRecent(profilesObserved.col(i)),
-                             matrixF, E, T, S, nETS, nNonSeasonal, nSeasonal, nArima, nComponents, constant);
+            // Fill in the head of the series.
+            // if(nArima==0){
+                for (int i=lagsModelMax-1; i>=0; i=i-1) {
+                    profilesRecent(indexLookupTable.col(i)) = adamFvalue(profilesRecent(indexLookupTable.col(i)),
+                                   matrixF, E, T, S, nETS, nNonSeasonal, nSeasonal, nArima, nComponents, constant);
 
-                matrixVt.col(i) = profilesRecent(profilesObserved.col(i));
-            }
+                    matrixVt.col(i) = profilesRecent(indexLookupTable.col(i));
+                }
+            // }
 
             // Change back the specific element in the state vector
             if(T=='A'){
                 profilesRecent(1) = -profilesRecent(1);
+                // Write down correct initials
+                matrixVt.col(0) = profilesRecent(indexLookupTable.col(0));
             }
             else if(T=='M'){
                 profilesRecent(1) = 1/profilesRecent(1);
+                // Write down correct initials
+                matrixVt.col(0) = profilesRecent(indexLookupTable.col(0));
             }
         }
     }
@@ -131,7 +138,7 @@ List adamFitter(arma::mat &matrixVt, arma::mat const &matrixWt, arma::mat &matri
 /* # Wrapper for fitter */
 // [[Rcpp::export]]
 RcppExport SEXP adamFitterWrap(arma::mat matrixVt, arma::mat &matrixWt, arma::mat &matrixF, arma::vec &vectorG,
-                               arma::uvec &lags, arma::umat &profilesObserved, arma::mat &profilesRecent,
+                               arma::uvec &lags, arma::umat &indexLookupTable, arma::mat &profilesRecent,
                                char const &Etype, char const &Ttype, char const &Stype,
                                unsigned int const &componentsNumberETS, unsigned int const &nSeasonal,
                                unsigned int const &nArima, unsigned int const &nXreg, bool const &constant,
@@ -140,30 +147,30 @@ RcppExport SEXP adamFitterWrap(arma::mat matrixVt, arma::mat &matrixWt, arma::ma
     unsigned int nNonSeasonal = componentsNumberETS - nSeasonal;
 
     return wrap(adamFitter(matrixVt, matrixWt, matrixF, vectorG,
-                           lags, profilesObserved, profilesRecent, Etype, Ttype, Stype,
+                           lags, indexLookupTable, profilesRecent, Etype, Ttype, Stype,
                            nNonSeasonal, nSeasonal, nArima, nXreg, constant,
                            vectorYt, vectorOt, backcast));
 }
 
 /* # Function produces the point forecasts for the specified model */
 arma::vec adamForecaster(arma::mat const &matrixWt, arma::mat const &matrixF,
-                         arma::uvec lags, arma::umat const &profilesObserved, arma::mat profilesRecent,
+                         arma::uvec lags, arma::umat const &indexLookupTable, arma::mat profilesRecent,
                          char const &E, char const &T, char const &S,
                          unsigned int const &nNonSeasonal, unsigned int const &nSeasonal,
                          unsigned int const &nArima, unsigned int const &nXreg, bool const &constant,
                          unsigned int const &horizon){
     // unsigned int lagslength = lags.n_rows;
     unsigned int nETS = nNonSeasonal + nSeasonal;
-    unsigned int nComponents = profilesObserved.n_rows;
+    unsigned int nComponents = indexLookupTable.n_rows;
 
     arma::vec vecYfor(horizon, arma::fill::zeros);
 
     /* # Fill in the new xt matrix using F. Do the forecasts. */
     for (unsigned int i=0; i(armaParameters);
+    }
+
+// Form matrices with parameters, that are then used for polynomial multiplication
+    arma::mat arParameters(max(arOrders % lags)+1, arOrders.n_elem, arma::fill::zeros);
+    arma::mat iParameters(max(iOrders % lags)+1, iOrders.n_elem, arma::fill::zeros);
+    arma::mat maParameters(max(maOrders % lags)+1, maOrders.n_elem, arma::fill::zeros);
+
+    arParameters.row(0).fill(1);
+    iParameters.row(0).fill(1);
+    maParameters.row(0).fill(1);
+
+    int lagsModelMax = max(lags);
+
+    int nParam = 0;
+    int armanParam = 0;
+    for(unsigned int i=0; i1){
+            for(unsigned int j=1; j
 #include 
 #include 
-// [[Rcpp::depends(RcppArmadillo)]]
 
 using namespace Rcpp;
 
diff --git a/src/adamRefitter.cpp b/src/adamRefitter.cpp
index eac04fa..36787d7 100644
--- a/src/adamRefitter.cpp
+++ b/src/adamRefitter.cpp
@@ -11,7 +11,7 @@ using namespace Rcpp;
 List adamRefitter(arma::mat const &matrixYt, arma::mat const &matrixOt, arma::cube &arrayVt, arma::cube const &arrayF,
                   arma::cube const &arrayWt, arma::mat const &matrixG,
                   char const &E, char const &T, char const &S, arma::uvec &lags,
-                  arma::umat const &profilesObserved, arma::cube arrayProfilesRecent,
+                  arma::umat const &indexLookupTable, arma::cube arrayProfilesRecent,
                   unsigned int const &nNonSeasonal, unsigned int const &nSeasonal,
                   unsigned int const &nArima, unsigned int const &nXreg, bool const &constant) {
 
@@ -28,15 +28,15 @@ List adamRefitter(arma::mat const &matrixYt, arma::mat const &matrixOt, arma::cu
     for(unsigned int i=0; i1e+100)){
             //     matrixVt.col(j) = matrixVt(lagrows);
             // }
-            matrixVt.col(j) = profilesRecent(profilesObserved.col(j-lagsModelMax));
+            matrixVt.col(j) = profilesRecent(indexLookupTable.col(j-lagsModelMax));
         }
         arrayVt.slice(i) = matrixVt;
     }
@@ -84,13 +84,13 @@ List adamSimulator(arma::cube &arrayVt, arma::mat const &matrixErrors, arma::mat
 RcppExport SEXP adamSimulatorWrap(arma::cube arrayVt, arma::mat matrixErrors, arma::mat matrixOt,
                                   arma::cube arrayF, arma::mat matrixWt, arma::mat matrixG,
                                   char const &E, char const &T, char const &S, arma::uvec lags,
-                                  arma::umat profilesObserved, arma::mat profilesRecent,
+                                  arma::umat indexLookupTable, arma::mat profilesRecent,
                                   unsigned int const &nSeasonal, unsigned int const &componentsNumber,
                                   unsigned int const &nArima, unsigned int const &nXreg, bool const &constant){
 
     unsigned int nNonSeasonal = componentsNumber - nSeasonal;
 
     return wrap(adamSimulator(arrayVt, matrixErrors, matrixOt, arrayF, matrixWt, matrixG,
-                              E, T, S, lags, profilesObserved, profilesRecent,
+                              E, T, S, lags, indexLookupTable, profilesRecent,
                               nNonSeasonal, nSeasonal, nArima, nXreg, constant));
 }
diff --git a/src/ssGeneral.cpp b/src/ssGeneral.cpp
index 74a2630..747af30 100644
--- a/src/ssGeneral.cpp
+++ b/src/ssGeneral.cpp
@@ -6,23 +6,6 @@
 
 using namespace Rcpp;
 
-/* # Function allows to multiply polinomails */
-arma::vec polyMult(arma::vec const &poly1, arma::vec const &poly2){
-
-    int poly1Nonzero = arma::as_scalar(find(poly1,1,"last"));
-    int poly2Nonzero = arma::as_scalar(find(poly2,1,"last"));
-
-    arma::vec poly3(poly1Nonzero + poly2Nonzero + 1, arma::fill::zeros);
-
-    for(int i = 0; i <= poly1Nonzero; ++i){
-        for(int j = 0; j <= poly2Nonzero; ++j){
-            poly3(i+j) += poly1(i) * poly2(j);
-        }
-    }
-
-    return poly3;
-}
-
 /* # Function returns value of CDF-based likelihood function for the whole series */
 double cdf(arma::vec const &vecYt, arma::vec const &vecYfit, double const &errorSD, char const &E){
 
diff --git a/src/ssGeneral.h b/src/ssGeneral.h
index 8d0e273..c91c47f 100644
--- a/src/ssGeneral.h
+++ b/src/ssGeneral.h
@@ -1,7 +1,6 @@
 #include 
 #include 
 #include 
-// [[Rcpp::depends(RcppArmadillo)]]
 
 using namespace Rcpp;
 
@@ -46,6 +45,23 @@ inline arma::mat errorvf(arma::mat yact, arma::mat yfit, char const &E){
     }
 }
 
+/* # Function allows to multiply polinomails */
+inline arma::vec polyMult(arma::vec const &poly1, arma::vec const &poly2){
+
+    int poly1Nonzero = arma::as_scalar(find(poly1,1,"last"));
+    int poly2Nonzero = arma::as_scalar(find(poly2,1,"last"));
+
+    arma::vec poly3(poly1Nonzero + poly2Nonzero + 1, arma::fill::zeros);
+
+    for(int i = 0; i <= poly1Nonzero; ++i){
+        for(int j = 0; j <= poly2Nonzero; ++j){
+            poly3(i+j) += poly1(i) * poly2(j);
+        }
+    }
+
+    return poly3;
+}
+
 /* # Function returns value of w() -- y-fitted -- used in the measurement equation */
 inline double wvalue(arma::vec const &vecVt, arma::rowvec const &rowvecW, char const &E, char const &T, char const &S,
               arma::rowvec const &rowvecXt, arma::vec const &vecAt){
diff --git a/vignettes/smooth.Rmd b/vignettes/smooth.Rmd
index 8f5e69f..bd2d8a6 100644
--- a/vignettes/smooth.Rmd
+++ b/vignettes/smooth.Rmd
@@ -14,7 +14,7 @@ knitr::opts_chunk$set(fig.width=6, fig.height=4, fig.path='Figs/', fig.show='hol
                       warning=FALSE, message=FALSE)
 ```
 
-This vignette explains how to use functions in `smooth` package, what they produce, what each field in outputs and what returned values mean. Underlying statistical models are not discussed here, but if you want to know more about them, then there is a document "[Statistical models underlying functions of 'smooth' package for R](https://github.com/config-i1/smooth/blob/master/vignettes/smooth-Documentation.pdf)". Some of the features of the package are also explained in my [blog](https://forecasting.svetunkov.ru/en/tag/smooth/).
+This vignette explains how to use functions in `smooth` package, what they produce, what each field in outputs and what returned values mean. Underlying statistical models are not discussed here, but if you want to know more about them, then there is a working paper "[Smooth forecasting with the smooth package in R](https://doi.org/10.48550/arXiv.2301.01790)". Some of the features of the package are also explained in my [blog](https://openforecast.org/category/r-en/smooth/). Finally, I have a monograph [Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM)](https://openforecast.org/adam/), which explains the model underlying the majority of forecasting functions in the package.
 
 The package includes the following functions: