es()
is a part of smooth
package and is a wrapper for the ADAM
function with distribution="dnorm"
. It implements
Exponential Smoothing in the ETS form, selecting the most appropriate
model among 30 possible ones.
We will use some of the functions of the greybox
package
in this vignette for demonstrational purposes.
Let’s load the necessary packages:
require(smooth)
require(greybox)
The simplest call for the es()
function is:
<- es(BJsales, h=12, holdout=TRUE, silent=FALSE) ourModel
## Forming the pool of models based on... ANN , AAN , Estimation progress: 33 %44 %56 %67 %78 %89 %100 %... Done!
ourModel
## 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.8564
## Persistence vector g:
## 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.7127 488.3540 505.2763 506.8560
##
## Forecast errors:
## 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 brief information about the model, which contains:
holdout=TRUE
).The function has also produced a graph with actual values, fitted values and point forecasts.
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.555
## Persistence vector g:
## 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
## 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"
This handy function also works with ets()
from forecast
package.
If we need actual values from the model, we can use
actuals()
method from greybox
package:
actuals(ourModel)
## Time Series:
## Start = 1
## End = 138
## Frequency = 1
## [1] 200.1 199.5 199.4 198.9 199.0 200.2 198.6 200.0 200.3 201.2 201.6 201.5
## [13] 201.5 203.5 204.9 207.1 210.5 210.5 209.8 208.8 209.5 213.2 213.7 215.1
## [25] 218.7 219.8 220.5 223.8 222.8 223.8 221.7 222.3 220.8 219.4 220.1 220.6
## [37] 218.9 217.8 217.7 215.0 215.3 215.9 216.7 216.7 217.7 218.7 222.9 224.9
## [49] 222.2 220.7 220.0 218.7 217.0 215.9 215.8 214.1 212.3 213.9 214.6 213.6
## [61] 212.1 211.4 213.1 212.9 213.3 211.5 212.3 213.0 211.0 210.7 210.1 211.4
## [73] 210.0 209.7 208.8 208.8 208.8 210.6 211.9 212.8 212.5 214.8 215.3 217.5
## [85] 218.8 220.7 222.2 226.7 228.4 233.2 235.7 237.1 240.6 243.8 245.3 246.0
## [97] 246.3 247.7 247.6 247.8 249.4 249.0 249.9 250.5 251.5 249.0 247.6 248.8
## [109] 250.4 250.7 253.0 253.7 255.0 256.2 256.0 257.4 260.4 260.0 261.3 260.4
## [121] 261.6 260.8 259.8 259.0 258.9 257.4 257.7 257.9 257.4 257.3 257.6 258.9
## [133] 257.8 257.7 257.2 257.5 256.8 257.5
We can also use persistence or initials only from the model to construct the other one:
# Provided initials
es(BJsales, model=modelType(ourModel),
h=12, holdout=FALSE,
initial=ourModel$initial)
## 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.2509
## Persistence vector g:
## alpha beta
## 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.5019 518.7777 530.5444 531.2355
# Provided persistence
es(BJsales, model=modelType(ourModel),
h=12, holdout=FALSE,
persistence=ourModel$persistence)
## 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: 256.1161
## Persistence vector g:
## 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
## 520.2321 520.5080 532.2747 532.9658
or provide some arbitrary values:
es(BJsales, model=modelType(ourModel),
h=12, holdout=FALSE,
initial=200)
## 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.3409
## Persistence vector g:
## alpha beta
## 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.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.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.5365 0.0002
##
## Sample size: 138
## Number of estimated parameters: 4
## Number of degrees of freedom: 134
## Information criteria:
## AIC AICc BIC BICc
## 1021.898 1022.199 1033.607 1034.348
##
## Forecast errors:
## 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.23 seconds
## Model estimated: ETS(CCN)
## Loss function type: likelihood
##
## Number of models combined: 10
## Sample size: 138
## Average number of estimated parameters: 6.2884
## Average number of degrees of freedom: 131.7116
##
## Forecast errors:
## 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.11 seconds
## Model estimated using es() function: ETS(AAdN)
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 238.0241
## Persistence vector g:
## 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.0481 488.6893 505.6116 507.1914
##
## Forecast errors:
## 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)
## Time elapsed: 0.11 seconds
## Model estimated: ETS(CCN)
## Loss function type: likelihood
##
## Number of models combined: 6
## Sample size: 138
## Average number of estimated parameters: 6.6291
## Average number of degrees of freedom: 131.3709
##
## Forecast errors:
## 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:
<- BJsales.lead x
and fit an ETSX model with the exogenous variable first:
es(BJsales, model="ZZZ", h=12, holdout=TRUE,
xreg=x)
## 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.5664
## Persistence vector g (excluding xreg):
## alpha beta
## 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.1327 489.9943 509.6235 511.7460
##
## Forecast errors:
## 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.5 seconds
## Model estimated using es() function: ETSX(ANN)
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 252.7445
## Persistence vector g (excluding xreg):
## alpha
## 1
##
## Sample size: 138
## Number of estimated parameters: 6
## Number of degrees of freedom: 132
## Information criteria:
## AIC AICc BIC BICc
## 517.4889 518.1301 535.0524 536.6322
##
## Forecast errors:
## 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.99 seconds
## Model estimated using es() function: ETS(AMdN)
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 237.8564
## Persistence vector g:
## 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.7127 488.3540 505.2763 506.8560
##
## Forecast errors:
## 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:
es(Mcomp::M3$N2457, silent=FALSE)
This command has taken the data, split it into in-sample and holdout and produced the forecast of appropriate length to the holdout.