Smooth package

Description

Package contains functions implementing Single Source of Error state space models for purposes of time series analysis and forecasting.

Details

The following functions are included in the package:

• es - Exponential Smoothing in Single Source of Errors State Space form.

• ces - Complex Exponential Smoothing.

• gum - Generalised Exponential Smoothing.

• ssarima - SARIMA in state space framework.

• auto.ces - Automatic selection between seasonal and non-seasonal CES.

• auto.ssarima - Automatic selection of ARIMA orders.

• sma - Simple Moving Average in state space form.

• smoothCombine - the function that combines forecasts from es(), ces(), gum(), ssarima() and sma() functions.

• cma - Centered Moving Average. This is for smoothing time series, not for forecasting.

• sim.es - simulate time series using ETS as a model.

• sim.ces - simulate time series using CES as a model.

• sim.ssarima - simulate time series using SARIMA as a model.

• sim.gum - simulate time series using GUM as a model.

• sim.sma - simulate time series using SMA.

• sim.oes - simulate time series based on occurrence part of ETS model.

• oes - occurrence part of the intermittent state space model.

There are also several methods implemented in the package for the classes "smooth" and "smooth.sim":

• orders - extracts orders of the fitted model.

• lags - extracts lags of the fitted model.

• modelType - extracts type of the fitted model.

• forecast - produces forecast using provided model.

• multicov - returns covariance matrix of multiple steps ahead forecast errors.

• pls - returns Prediction Likelihood Score.

• nparam - returns number of the estimated parameters.

• fitted - extracts fitted values from provided model.

• getResponse - returns actual values from the provided model.

• residuals - extracts residuals of provided model.

• plot - plots either states of the model or produced forecast (depending on what object is passed).

• simulate - uses sim functions (sim.es, sim.ces, sim.ssarima, sim.gum, sim.sma and sim.oes) in order to simulate data using the provided object.

• summary - provides summary of the object.

• AICc, BICc - return, guess what...

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., Boylan, J.E., 2023a. iETS: State Space Model for Intermittent Demand Forecastings. International Journal of Production Economics. 109013. doi:10.1016/j.ijpe.2023.109013

• Teunter R., Syntetos A., Babai Z. (2011). Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214, 606-615.

• Croston, J. (1972) Forecasting and stock control for intermittent demands. Operational Research Quarterly, 23(3), 289-303.

• Svetunkov, I., Kourentzes, N., & Ord, J. K. (2022). Complex exponential smoothing. Naval Research Logistics, 69(8), 1108–1123. https://doi.org/10.1002/nav.22074

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

• Kolassa, S. (2011) Combining exponential smoothing forecasts using Akaike weights. International Journal of Forecasting, 27, pp 238 - 251.

• Svetunkov, I., Boylan, J.E., 2023b. Staying Positive: Challenges and Solutions in Using Pure Multiplicative ETS Models. IMA Journal of Management Mathematics. p. 403-425. doi:10.1093/imaman/dpad028

• Taylor, J.W. and Bunn, D.W. (1999) A Quantile Regression Approach to Generating Prediction Intervals. Management Science, Vol 45, No 2, pp 225-237.

• Lichtendahl Kenneth C., Jr., Grushka-Cockayne Yael, Winkler Robert L., (2013) Is It Better to Average Probabilities or Quantiles? Management Science 59(7):1594-1611. DOI: doi:10.1287/mnsc.1120.1667

See Also

forecast, es, ssarima, ces, gum

Examples

y <- ts(rnorm(100,10,3), frequency=12)

adam(y, h=20, holdout=TRUE)
es(y, h=20, holdout=TRUE)
gum(y, h=20, holdout=TRUE)
auto.ces(y, h=20, holdout=TRUE)
auto.ssarima(y, h=20, holdout=TRUE)

Error measures for an estimated model

Description

Function produces error measures for the provided object and the holdout values of the response variable. Note that instead of parameters x, test, the function accepts the vector of values in holdout. Also, the parameters d and D are not supported - MASE is always calculated via division by first differences.

Usage

## S3 method for class 'smooth'
accuracy(object, holdout = NULL, ...)

## S3 method for class 'smooth.forecast'
accuracy(object, holdout = NULL, ...)

Arguments

Details

The function is a wrapper for the measures function and is implemented for convenience.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

Examples

y <- rnorm(100, 100, 10)
ourModel <- adam(y, holdout=TRUE, h=10)
accuracy(ourModel)

ADAM is Augmented Dynamic Adaptive Model

Description

Function constructs an advanced Single Source of Error model, based on ETS taxonomy and ARIMA elements

Usage

adam(data, model = "ZXZ", lags = c(frequency(data)), orders = list(ar =
  c(0), i = c(0), ma = c(0), select = FALSE), constant = FALSE,
  formula = NULL, regressors = c("use", "select", "adapt"),
  occurrence = c("none", "auto", "fixed", "general", "odds-ratio",
  "inverse-odds-ratio", "direct"), distribution = c("default", "dnorm",
  "dlaplace", "ds", "dgnorm", "dlnorm", "dinvgauss", "dgamma"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "LASSO", "RIDGE", "MSEh",
  "TMSE", "GTMSE", "MSCE"), outliers = c("ignore", "use", "select"),
  level = 0.99, h = 0, holdout = FALSE, persistence = NULL,
  phi = NULL, initial = c("backcasting", "optimal", "two-stage",
  "complete"), arma = NULL, ic = c("AICc", "AIC", "BIC", "BICc"),
  bounds = c("usual", "admissible", "none"), silent = TRUE, ...)

## S3 method for class 'adam'
simulate(object, nsim = 1, seed = NULL,
  obs = nobs(object), ...)

auto.adam(data, model = "ZXZ", lags = c(frequency(data)),
  orders = list(ar = c(3, 3), i = c(2, 1), ma = c(3, 3), select = TRUE),
  formula = NULL, regressors = c("use", "select", "adapt"),
  occurrence = c("none", "auto", "fixed", "general", "odds-ratio",
  "inverse-odds-ratio", "direct"), distribution = c("dnorm", "dlaplace",
  "ds", "dgnorm", "dlnorm", "dinvgauss", "dgamma"), outliers = c("ignore",
  "use", "select"), level = 0.99, h = 0, holdout = FALSE,
  persistence = NULL, phi = NULL, initial = c("backcasting", "optimal",
  "two-stage", "complete"), arma = NULL, ic = c("AICc", "AIC", "BIC",
  "BICc"), bounds = c("usual", "admissible", "none"), silent = TRUE,
  parallel = FALSE, ...)

## S3 method for class 'adam'
sm(object, model = "YYY", lags = NULL, orders = list(ar =
  c(0), i = c(0), ma = c(0), select = FALSE), constant = FALSE,
  formula = NULL, regressors = c("use", "select", "adapt"), data = NULL,
  persistence = NULL, phi = NULL, initial = c("optimal", "backcasting"),
  arma = NULL, ic = c("AICc", "AIC", "BIC", "BICc"), bounds = c("usual",
  "admissible", "none"), silent = TRUE, ...)

Arguments

Details

Function estimates ADAM in a form of the Single Source of Error state space model of the following type:

y_{t} = o_t (w(v_{t-l}) + h(x_t, a_{t-1}) + r(v_{t-l}) \epsilon_{t})

v_{t} = f(v_{t-l}, a_{t-1}) + g(v_{t-l}, a_{t-1}, x_{t}) \epsilon_{t}

Where o_{t} is the Bernoulli distributed random variable (in case of normal data it equals to 1 for all observations), v_{t} is the state vector and l is the vector of lags, x_t is the vector of exogenous variables. w(.) is the measurement function, r(.) is the error function, f(.) is the transition function, g(.) is the persistence function and a_t is the vector of parameters for exogenous variables. Finally, \epsilon_{t} is the error term.

The implemented model allows introducing several seasonal states and supports intermittent data via the occurrence variable.

The error term \epsilon_t can follow different distributions, which are regulated via the distribution parameter. This includes:

1. default - Normal distribution is used for the Additive error models, Gamma is used for the Multiplicative error models.

2. dnorm - Normal distribution,

3. dlaplace - Laplace distribution,

4. ds - S distribution,

5. dgnorm - Generalised Normal distribution,

6. dlnorm - Log-Normal distribution,

7. dgamma - Gamma distribution,

8. dinvgauss - Inverse Gaussian distribution,

For some more information about the model and its implementation, see the vignette: vignette("adam","smooth"). The more detailed explanation of ADAM is provided by Svetunkov (2021).

The function auto.adam() tries out models with the specified distributions and returns the one with the most suitable one based on selected information criterion.

sm.adam method estimates the scale model for the already estimated adam. In order for ADAM to take the SM model into account, the latter needs to be recorded in the former, amending the likelihood and the number of degrees of freedom. This can be done using implant method.

Value

Object of class "adam" is returned. It contains the list of the following values:

• model - the name of the constructed model,

• timeElapsed - the time elapsed for the estimation of the model,

• data - the in-sample part of the data used for the training of the model. Includes the actual values in the first column,

• holdout - the holdout part of the data, excluded for purposes of model evaluation,

• fitted - the vector of fitted values,

• residuals - the vector of residuals,

• forecast - the point forecast for h steps ahead (by default NA is returned). NOTE that these do not always correspond to the conditional expectations for ETS models. See ADAM textbook, Section 6.4. for details (https://openforecast.org/adam/ETSTaxonomyMaths.html),

• states - the matrix of states with observations in rows and states in columns,

• persisten - the vector of smoothing parameters,

• phi - the value of damping parameter,

• transition - the transition matrix,

• measurement - the measurement matrix with observations in rows and state elements in columns,

• initial - the named list of initial values, including level, trend, seasonal, ARIMA and xreg components,

• initialEstimated - the named vector, defining which of the initials were estimated in the model,

• initialType - the type of initialisation used (backcasting/optimal/two-stage/complete/provided),

• orders - the orders of ARIMA used in the estimation,

• constant - the value of the constant (if it was included),

• arma - the list of AR / MA parameters used in the model,

• nParam - the matrix of the estimated / provided parameters,

• occurrence - the oes model used for the occurrence part of the model,

• formula - the formula used for the explanatory variables expansion,

• loss - the type of loss function used in the estimation,

• lossValue - the value of that loss function,

• logLik - the value of the log-likelihood,

• distribution - the distribution function used in the calculation of the likelihood,

• scale - the value of the scale parameter,

• lambda - the value of the parameter used in LASSO / dalaplace / dt,

• B - the vector of all estimated parameters,

• lags - the vector of lags used in the model construction,

• lagsAll - the vector of the internal lags used in the model,

• profile - the matrix with the profile used in the construction of the model,

• profileInitial - the matrix with the initial profile (for the before the sample values),

• call - the call used in the evaluation,

• bounds - the type of bounds used in the process,

• res - result of the model estimation, the output of the nloptr() function, explaining how optimisation went,

• other - the list with other parameters, such as shape for distributions or ARIMA polynomials.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

• 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: https://openforecast.org/adam/.

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., Boylan, J.E., 2023a. iETS: State Space Model for Intermittent Demand Forecastings. International Journal of Production Economics. 109013. doi:10.1016/j.ijpe.2023.109013

• Teunter R., Syntetos A., Babai Z. (2011). Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214, 606-615.

• Croston, J. (1972) Forecasting and stock control for intermittent demands. Operational Research Quarterly, 23(3), 289-303.

• Kolassa, S. (2011) Combining exponential smoothing forecasts using Akaike weights. International Journal of Forecasting, 27, pp 238 - 251.

• Svetunkov, I., Boylan, J.E., 2023b. Staying Positive: Challenges and Solutions in Using Pure Multiplicative ETS Models. IMA Journal of Management Mathematics. p. 403-425. doi:10.1093/imaman/dpad028

• Taylor, J.W. and Bunn, D.W. (1999) A Quantile Regression Approach to Generating Prediction Intervals. Management Science, Vol 45, No 2, pp 225-237.

• Lichtendahl Kenneth C., Jr., Grushka-Cockayne Yael, Winkler Robert L., (2013) Is It Better to Average Probabilities or Quantiles? Management Science 59(7):1594-1611. DOI: doi:10.1287/mnsc.1120.1667

See Also

es, msarima

Examples

### The main examples are provided in the adam vignette, check it out via:
## Not run: vignette("adam","smooth")

# Model selection using a specified pool of models
ourModel <- adam(rnorm(100,100,10), model=c("ANN","ANA","AAA"), lags=c(5,10))
adamSummary <- summary(ourModel)
xtable(adamSummary)

forecast(ourModel)
par(mfcol=c(3,4))
plot(ourModel, c(1:11))

# Model combination using a specified pool
ourModel <- adam(rnorm(100,100,10), model=c("ANN","AAN","MNN","CCC"),
                          lags=c(5,10))

# ADAM ARIMA
ourModel <- adam(rnorm(100,100,10), model="NNN",
                          lags=c(1,4), orders=list(ar=c(1,0),i=c(1,0),ma=c(1,1)))

# Fit ADAM to the data
ourModel <- adam(rnorm(100,100,10), model="AAdN")
# Simulate the data
x <- simulate(ourModel)

# Automatic selection of appropriate distribution and orders of ADAM ETS+ARIMA
ourModel <- auto.adam(rnorm(100,100,10), model="ZZN", lags=c(1,4),
                      orders=list(ar=c(2,2),ma=c(2,2),select=TRUE))

Complex Exponential Smoothing

Description

Function estimates CES in state space form with information potential equal to errors and returns several variables.

Usage

ces(y, seasonality = c("none", "simple", "partial", "full"),
  lags = c(frequency(y)), initial = c("backcasting", "optimal",
  "two-stage", "complete"), a = NULL, b = NULL, loss = c("likelihood",
  "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 0,
  holdout = FALSE, bounds = c("admissible", "none"), silent = TRUE,
  model = NULL, xreg = NULL, regressors = c("use", "select", "adapt"),
  initialX = NULL, ...)

auto.ces(y, seasonality = c("none", "simple", "partial", "full"),
  lags = c(frequency(y)), initial = c("backcasting", "optimal",
  "two-stage", "complete"), ic = c("AICc", "AIC", "BIC", "BICc"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE",
  "MSCE"), h = 0, holdout = FALSE, bounds = c("admissible", "none"),
  silent = TRUE, xreg = NULL, regressors = c("use", "select", "adapt"),
  ...)

ces_old(y, seasonality = c("none", "simple", "partial", "full"),
  initial = c("backcasting", "optimal"), a = NULL, b = NULL,
  ic = c("AICc", "AIC", "BIC", "BICc"), loss = c("likelihood", "MSE",
  "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10, holdout = FALSE,
  bounds = c("admissible", "none"), silent = c("all", "graph", "legend",
  "output", "none"), xreg = NULL, regressors = c("use", "select"),
  initialX = NULL, ...)

Arguments

Details

The function estimates Complex Exponential Smoothing in the state space form described in Svetunkov et al. (2022) with the information potential equal to the approximation error.

The auto.ces() function implements the automatic seasonal component selection based on information criteria.

ces_old() is the old implementation of the model and will be discontinued starting from smooth v4.5.0.

ces() uses two optimisers to get good estimates of parameters. By default these are BOBYQA and then Nelder-Mead. This can be regulated via ... - see details below.

For some more information about the model and its implementation, see the vignette: vignette("ces","smooth")

Value

Object of class "adam" is returned with similar elements to the adam function.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

• 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: https://openforecast.org/adam/.

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., Kourentzes, N., & Ord, J. K. (2022). Complex exponential smoothing. Naval Research Logistics, 69(8), 1108–1123. https://doi.org/10.1002/nav.22074

See Also

adam, es

Examples

y <- rnorm(100,10,3)
ces(y, h=20, holdout=FALSE)

y <- 500 - c(1:100)*0.5 + rnorm(100,10,3)
ces(y, h=20, holdout=TRUE)

ces(BJsales, h=8, holdout=TRUE)

ces(AirPassengers, h=18, holdout=TRUE, seasonality="s")
ces(AirPassengers, h=18, holdout=TRUE, seasonality="p")
ces(AirPassengers, h=18, holdout=TRUE, seasonality="f")

y <- ts(rnorm(100,10,3),frequency=12)
# CES with and without holdout
auto.ces(y,h=20,holdout=TRUE)
auto.ces(y,h=20,holdout=FALSE)


# Selection between "none" and "full" seasonalities
auto.ces(AirPassengers, h=12, holdout=TRUE,
                   seasonality=c("n","f"), ic="AIC")

ourModel <- auto.ces(AirPassengers)

summary(ourModel)
forecast(ourModel, h=12)

Centered Moving Average

Description

Function constructs centered moving average based on state space SMA

Usage

cma(y, order = NULL, silent = TRUE, ...)

Arguments

Details

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

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 monthly data and we use order=12, then half of the first January and half of the new one is taken.

This is not a forecasting tool. This is supposed to smooth the time series in order to find trend. So don't expect any forecasts from this function!

Value

Object of class "smooth" is returned. It contains the list of the following values:

• model - the name of the estimated model.

• timeElapsed - time elapsed for the construction of the model.

• order - order of the moving average.

• nParam - table with the number of estimated / provided parameters. If a previous model was reused, then its initials are reused and the number of provided parameters will take this into account.

• fitted - the fitted values, shifted in time.

• forecast - NAs, because this function does not produce forecasts.

• residuals - the residuals of the SMA / AR model.

• s2 - variance of the residuals (taking degrees of freedom into account) of the SMA / AR model.

• y - the original data.

• ICs - values of information criteria from the respective SMA or AR model. Includes AIC, AICc, BIC and BICc.

• logLik - log-likelihood of the SMA / AR model.

• lossValue - Cost function value (for the SMA / AR model).

• loss - Type of loss function used in the estimation.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

See Also

es, ssarima

Examples

# CMA of specific order
ourModel <- cma(rnorm(118,100,3),order=12)

# CMA of arbitrary order
ourModel <- cma(rnorm(118,100,3))

summary(ourModel)

Exponential Smoothing in SSOE state space model

Description

Function constructs ETS model and returns forecast, fitted values, errors and matrix of states. It is a wrapper of adam function.

Usage

es(y, model = "ZZZ", lags = c(frequency(y)), persistence = NULL,
  phi = NULL, initial = c("backcasting", "optimal", "two-stage",
  "complete"), initialSeason = NULL, ic = c("AICc", "AIC", "BIC", "BICc"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE",
  "MSCE"), h = 10, holdout = FALSE, bounds = c("usual", "admissible",
  "none"), silent = TRUE, xreg = NULL, regressors = c("use", "select"),
  initialX = NULL, ...)

es_old(y, model = "ZZZ", persistence = NULL, phi = NULL,
  initial = c("optimal", "backcasting"), initialSeason = NULL,
  ic = c("AICc", "AIC", "BIC", "BICc"), loss = c("likelihood", "MSE",
  "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10, holdout = FALSE,
  cumulative = FALSE, interval = c("none", "parametric", "likelihood",
  "semiparametric", "nonparametric"), level = 0.95, bounds = c("usual",
  "admissible", "none"), silent = c("all", "graph", "legend", "output",
  "none"), xreg = NULL, regressors = c("use", "select"), initialX = NULL,
  ...)

Arguments

Details

Function estimates ETS in a form of the Single Source of Error state space model of the following type:

y_{t} = o_t (w(v_{t-l}) + h(x_t, a_{t-1}) + r(v_{t-l}) \epsilon_{t})

v_{t} = f(v_{t-l}) + g(v_{t-l}) \epsilon_{t}

a_{t} = F_{X} a_{t-1} + g_{X} \epsilon_{t} / x_{t}

Where o_{t} is the Bernoulli distributed random variable (in case of normal data it equals to 1 for all observations), v_{t} is the state vector and l is the vector of lags, x_t is the vector of exogenous variables. w(.) is the measurement function, r(.) is the error function, f(.) is the transition function, g(.) is the persistence function and h(.) is the explanatory variables function. a_t is the vector of parameters for exogenous variables, F_{X} is the transitionX matrix and g_{X} is the persistenceX matrix. Finally, \epsilon_{t} is the error term.

For the details see Hyndman et al.(2008).

For some more information about the model and its implementation, see the vignette: vignette("es","smooth").

Also, there are posts about the functions of the package smooth on the website of Ivan Svetunkov: https://openforecast.org/category/r-en/smooth/ - they explain the underlying models and how to use the functions.

Value

Object of class "adam" is returned. It contains the list of the following values for classical ETS models:

• model - type of constructed model.

• formula - mathematical formula, describing interactions between components of es() and exogenous variables.

• timeElapsed - time elapsed for the construction of the model.

• states - matrix of the components of ETS.

• persistence - persistence vector. This is the place, where smoothing parameters live.

• phi - value of damping parameter.

• transition - transition matrix of the model.

• measurement - measurement vector of the model.

• initialType - type of the initial values used.

• initial - initial values of the state vector (non-seasonal).

• initialSeason - initial values of the seasonal part of state vector.

• nParam - table with the number of estimated / provided parameters. If a previous model was reused, then its initials are reused and the number of provided parameters will take this into account.

• fitted - fitted values of ETS. In case of the intermittent model, the fitted are multiplied by the probability of occurrence.

• forecast - the point forecast for h steps ahead (by default NA is returned). NOTE that these do not always correspond to the conditional expectations. See ADAM textbook, Section 4.4. for details (https://openforecast.org/adam/ETSTaxonomyMaths.html),

• lower - lower bound of prediction interval. When interval="none" then NA is returned.

• upper - higher bound of prediction interval. When interval="none" then NA is returned.

• residuals - residuals of the estimated model.

• errors - trace forecast in-sample errors, returned as a matrix. Only returned when the multistep losses are used and semiparametric interval is needed.

• s2 - variance of the residuals (taking degrees of freedom into account). This is an unbiased estimate of variance.

• interval - type of interval asked by user.

• level - confidence level for interval.

• cumulative - whether the produced forecast was cumulative or not.

• y - original data.

• holdout - holdout part of the original data.

• xreg - provided vector or matrix of exogenous variables. If regressors="s", then this value will contain only selected exogenous variables.

• initialX - initial values for parameters of exogenous variables.

• ICs - values of information criteria of the model. Includes AIC, AICc, BIC and BICc.

• logLik - concentrated log-likelihood of the function.

• lossValue - loss function value.

• loss - type of loss function used in the estimation.

• FI - Fisher Information. Equal to NULL if FI=FALSE or when FI is not provided at all.

• accuracy - vector of accuracy measures for the holdout sample. In case of non-intermittent data includes: MPE, MAPE, SMAPE, MASE, sMAE, RelMAE, sMSE and Bias coefficient (based on complex numbers). In case of intermittent data the set of errors will be: sMSE, sPIS, sCE (scaled cumulative error) and Bias coefficient. This is available only when holdout=TRUE.

• B - the vector of all the estimated parameters.

If combination of forecasts is produced (using model="CCC"), then a shorter list of values is returned:

• model,

• timeElapsed,

• initialType,

• fitted,

• forecast,

• lower,

• upper,

• residuals,

• s2 - variance of additive error of combined one-step-ahead forecasts,

• interval,

• level,

• cumulative,

• y,

• holdout,

• ICs - combined ic,

• ICw - ic weights used in the combination,

• loss,

• xreg,

• accuracy.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., Boylan, J.E., 2023a. iETS: State Space Model for Intermittent Demand Forecastings. International Journal of Production Economics. 109013. doi:10.1016/j.ijpe.2023.109013

• Teunter R., Syntetos A., Babai Z. (2011). Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214, 606-615.

• Croston, J. (1972) Forecasting and stock control for intermittent demands. Operational Research Quarterly, 23(3), 289-303.

• Kolassa, S. (2011) Combining exponential smoothing forecasts using Akaike weights. International Journal of Forecasting, 27, pp 238 - 251.

• Svetunkov, I., Boylan, J.E., 2023b. Staying Positive: Challenges and Solutions in Using Pure Multiplicative ETS Models. IMA Journal of Management Mathematics. p. 403-425. doi:10.1093/imaman/dpad028

• Taylor, J.W. and Bunn, D.W. (1999) A Quantile Regression Approach to Generating Prediction Intervals. Management Science, Vol 45, No 2, pp 225-237.

• Lichtendahl Kenneth C., Jr., Grushka-Cockayne Yael, Winkler Robert L., (2013) Is It Better to Average Probabilities or Quantiles? Management Science 59(7):1594-1611. DOI: doi:10.1287/mnsc.1120.1667

See Also

adam, forecast, ts, sim.es

Examples

# See how holdout and trace parameters influence the forecast
es(BJsales,model="AAdN",h=8,holdout=FALSE,loss="MSE")
es(AirPassengers,model="MAM",h=18,holdout=TRUE,loss="TMSE")

# Model selection example
es(BJsales,model="ZZN",ic="AIC",h=8,holdout=FALSE,bounds="a")

# Model selection. Compare AICc of these two models:
es(AirPassengers,"ZZZ",h=10,holdout=TRUE)
es(AirPassengers,"MAdM",h=10,holdout=TRUE)

# Model selection, excluding multiplicative trend
es(AirPassengers,model="ZXZ",h=8,holdout=TRUE)

# Combination example
es(BJsales,model="CCN",h=8,holdout=TRUE)

# Model selection using a specified pool of models
ourModel <- es(AirPassengers,model=c("ANN","AAM","AMdA"),h=18)

# Produce forecast and prediction interval
forecast(ourModel, h=18, interval="parametric")

# Semiparametric interval example
forecast(ourModel, h=18, interval="semiparametric")

# This will be the same model as in previous line but estimated on new portion of data
es(BJsales,model=ourModel,h=18,holdout=FALSE)

Forecasting time series using smooth functions

Description

Function produces conditional expectation (point forecasts) and prediction intervals for the estimated model.

Usage

## S3 method for class 'adam'
forecast(object, h = 10, newdata = NULL,
  occurrence = NULL, interval = c("none", "prediction", "confidence",
  "simulated", "approximate", "semiparametric", "nonparametric", "empirical",
  "complete"), level = 0.95, side = c("both", "upper", "lower"),
  cumulative = FALSE, nsim = NULL, scenarios = FALSE, ...)

## S3 method for class 'smooth'
forecast(object, h = 10, interval = c("parametric",
  "semiparametric", "nonparametric", "none"), level = 0.95,
  side = c("both", "upper", "lower"), ...)

## S3 method for class 'oes'
forecast(object, h = 10, ...)

## S3 method for class 'msdecompose'
forecast(object, h = 10, interval = c("parametric",
  "semiparametric", "nonparametric", "none"), level = 0.95, model = NULL,
  ...)

Arguments

Details

By default the function will generate conditional expectations from the estimated model and will also produce a variety of prediction intervals based on user preferences.

Value

Returns object of class "smooth.forecast", which contains:

• model - the estimated model (ES / CES / GUM / SSARIMA).

• method - the name of the estimated model (ES / CES / GUM / SSARIMA).

• forecast aka mean - point forecasts of the model (conditional mean).

• lower - lower bound of prediction interval.

• upper - upper bound of prediction interval.

• level - confidence level.

• interval - binary variable (whether interval were produced or not).

• scenarios - in case of forecast.adam() and interval="simulated" returns matrix with scenarios (future paths) that were used in simulations.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag.

See Also

forecast

Examples

ourModel <- es(rnorm(100,0,1), h=10)
forecast(ourModel, h=10, interval="parametric")

Generalised Univariate Model

Description

Function constructs Generalised Univariate Model, estimating matrices F, w, vector g and initial parameters.

Usage

gum(y, orders = c(1, 1), lags = c(1, frequency(y)), type = c("additive",
  "multiplicative"), initial = c("backcasting", "optimal", "two-stage",
  "complete"), persistence = NULL, transition = NULL,
  measurement = rep(1, sum(orders)), loss = c("likelihood", "MSE", "MAE",
  "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 0, holdout = FALSE,
  bounds = c("admissible", "none"), silent = TRUE, model = NULL,
  xreg = NULL, regressors = c("use", "select", "adapt", "integrate"),
  initialX = NULL, ...)

auto.gum(y, orders = 3, lags = frequency(y), type = c("additive",
  "multiplicative", "select"), initial = c("backcasting", "optimal",
  "two-stage", "complete"), ic = c("AICc", "AIC", "BIC", "BICc"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE",
  "MSCE"), h = 0, holdout = FALSE, bounds = c("admissible", "none"),
  silent = TRUE, xreg = NULL, regressors = c("use", "select", "adapt",
  "integrate"), ...)

gum_old(y, orders = c(1, 1), lags = c(1, frequency(y)),
  type = c("additive", "multiplicative"), persistence = NULL,
  transition = NULL, measurement = rep(1, sum(orders)),
  initial = c("optimal", "backcasting"), loss = c("likelihood", "MSE",
  "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10, holdout = FALSE,
  bounds = c("restricted", "admissible", "none"), silent = c("all",
  "graph", "legend", "output", "none"), xreg = NULL, regressors = c("use",
  "select"), initialX = NULL, ...)

ges(...)

Arguments

Details

The function estimates the Single Source of Error state space model of the following type:

y_{t} = w_t' v_{t-l} + \epsilon_{t}

v_{t} = F v_{t-l} + g_{t} \epsilon_{t}

where v_{t} is the state vector (defined using orders) and l is the vector of lags, w_t is the measurement vector (which includes fixed elements and explanatory variables), F is the transition matrix, g_t is the persistence vector (includes explanatory variables as well if provided), finally, \epsilon_{t} is the error term.

For some more information about the model and its implementation, see the vignette: vignette("gum","smooth")

Value

Object of class "adam" is returned with similar elements to the adam function.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

See Also

adam, es, ces

gum, es, ces, sim.es, ssarima

Examples

gum(BJsales, h=8, holdout=TRUE)

ourModel <- gum(rnorm(118,100,3), orders=c(2,1), lags=c(1,4), h=18, holdout=TRUE)

# Redo previous model on a new data and produce prediction interval
gum(rnorm(118,100,3), model=ourModel, h=18)

# Produce something crazy with optimal initials (not recommended)
gum(rnorm(118,100,3), orders=c(1,1,1), lags=c(1,3,5), h=18, holdout=TRUE, initial="o")

# Simpler model estimated using trace forecast error loss function and its analytical analogue
gum(rnorm(118,100,3), orders=c(1), lags=c(1), h=18, holdout=TRUE, bounds="n", loss="TMSE")


x <- rnorm(50,100,3)

# The best GUM model for the data
ourModel <- auto.gum(x, orders=2, lags=4, h=18, holdout=TRUE)

summary(ourModel)

Smooth classes checkers

Description

Functions to check if an object is of the specified class

Functions to check if an object is of the specified class

Usage

is.smooth(x)

is.smoothC(x)

is.msarima(x)

is.oes(x)

is.oesg(x)

is.smooth.sim(x)

is.smooth.forecast(x)

is.adam(x)

is.adam.sim(x)

is.msdecompose(x)

is.msdecompose.forecast(x)

Arguments

Details

The list of functions includes:

• is.smooth() tests if the object was produced by a smooth function (e.g. es / ces / ssarima / gum / sma / msarima);

• is.msarima() tests if the object was produced by the msarima function;

• is.smoothC() tests if the object was produced by a combination function (currently applies only to smoothCombine);

• is.oes() tests if the object was produced by oes function;

• is.smooth.sim() tests if the object was produced by simulate functions (e.g. sim.es / sim.ces / sim.ssarima / sim.sma / sim.gum);

• is.smooth.forecast() checks if the forecast was produced from a smooth function using forecast() function.

The list of functions includes:

• is.adam() tests if the object was produced by a adam function

• is.adam.sim() tests if the object was produced by sim.adam() function (not implemented yet);

Value

TRUE if this is the specified class and FALSE otherwise.

TRUE if this is the specified class and FALSE otherwise.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

Examples

ourModel <- msarima(rnorm(100,100,10))

is.smooth(ourModel)
is.msarima(ourModel)


ourModel <- adam(rnorm(100,100,10))
is.adam(ourModel)

Multiple Seasonal ARIMA

Description

Function constructs Multiple Seasonal State Space ARIMA, estimating AR, MA terms and initial states. It is a wrapper of adam function.

Usage

msarima(y, orders = list(ar = c(0), i = c(1), ma = c(1)), lags = c(1),
  constant = FALSE, AR = NULL, MA = NULL, model = NULL,
  initial = c("backcasting", "optimal", "two-stage", "complete"),
  ic = c("AICc", "AIC", "BIC", "BICc"), loss = c("likelihood", "MSE",
  "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10, holdout = FALSE,
  bounds = c("usual", "admissible", "none"), silent = TRUE, xreg = NULL,
  regressors = c("use", "select", "adapt"), initialX = NULL, ...)

auto.msarima(y, orders = list(ar = c(3, 3), i = c(2, 1), ma = c(3, 3)),
  lags = c(1, frequency(y)), initial = c("backcasting", "optimal",
  "two-stage", "complete"), ic = c("AICc", "AIC", "BIC", "BICc"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE",
  "MSCE"), h = 10, holdout = FALSE, bounds = c("usual", "admissible",
  "none"), silent = TRUE, xreg = NULL, regressors = c("use", "select",
  "adapt"), initialX = NULL, ...)

msarima_old(y, orders = list(ar = c(0), i = c(1), ma = c(1)), lags = c(1),
  constant = FALSE, AR = NULL, MA = NULL, initial = c("backcasting",
  "optimal"), ic = c("AICc", "AIC", "BIC", "BICc"), loss = c("likelihood",
  "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10,
  holdout = FALSE, cumulative = FALSE, interval = c("none", "parametric",
  "likelihood", "semiparametric", "nonparametric"), level = 0.95,
  bounds = c("admissible", "none"), silent = c("all", "graph", "legend",
  "output", "none"), xreg = NULL, regressors = c("use", "select"),
  initialX = NULL, ...)

Arguments

Details

The model, implemented in this function differs from the one in ssarima function (Svetunkov & Boylan, 2019), but it is more efficient and better fitting the data (which might be a limitation).

The basic ARIMA(p,d,q) used in the function has the following form:

(1 - B)^d (1 - a_1 B - a_2 B^2 - ... - a_p B^p) y_[t] = (1 + b_1 B + b_2 B^2 + ... + b_q B^q) \epsilon_[t] + c

where y_[t] is the actual values, \epsilon_[t] is the error term, a_i, b_j are the parameters for AR and MA respectively and c is the constant. In case of non-zero differences c acts as drift.

This model is then transformed into ARIMA in the Single Source of Error State space form (based by Snyder, 1985, but in a slightly different formulation):

y_{t} = o_{t} (w' v_{t-l} + x_t a_{t-1} + \epsilon_{t})

v_{t} = F v_{t-l} + g \epsilon_{t}

a_{t} = F_{X} a_{t-1} + g_{X} \epsilon_{t} / x_{t}

Where o_{t} is the Bernoulli distributed random variable (in case of normal data equal to 1), v_{t} is the state vector (defined based on orders) and l is the vector of lags, x_t is the vector of exogenous parameters. w is the measurement vector, F is the transition matrix, g is the persistence vector, a_t is the vector of parameters for exogenous variables, F_{X} is the transitionX matrix and g_{X} is the persistenceX matrix. The main difference from ssarima function is that this implementation skips zero polynomials, substantially decreasing the dimension of the transition matrix. As a result, this function works faster than ssarima on high frequency data, and it is more accurate.

Due to the flexibility of the model, multiple seasonalities can be used. For example, something crazy like this can be constructed: SARIMA(1,1,1)(0,1,1)[24](2,0,1)[24*7](0,0,1)[24*30], but the estimation may take some time... Still this should be estimated in finite time (not like with ssarima).

The auto.msarima function constructs several ARIMA models in Single Source of Error state space form based on adam function (see adam documentation) and selects the best one based on the selected information criterion.

For some additional details see the vignettes: vignette("adam","smooth") and vignette("ssarima","smooth")

Value

Object of class "adam" is returned. It contains the list of the following values:

• model - the name of the estimated model.

• timeElapsed - time elapsed for the construction of the model.

• states - the matrix of the fuzzy components of msarima, where rows correspond to time and cols to states.

• transition - matrix F.

• persistence - the persistence vector. This is the place, where smoothing parameters live.

• measurement - measurement vector of the model.

• AR - the matrix of coefficients of AR terms.

• I - the matrix of coefficients of I terms.

• MA - the matrix of coefficients of MA terms.

• constant - the value of the constant term.

• initialType - Type of the initial values used.

• initial - the initial values of the state vector (extracted from states).

• nParam - table with the number of estimated / provided parameters. If a previous model was reused, then its initials are reused and the number of provided parameters will take this into account.

• fitted - the fitted values.

• forecast - the point forecast.

• lower - the lower bound of prediction interval. When interval="none" then NA is returned.

• upper - the higher bound of prediction interval. When interval="none" then NA is returned.

• residuals - the residuals of the estimated model.

• errors - The matrix of 1 to h steps ahead errors. Only returned when the multistep losses are used and semiparametric interval is needed.

• s2 - variance of the residuals (taking degrees of freedom into account).

• interval - type of interval asked by user.

• level - confidence level for interval.

• cumulative - whether the produced forecast was cumulative or not.

• y - the original data.

• holdout - the holdout part of the original data.

• xreg - provided vector or matrix of exogenous variables. If regressors="s", then this value will contain only selected exogenous variables.

• initialX - initial values for parameters of exogenous variables.

• ICs - values of information criteria of the model. Includes AIC, AICc, BIC and BICc.

• logLik - log-likelihood of the function.

• lossValue - Cost function value.

• loss - Type of loss function used in the estimation.

• FI - Fisher Information. Equal to NULL if FI=FALSE or when FI is not provided at all.

• accuracy - vector of accuracy measures for the holdout sample. In case of non-intermittent data includes: MPE, MAPE, SMAPE, MASE, sMAE, RelMAE, sMSE and Bias coefficient (based on complex numbers). In case of intermittent data the set of errors will be: sMSE, sPIS, sCE (scaled cumulative error) and Bias coefficient. This is available only when holdout=TRUE.

• B - the vector of all the estimated parameters.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., & Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. International Journal of Production Research, 0(0), 1–10. doi:10.1080/00207543.2019.1600764

See Also

adam, orders, es, auto.ssarima

Examples

x <- rnorm(118,100,3)

# The simple call of ARIMA(1,1,1):
ourModel <- msarima(x,orders=c(1,1,1),lags=1,h=18,holdout=TRUE)

# Example of SARIMA(2,0,0)(1,0,0)[4]
msarima(x,orders=list(ar=c(2,1)),lags=c(1,4),h=18,holdout=TRUE)

# SARIMA of a peculiar order on AirPassengers data
ourModel <- msarima(AirPassengers,orders=list(ar=c(1,0,3),i=c(1,0,1),ma=c(0,1,2)),
                    lags=c(1,6,12),h=10,holdout=TRUE)

# ARIMA(1,1,1) with Mean Squared Trace Forecast Error
msarima(x,orders=c(1,1,1),lags=1,h=18,holdout=TRUE,loss="TMSE")

plot(forecast(ourModel, h=18, interval="approximate"))


x <- rnorm(118,100,3)
# The best ARIMA for the data
ourModel <- auto.msarima(x,orders=list(ar=c(2,1),i=c(1,1),ma=c(2,1)),lags=c(1,12),
                     h=18,holdout=TRUE)

# The other one using optimised states
auto.msarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
                       h=18,holdout=TRUE)

# And now combined ARIMA
auto.msarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
                       combine=TRUE,h=18,holdout=TRUE)

plot(forecast(ourModel, h=18, interval="simulated"))

Multiple seasonal classical decomposition

Description

Function decomposes multiple seasonal time series into components using the principles of classical decomposition.

Usage

msdecompose(y, lags = c(12), type = c("additive", "multiplicative"),
  smoother = c("ma", "lowess", "supsmu"), ...)

Arguments

Details

The function applies centred moving averages based on filter function and order specified in lags variable in order to smooth the original series and obtain level, trend and seasonal components of the series.

Value

The object of the class "msdecompose" is return, containing:

• y - the original time series.

• initial - the estimates of the initial level and trend.

• trend - the long term trend in the data.

• seasonal - the list of seasonal parameters.

• lags - the provided lags.

• type - the selected type of the decomposition.

• yName - the name of the provided data.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

See Also

filter

Examples

# Decomposition of multiple frequency data
## Not run: ourModel <- msdecompose(forecast::taylor, lags=c(48,336), type="m")
ourModel <- msdecompose(AirPassengers, lags=c(12), type="m")

plot(ourModel)
plot(forecast(ourModel, model="AAN", h=12))

Function returns the multiple steps ahead covariance matrix of forecast errors

Description

This function extracts covariance matrix of 1 to h steps ahead forecast errors for adam(), ssarima(), gum(), sma(), es() and ces() models.

Usage

multicov(object, type = c("analytical", "empirical", "simulated"), h = 10,
  nsim = 1000, ...)

## S3 method for class 'smooth'
multicov(object, type = c("analytical", "empirical",
  "simulated"), h = 10, nsim = 1000, ...)

Arguments

Details

The function returns either scalar (if it is a non-smooth model) or the matrix of (h x h) size with variances and covariances of 1 to h steps ahead forecast errors.

Value

Scalar in cases of non-smooth functions. (h x h) matrix otherwise.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

See Also

orders

Examples

x <- rnorm(100,0,1)

# A simple example with a 5x5 covariance matrix
ourModel <- ces(x, h=5)
multicov(ourModel)

Occurrence ETS model

Description

Function returns the occurrence part of iETS model with the specified probability update and model types.

Usage

oes(y, model = "MNN", persistence = NULL, initial = "o",
  initialSeason = NULL, phi = NULL, occurrence = c("fixed", "general",
  "odds-ratio", "inverse-odds-ratio", "direct", "auto", "none"),
  ic = c("AICc", "AIC", "BIC", "BICc"), h = 10, holdout = FALSE,
  bounds = c("usual", "admissible", "none"), silent = c("all", "graph",
  "legend", "output", "none"), xreg = NULL, regressors = c("use",
  "select"), initialX = NULL, ...)

Arguments

Details

The function estimates probability of demand occurrence, using the selected ETS state space models.

For the details about the model and its implementation, see the respective vignette: vignette("oes","smooth")

Value

The object of class "occurrence" is returned. It contains following list of values:

• model - the type of the estimated ETS model;

• timeElapsed - the time elapsed for the construction of the model;

• fitted - the fitted values for the probability;

• fittedModel - the fitted values of the underlying ETS model, where applicable (only for occurrence=c("o","i","d"));

• forecast - the forecast of the probability for h observations ahead;

• forecastModel - the forecast of the underlying ETS model, where applicable (only for occurrence=c("o","i","d"));

• lower - the lower bound of the interval if interval!="none";

• upper - the upper bound of the interval if interval!="none";

• lowerModel - the lower bound of the interval of the underlying ETS model if interval!="none";

• upperModel - the upper bound of the interval of the underlying ETS model if interval!="none";

• states - the values of the state vector;

• logLik - the log-likelihood value of the model;

• nParam - the number of parameters in the model (the matrix is returned);

• residuals - the residuals of the model;

• y - actual values of occurrence (zeros and ones).

• persistence - the vector of smoothing parameters;

• phi - the value of the damped trend parameter;

• initial - initial values of the state vector;

• initialSeason - the matrix of initials seasonal states;

• occurrence - the type of the occurrence model;

• updateX - boolean, defining, if the states of exogenous variables were estimated as well.

• initialX - initial values for parameters of exogenous variables.

• persistenceX - persistence vector g for exogenous variables.

• transitionX - transition matrix F for exogenous variables.

• accuracy - The error measures for the forecast (in case of holdout=TRUE).

• B - the vector of all the estimated parameters (in case of "odds-ratio", "inverse-odds-ratio" and "direct" models).

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., Boylan, J.E., 2023a. iETS: State Space Model for Intermittent Demand Forecastings. International Journal of Production Economics. 109013. doi:10.1016/j.ijpe.2023.109013

• Teunter R., Syntetos A., Babai Z. (2011). Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214, 606-615.

• Croston, J. (1972) Forecasting and stock control for intermittent demands. Operational Research Quarterly, 23(3), 289-303.

See Also

adam, oesg, es

Examples

y <- rpois(100,0.1)
oes(y, occurrence="auto")

oes(y, occurrence="f")

Occurrence ETS, general model

Description

Function returns the general occurrence model of the of iETS model.

Usage

oesg(y, modelA = "MNN", modelB = "MNN", persistenceA = NULL,
  persistenceB = NULL, phiA = NULL, phiB = NULL, initialA = "o",
  initialB = "o", initialSeasonA = NULL, initialSeasonB = NULL,
  ic = c("AICc", "AIC", "BIC", "BICc"), h = 10, holdout = FALSE,
  bounds = c("usual", "admissible", "none"), silent = c("all", "graph",
  "legend", "output", "none"), xregA = NULL, xregB = NULL,
  initialXA = NULL, initialXB = NULL, regressorsA = c("use", "select"),
  regressorsB = c("use", "select"), ...)

Arguments

Details

The function estimates probability of demand occurrence, based on the iETS_G state-space model. It involves the estimation and modelling of the two simultaneous state space equations. Thus two parts for the model type, persistence, initials etc.

For the details about the model and its implementation, see the respective vignette: vignette("oes","smooth")

The model is based on:

o_t \sim Bernoulli(p_t)

p_t = \frac{a_t}{a_t+b_t}

,

where a_t and b_t are the parameters of the Beta distribution and are modelled using separate ETS models.

Value

The object of class "occurrence" is returned. It contains following list of values:

• modelA - the model A of the class oes, that contains the output similar to the one from the oes() function;

• modelB - the model B of the class oes, that contains the output similar to the one from the oes() function.

• B - the vector of all the estimated parameters.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

See Also

es, oes

Examples

y <- rpois(100,0.1)
oesg(y, modelA="MNN", modelB="ANN")

Functions that extract values from the fitted model

Description

These functions allow extracting orders and lags for ssarima(), gum() and sma(), type of model from es() and ces() and name of model.

Usage

orders(object, ...)

lags(object, ...)

modelName(object, ...)

modelType(object, ...)

Arguments

Details

orders() and lags() are useful only for SSARIMA, GUM and SMA. They return NA for other functions. This can also be applied to arima(), Arima() and auto.arima() functions from stats and forecast packages. modelType() is useful only for ETS and CES. They return NA for other functions. This can also be applied to ets() function from forecast package. errorType extracts the type of error from the model (either additive or multiplicative). Finally, modelName returns the name of the fitted model. For example, "ARIMA(0,1,1)". This is purely descriptive and can also be applied to non-smooth classes, so that, for example, you can easily extract the name of the fitted AR model from ar() function from stats package.

Value

Either vector, scalar or list with values is returned. orders() in case of ssarima returns list of values:

• ar - AR orders.

• i - I orders.

• ma - MA orders.

lags() returns the vector of lags of the model. All the other functions return strings of character.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

See Also

ssarima, msarima

Examples

x <- rnorm(100,0,1)

# Just as example. orders and lags do not return anything for ces() and es(). But modelType() does.
ourModel <- ces(x, h=10)
orders(ourModel)
lags(ourModel)
modelType(ourModel)
modelName(ourModel)

# And as another example it does the opposite for gum() and ssarima()
ourModel <- gum(x, h=10, orders=c(1,1), lags=c(1,4))
orders(ourModel)
lags(ourModel)
modelType(ourModel)
modelName(ourModel)

# Finally these values can be used for simulate functions or original functions.
ourModel <- auto.ssarima(x)
ssarima(x, orders=orders(ourModel), lags=lags(ourModel), constant=ourModel$constant)
sim.ssarima(orders=orders(ourModel), lags=lags(ourModel), constant=ourModel$constant)

ourModel <- es(x)
es(x, model=modelType(ourModel))
sim.es(model=modelType(ourModel))

Plots for the fit and states

Description

The function produces diagnostics plots for a smooth model

Usage

## S3 method for class 'adam'
plot(x, which = c(1, 2, 4, 6), level = 0.95,
  legend = FALSE, ask = prod(par("mfcol")) < length(which) &&
  dev.interactive(), lowess = TRUE, ...)

## S3 method for class 'smooth'
plot(x, which = c(1, 2, 4, 6), level = 0.95,
  legend = FALSE, ask = prod(par("mfcol")) < length(which) &&
  dev.interactive(), lowess = TRUE, ...)

## S3 method for class 'msdecompose'
plot(x, which = c(1, 2, 4, 6), level = 0.95,
  legend = FALSE, ask = prod(par("mfcol")) < length(which) &&
  dev.interactive(), lowess = TRUE, ...)

Arguments

Details

The list of produced plots includes:

1. Actuals vs Fitted values. Allows analysing, whether there are any issues in the fit. Does the variability of actuals increase with the increase of fitted values? Is the relation well captured? They grey line on the plot corresponds to the perfect fit of the model.

2. Standardised residuals vs Fitted. Plots the points and the confidence bounds (red lines) for the specified confidence level. Useful for the analysis of outliers;

3. Studentised residuals vs Fitted. This is similar to the previous plot, but with the residuals divided by the scales with the leave-one-out approach. Should be more sensitive to outliers;

4. Absolute residuals vs Fitted. Useful for the analysis of heteroscedasticity;

5. Squared residuals vs Fitted - similar to (3), but with squared values;

6. Q-Q plot with the specified distribution. Can be used in order to see if the residuals follow the assumed distribution. The type of distribution depends on the one used in the estimation (see distribution parameter in alm);

7. ACF of the residuals. Are the residuals autocorrelated? See acf for details;

8. Fitted over time. Plots actuals (black line), fitted values (purple line), point forecast (blue line) and prediction interval (grey lines). Can be used in order to make sure that the model did not miss any important events over time;

9. Standardised residuals vs Time. Useful if you want to see, if there is autocorrelation or if there is heteroscedasticity in time. This also shows, when the outliers happen;

10. Studentised residuals vs Time. Similar to previous, but with studentised residuals;

11. PACF of the residuals. No, really, are they autocorrelated? See pacf function from stats package for details;

12. Plot of the states of the model. It is not recommended to produce this plot together with the others, because there might be several states, which would cause the creation of a different canvas. In case of "msdecompose", this will produce the decomposition of the series into states on a different canvas;

13. Absolute standardised residuals vs Fitted. Similar to the previous, but with absolute values. This is more relevant to the models where scale is calculated as an absolute value of something (e.g. Laplace);

14. Squared standardised residuals vs Fitted. This is an additional plot needed to diagnose heteroscedasticity in a model with varying scale. The variance on this plot will be constant if the adequate model for scale was constructed. This is more appropriate for normal and the related distributions.

Which of the plots to produce, is specified via the which parameter.

Value

The function produces the number of plots, specified in the parameter which.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

See Also

plot.greybox

Examples

ourModel <- es(c(rnorm(50,100,10),rnorm(50,120,10)), "ANN", h=10)
par(mfcol=c(3,4))
plot(ourModel, c(1:11))
plot(ourModel, 12)

Prediction Likelihood Score

Description

Function estimates Prediction Likelihood Score for the provided model

Usage

pls(object, holdout = NULL, ...)

## S3 method for class 'smooth'
pls(object, holdout = NULL, ...)

Arguments

Details

Prediction likelihood score (PLS) is based on either normal or log-normal distribution of errors. This is extracted from the provided model. The likelihood based on the distribution of 1 to h steps ahead forecast errors is used in the process.

Value

A value of the log-likelihood.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

distribution. IEEE Signal Processing Letters. 13 (5): 300-303. doi:10.1109/LSP.2006.870353 - this is not yet used in the function.

• Snyder, R. D., Ord, J. K., Beaumont, A., 2012. Forecasting the intermittent demand for slow-moving inventories: A modelling approach. International Journal of Forecasting 28 (2), 485-496.

• Kolassa, S., 2016. Evaluating predictive count data distributions in retail sales forecasting. International Journal of Forecasting 32 (3), 788-803..

Examples

# Generate data, apply es() with the holdout parameter and calculate PLS
x <- rnorm(100,0,1)
ourModel <- es(x, h=10, holdout=TRUE)
pls(ourModel, type="a")
pls(ourModel, type="e")
pls(ourModel, type="s", obs=100, nsim=100)

Reapply the model with randomly generated initial parameters and produce forecasts

Description

reapply function generates the parameters based on the values in the provided object and then reapplies the same model with those parameters to the data, getting the fitted paths and updated states. reforecast function uses those values in order to produce forecasts for the h steps ahead.

Usage

reapply(object, nsim = 1000, bootstrap = FALSE, heuristics = NULL, ...)

reforecast(object, h = 10, newdata = NULL, occurrence = NULL,
  interval = c("prediction", "confidence", "none"), level = 0.95,
  side = c("both", "upper", "lower"), cumulative = FALSE, nsim = 100,
  ...)

Arguments

Details

The main motivation of the function is to take the randomness due to the in-sample estimation of parameters into account when fitting the model and to propagate this randomness to the forecasts. The methods can be considered as a special case of recursive bootstrap.

Value

reapply() returns object of the class "reapply", which contains:

• timeElapsed - Time elapsed for the code execution;

• y - The actual values;

• states - The array of states of the model;

• refitted - The matrix with fitted values, where columns correspond to different paths;

• fitted - The vector of fitted values (conditional mean);

• model - The name of the constructed model;

• transition - The array of transition matrices;

• measurement - The array of measurement matrices;

• persistence - The matrix of persistence vectors (paths in columns);

• profile - The array of profiles obtained by the end of each fit.

reforecast() returns the object of the class forecast.smooth, which contains in addition to the standard list the variable paths - all simulated trajectories with h in rows, simulated future paths for each state in columns and different states (obtained from reapply() function) in the third dimension.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

See Also

forecast.smooth

Examples

x <- rnorm(100,0,1)

# Just as example. orders and lags do not return anything for ces() and es(). But modelType() does.
ourModel <- adam(x, "ANN")
refittedModel <- reapply(ourModel, nsim=50)
plot(refittedModel)

ourForecast <- reforecast(ourModel, nsim=50)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

accuracy

xtable

xtable

Multiple steps ahead forecast errors

Description

The function extracts 1 to h steps ahead forecast errors from the model.

Usage

rmultistep(object, h = 10, ...)

Arguments

Details

The errors correspond to the error term epsilon_t in the ETS models. Don't forget that different models make different assumptions about epsilon_t and / or 1+epsilon_t.

Value

The matrix with observations in rows and h steps ahead values in columns. So, the first row corresponds to the forecast produced from the 0th observation from 1 to h steps ahead.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

See Also

residuals,

Examples

x <- rnorm(100,0,1)
ourModel <- adam(x)
rmultistep(ourModel, h=13)

Simulate Complex Exponential Smoothing

Description

Function generates data using CES with Single Source of Error as a data generating process.

Usage

sim.ces(seasonality = c("none", "simple", "partial", "full"), obs = 10,
  nsim = 1, frequency = 1, a = NULL, b = NULL, initial = NULL,
  randomizer = c("rnorm", "rt", "rlaplace", "rs"), probability = 1, ...)

Arguments

Details

For the information about the function, see the vignette: vignette("simulate","smooth")

Value

List of the following values is returned:

• model - Name of CES model.

• a - Value of complex smoothing parameter a. If nsim>1, then this is a vector.

• b - Value of complex smoothing parameter b. If seasonality="n" or seasonality="s", then this is equal to NULL. If nsim>1, then this is a vector.

• initial - Initial values of CES in a form of matrix. If nsim>1, then this is an array.

• data - Time series vector (or matrix if nsim>1) of the generated series.

• states - Matrix (or array if nsim>1) of states. States are in columns, time is in rows.

• residuals - Error terms used in the simulation. Either vector or matrix, depending on nsim.

• occurrence - Values of occurrence variable. Once again, can be either a vector or a matrix...

• logLik - Log-likelihood of the constructed model.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., Kourentzes, N., & Ord, J. K. (2022). Complex exponential smoothing. Naval Research Logistics, 69(8), 1108–1123. https://doi.org/10.1002/nav.22074

See Also

sim.es, sim.ssarima, ces, Distributions

Examples

# Create 120 observations from CES(n). Generate 100 time series of this kind.
x <- sim.ces("n",obs=120,nsim=100)

# Generate similar thing for seasonal series of CES(s)_4
x <- sim.ces("s",frequency=4,obs=80,nsim=100)

# Estimate model and then generate 10 time series from it
ourModel <- ces(rnorm(100,100,5))
simulate(ourModel,nsim=10)

Simulate Exponential Smoothing

Description

Function generates data using ETS with Single Source of Error as a data generating process.

Usage

sim.es(model = "ANN", obs = 10, nsim = 1, frequency = 1,
  persistence = NULL, phi = 1, initial = NULL, initialSeason = NULL,
  bounds = c("usual", "admissible", "restricted"), randomizer = c("rnorm",
  "rlnorm", "rt", "rlaplace", "rs"), probability = 1, ...)

Arguments

Details

For the information about the function, see the vignette: vignette("simulate","smooth")

Value

List of the following values is returned:

• model - Name of ETS model.

• data - Time series vector (or matrix if nsim>1) of the generated series.

• states - Matrix (or array if nsim>1) of states. States are in columns, time is in rows.

• persistence - Vector (or matrix if nsim>1) of smoothing parameters used in the simulation.

• phi - Value of damping parameter used in time series generation.

• initial - Vector (or matrix) of initial values.

• initialSeason - Vector (or matrix) of initial seasonal coefficients.

• probability - vector of probabilities used in the simulation.

• intermittent - type of the intermittent model used.

• residuals - Error terms used in the simulation. Either vector or matrix, depending on nsim.

• occurrence - Values of occurrence variable. Once again, can be either a vector or a matrix...

• logLik - Log-likelihood of the constructed model.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

See Also

es, ts, Distributions

Examples

# Create 40 observations of quarterly data using AAA model with errors from normal distribution
ETSAAA <- sim.es(model="AAA",frequency=4,obs=40,randomizer="rnorm",mean=0,sd=100)

# Create 50 series of quarterly data using AAA model
# with 40 observations each with errors from normal distribution
ETSAAA <- sim.es(model="AAA",frequency=4,obs=40,randomizer="rnorm",mean=0,sd=100,nsim=50)

# Create 50 series of quarterly data using AAdA model
# with 40 observations each with errors from normal distribution
# and smoothing parameters lying in the "admissible" range.
ETSAAA <- sim.es(model="AAA",phi=0.9,frequency=4,obs=40,bounds="admissible",
                  randomizer="rnorm",mean=0,sd=100,nsim=50)

# Create 60 observations of monthly data using ANN model
# with errors from beta distribution
ETSANN <- sim.es(model="ANN",persistence=c(1.5),frequency=12,obs=60,
                  randomizer="rbeta",shape1=1.5,shape2=1.5)
plot(ETSANN$states)

# Create 60 observations of monthly data using MAM model
# with errors from uniform distribution
ETSMAM <- sim.es(model="MAM",persistence=c(0.3,0.2,0.1),initial=c(2000,50),
           phi=0.8,frequency=12,obs=60,randomizer="runif",min=-0.5,max=0.5)

# Create 80 observations of quarterly data using MMM model
# with predefined initial values and errors from the normal distribution
ETSMMM <- sim.es(model="MMM",persistence=c(0.1,0.1,0.1),initial=c(2000,1),
           initialSeason=c(1.1,1.05,0.9,.95),frequency=4,obs=80,mean=0,sd=0.01)

# Generate intermittent data using AAdN
iETSAAdN <- sim.es("AAdN",obs=30,frequency=1,probability=0.1,initial=c(3,0),phi=0.8)

# Generate iETS(MNN) with TSB style probabilities
oETSMNN <- sim.oes("MNN",obs=50,occurrence="d",persistence=0.2,initial=1,
                   randomizer="rlnorm",meanlog=0,sdlog=0.3)
iETSMNN <- sim.es("MNN",obs=50,frequency=12,persistence=0.2,initial=4,
                  probability=oETSMNN$probability)

Simulate Generalised Exponential Smoothing

Description

Function generates data using GUM with Single Source of Error as a data generating process.

Usage

sim.gum(orders = c(1), lags = c(1), obs = 10, nsim = 1,
  frequency = 1, measurement = NULL, transition = NULL,
  persistence = NULL, initial = NULL, randomizer = c("rnorm", "rt",
  "rlaplace", "rs"), probability = 1, ...)

Arguments

Details

For the information about the function, see the vignette: vignette("simulate","smooth")

Value

List of the following values is returned:

• model - Name of GUM model.

• measurement - Matrix w.

• transition - Matrix F.

• persistence - Persistence vector. This is the place, where smoothing parameters live.

• initial - Initial values of GUM in a form of matrix. If nsim>1, then this is an array.

• data - Time series vector (or matrix if nsim>1) of the generated series.

• states - Matrix (or array if nsim>1) of states. States are in columns, time is in rows.

• residuals - Error terms used in the simulation. Either vector or matrix, depending on nsim.

• occurrence - Values of occurrence variable. Once again, can be either a vector or a matrix...

• logLik - Log-likelihood of the constructed model.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

See Also

sim.es, sim.ssarima, sim.ces, gum, Distributions

Examples

# Create 120 observations from GUM(1[1]). Generate 100 time series of this kind.
x <- sim.gum(orders=c(1),lags=c(1),obs=120,nsim=100)

# Generate similar thing for seasonal series of GUM(1[1],1[4]])
x <- sim.gum(orders=c(1,1),lags=c(1,4),frequency=4,obs=80,nsim=100,transition=c(1,0,0.9,0.9))

# Estimate model and then generate 10 time series from it
ourModel <- gum(rnorm(100,100,5))
simulate(ourModel,nsim=10)

Simulate Occurrence Part of ETS model

Description

Function generates data using ETS with Single Source of Error as a data generating process for the demand occurrence. As the main output it produces probabilities of occurrence.

Usage

sim.oes(model = "MNN", obs = 10, nsim = 1, frequency = 1,
  occurrence = c("odds-ratio", "inverse-odds-ratio", "direct", "general"),
  bounds = c("usual", "admissible", "restricted"), randomizer = c("rlnorm",
  "rinvgauss", "rgamma", "rnorm"), persistence = NULL, phi = 1,
  initial = NULL, initialSeason = NULL, modelB = model,
  persistenceB = persistence, phiB = phi, initialB = initial,
  initialSeasonB = initialSeason, ...)

Arguments

Details

For the information about the function, see the vignette: vignette("simulate","smooth")

Value

List of the following values is returned:

• model - Name of ETS model.

• modelA - Model A, generated using sim.es() function;

• modelB - Model B, generated using sim.es() function;

• probability - The value of probability generated by the model;

• occurrence - Type of occurrence used in the model;

• logLik - Log-likelihood of the constructed model.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

See Also

oes, sim.es, Distributions

Examples

# This example uses rinvgauss function from statmod package.
oETSMNNIG <- sim.oes(model="MNN",frequency=12,obs=60,
                               randomizer="rinvgauss",mean=1,dispersion=0.5)

# A simpler example with log normal distribution
oETSMNNlogN <- sim.oes(model="MNN",frequency=12,obs=60,initial=1,
                       randomizer="rlnorm",meanlog=0,sdlog=0.1)

Simulate Simple Moving Average

Description

Function generates data using SMA in a Single Source of Error state space model as a data generating process.

Usage

sim.sma(order = NULL, obs = 10, nsim = 1, frequency = 1,
  initial = NULL, randomizer = c("rnorm", "rt", "rlaplace", "rs"),
  probability = 1, ...)

Arguments

Details

For the information about the function, see the vignette: vignette("simulate","smooth")

Value

List of the following values is returned:

• model - Name of SMA model.

• data - Time series vector (or matrix if nsim>1) of the generated series.

• states - Matrix (or array if nsim>1) of states. States are in columns, time is in rows.

• initial - Vector (or matrix) of initial values.

• probability - vector of probabilities used in the simulation.

• intermittent - type of the intermittent model used.

• residuals - Error terms used in the simulation. Either vector or matrix, depending on nsim.

• occurrence - Values of occurrence variable. Once again, can be either a vector or a matrix...

• logLik - Log-likelihood of the constructed model.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

See Also

es, ts, Distributions

Examples

# Create 40 observations of quarterly data using AAA model with errors from normal distribution
sma10 <- sim.sma(order=10,frequency=4,obs=40,randomizer="rnorm",mean=0,sd=100)

Simulate SSARIMA

Description

Function generates data using SSARIMA with Single Source of Error as a data generating process.

Usage

sim.ssarima(orders = list(ar = 0, i = 1, ma = 1), lags = 1, obs = 10,
  nsim = 1, frequency = 1, AR = NULL, MA = NULL, constant = FALSE,
  initial = NULL, bounds = c("admissible", "none"),
  randomizer = c("rnorm", "rt", "rlaplace", "rs"), probability = 1, ...)

Arguments

Details

For the information about the function, see the vignette: vignette("simulate","smooth")

Value

List of the following values is returned:

• model - Name of SSARIMA model.

• AR - Value of AR parameters. If nsim>1, then this is a matrix.

• MA - Value of MA parameters. If nsim>1, then this is a matrix.

• constant - Value of constant term. If nsim>1, then this is a vector.

• initial - Initial values of SSARIMA. If nsim>1, then this is a matrix.

• data - Time series vector (or matrix if nsim>1) of the generated series.

• states - Matrix (or array if nsim>1) of states. States are in columns, time is in rows.

• residuals - Error terms used in the simulation. Either vector or matrix, depending on nsim.

• occurrence - Values of occurrence variable. Once again, can be either a vector or a matrix...

• logLik - Log-likelihood of the constructed model.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., & Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. International Journal of Production Research, 0(0), 1–10. doi:10.1080/00207543.2019.1600764

See Also

sim.es, ssarima, Distributions, orders

Examples

# Create 120 observations from ARIMA(1,1,1) with drift. Generate 100 time series of this kind.
x <- sim.ssarima(ar.orders=1,i.orders=1,ma.orders=1,obs=120,nsim=100,constant=TRUE)

# Generate similar thing for seasonal series of SARIMA(1,1,1)(0,0,2)_4
x <- sim.ssarima(ar.orders=c(1,0),i.orders=c(1,0),ma.orders=c(1,2),lags=c(1,4),
                 frequency=4,obs=80,nsim=100,constant=FALSE)

# Generate 10 series of high frequency data from SARIMA(1,0,2)_1(0,1,1)_7(1,0,1)_30
x <- sim.ssarima(ar.orders=c(1,0,1),i.orders=c(0,1,0),ma.orders=c(2,1,1),lags=c(1,7,30),
                 obs=360,nsim=10)

Simple Moving Average

Description

Function constructs state space simple moving average of predefined order

Usage

sma(y, order = NULL, ic = c("AICc", "AIC", "BIC", "BICc"), h = 10,
  holdout = FALSE, silent = TRUE, fast = TRUE, ...)

sma_old(y, order = NULL, ic = c("AICc", "AIC", "BIC", "BICc"), h = 10,
  holdout = FALSE, cumulative = FALSE, interval = c("none", "parametric",
  "likelihood", "semiparametric", "nonparametric"), level = 0.95,
  silent = c("all", "graph", "legend", "output", "none"), ...)

Arguments

Details

The function constructs AR model in the Single Source of Error state space form based on the idea that:

y_{t} = \frac{1}{n} \sum_{j=1}^n y_{t-j}

which is AR(n) process, that can be modelled using:

y_{t} = w' v_{t-1} + \epsilon_{t}

v_{t} = F v_{t-1} + g \epsilon_{t}

Where v_{t} is a state vector.

For some more information about the model and its implementation, see the vignette: vignette("sma","smooth")

Value

Object of class "smooth" is returned. It contains the list of the following values:

• model - the name of the estimated model.

• timeElapsed - time elapsed for the construction of the model.

• states - the matrix of the fuzzy components of ssarima, where rows correspond to time and cols to states.

• transition - matrix F.

• persistence - the persistence vector. This is the place, where smoothing parameters live.

• measurement - measurement vector of the model.

• order - order of moving average.

• initial - Initial state vector values.

• initialType - Type of initial values used.

• nParam - table with the number of estimated / provided parameters. If a previous model was reused, then its initials are reused and the number of provided parameters will take this into account.

• fitted - the fitted values.

• forecast - the point forecast.

• lower - the lower bound of prediction interval. When interval=FALSE then NA is returned.

• upper - the higher bound of prediction interval. When interval=FALSE then NA is returned.

• residuals - the residuals of the estimated model.

• errors - The matrix of 1 to h steps ahead errors. Only returned when the multistep losses are used and semiparametric interval is needed.

• s2 - variance of the residuals (taking degrees of freedom into account).

• interval - type of interval asked by user.

• level - confidence level for interval.

• cumulative - whether the produced forecast was cumulative or not.

• y - the original data.

• holdout - the holdout part of the original data.

• ICs - values of information criteria of the model. Includes AIC, AICc, BIC and BICc.

• logLik - log-likelihood of the function.

• lossValue - Cost function value.

• loss - Type of loss function used in the estimation.

• accuracy - vector of accuracy measures for the holdout sample. Includes: MPE, MAPE, SMAPE, MASE, sMAE, RelMAE, sMSE and Bias coefficient (based on complex numbers). This is available only when holdout=TRUE.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

• Svetunkov, I., & Petropoulos, F. (2017). Old dog, new tricks: a modelling view of simple moving averages. International Journal of Production Research, 7543(January), 1-14. doi:10.1080/00207543.2017.1380326

See Also

filter, adam, msarima

Examples

# SMA of specific order
ourModel <- sma(rnorm(118,100,3), order=12, h=18, holdout=TRUE)

# SMA of arbitrary order
ourModel <- sma(rnorm(118,100,3), h=18, holdout=TRUE)

plot(forecast(ourModel, h=18, interval="empirical"))

Combination of forecasts of state space models

Description

Function constructs ETS, SSARIMA, CES, GUM and SMA and combines their forecasts using IC weights.

Usage

smoothCombine(y, models = NULL, initial = c("backcasting", "optimal",
  "two-stage", "complete"), ic = c("AICc", "AIC", "BIC", "BICc"),
  loss = c("MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10,
  holdout = FALSE, cumulative = FALSE, interval = c("none", "prediction",
  "confidence", "simulated", "approximate", "semiparametric", "nonparametric",
  "empirical", "complete"), level = 0.95, bins = 200,
  intervalCombine = c("quantile", "probability"), bounds = c("usual",
  "admissible", "none"), silent = TRUE, xreg = NULL,
  regressors = c("use", "select"), initialX = NULL, ...)

Arguments

Details

The combination of these models using information criteria weights is possible because they are all formulated in Single Source of Error framework. Due to the the complexity of some of the models, the estimation process may take some time. So be patient.

The prediction interval are combined either probability-wise or quantile-wise (Lichtendahl et al., 2013), which may take extra time, because we need to produce all the distributions for all the models. This can be sped up with the smaller value for bins parameter, but the resulting interval may be imprecise.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Kolassa, S. (2011) Combining exponential smoothing forecasts using Akaike weights. International Journal of Forecasting, 27, pp 238 - 251.

• Svetunkov, I., Boylan, J.E., 2023b. Staying Positive: Challenges and Solutions in Using Pure Multiplicative ETS Models. IMA Journal of Management Mathematics. p. 403-425. doi:10.1093/imaman/dpad028

• Taylor, J.W. and Bunn, D.W. (1999) A Quantile Regression Approach to Generating Prediction Intervals. Management Science, Vol 45, No 2, pp 225-237.

• Lichtendahl Kenneth C., Jr., Grushka-Cockayne Yael, Winkler Robert L., (2013) Is It Better to Average Probabilities or Quantiles? Management Science 59(7):1594-1611. DOI: doi:10.1287/mnsc.1120.1667

See Also

es, auto.ssarima, auto.ces, auto.gum, sma

Examples

## Not run: ourModel <- smoothCombine(BJsales,interval="p")
plot(ourModel)
## End(Not run)

# models parameter accepts either previously estimated smoothCombine
# or a manually formed list of smooth models estimated in sample:
## Not run: smoothCombine(BJsales,models=ourModel)

## Not run: models <- list(es(BJsales), sma(BJsales))
smoothCombine(BJsales,models=models)
## End(Not run)

Function returns the ultimate answer to any question

Description

You need a description? So what?

Usage

sowhat(...)

Arguments

Details

You need details? So what?

Value

It doesn't return any value, only messages. So what?

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Sowhat?

• Sowhat?

• 42

See Also

Nowwhat (to be implemented),

Examples

x <- rnorm(10000,0,1);
sowhat(x);

sowhat("What's the meaning of life?")

sowhat("I don't have a girlfriend.")

State Space ARIMA

Description

Function constructs State Space ARIMA, estimating AR, MA terms and initial states.

Function selects the best State Space ARIMA based on information criteria, using fancy branch and bound mechanism. The resulting model can be not optimal in IC meaning, but it is usually reasonable.

Function constructs State Space ARIMA, estimating AR, MA terms and initial states.

Usage

ssarima(y, orders = list(ar = c(0), i = c(1), ma = c(1)), lags = c(1),
  constant = FALSE, arma = NULL, model = NULL,
  initial = c("backcasting", "optimal", "two-stage", "complete"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE",
  "MSCE"), h = 0, holdout = FALSE, bounds = c("admissible", "usual",
  "none"), silent = TRUE, xreg = NULL, regressors = c("use", "select",
  "adapt"), initialX = NULL, ...)

auto.ssarima(y, orders = list(ar = c(3, 3), i = c(2, 1), ma = c(3, 3)),
  lags = c(1, frequency(y)), fast = TRUE, constant = NULL,
  initial = c("backcasting", "optimal", "two-stage", "complete"),
  loss = c("likelihood", "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE",
  "MSCE"), ic = c("AICc", "AIC", "BIC", "BICc"), h = 0, holdout = FALSE,
  bounds = c("admissible", "usual", "none"), silent = TRUE, xreg = NULL,
  regressors = c("use", "select", "adapt"), ...)

ssarima_old(y, orders = list(ar = c(0), i = c(1), ma = c(1)), lags = c(1),
  constant = FALSE, AR = NULL, MA = NULL, initial = c("backcasting",
  "optimal"), ic = c("AICc", "AIC", "BIC", "BICc"), loss = c("likelihood",
  "MSE", "MAE", "HAM", "MSEh", "TMSE", "GTMSE", "MSCE"), h = 10,
  holdout = FALSE, bounds = c("admissible", "none"), silent = c("all",
  "graph", "legend", "output", "none"), xreg = NULL, regressors = c("use",
  "select"), initialX = NULL, ...)

Arguments

Details

The model, implemented in this function, is discussed in Svetunkov & Boylan (2019).

The basic ARIMA(p,d,q) used in the function has the following form:

(1 - B)^d (1 - a_1 B - a_2 B^2 - ... - a_p B^p) y_[t] = (1 + b_1 B + b_2 B^2 + ... + b_q B^q) \epsilon_[t] + c

where y_[t] is the actual values, \epsilon_[t] is the error term, a_i, b_j are the parameters for AR and MA respectively and c is the constant. In case of non-zero differences c acts as drift.

This model is then transformed into ARIMA in the Single Source of Error State space form (proposed in Snyder, 1985):

y_{t} = w' v_{t-l} + \epsilon_{t}

v_{t} = F v_{t-l} + g_t \epsilon_{t}

where v_{t} is the state vector (defined based on orders) and l is the vector of lags, w_t is the measurement vector (with explanatory variables if provided), F is the transition matrix, g_t is the persistence vector (which includes explanatory variables if they were used).

Due to the flexibility of the model, multiple seasonalities can be used. For example, something crazy like this can be constructed: SARIMA(1,1,1)(0,1,1)[24](2,0,1)[24*7](0,0,1)[24*30], but the estimation may take a lot of time... If you plan estimating a model with more than one seasonality, it is recommended to use msarima instead.

The model selection for SSARIMA is done by the auto.ssarima function.

For some more information about the model and its implementation, see the vignette: vignette("ssarima","smooth")

The function constructs bunch of ARIMAs in Single Source of Error state space form (see ssarima documentation) and selects the best one based on information criterion. The mechanism is described in Svetunkov & Boylan (2019).

Due to the flexibility of the model, multiple seasonalities can be used. For example, something crazy like this can be constructed: SARIMA(1,1,1)(0,1,1)[24](2,0,1)[24*7](0,0,1)[24*30], but the estimation may take a lot of time... It is recommended to use auto.msarima in cases with more than one seasonality and high frequencies.

For some more information about the model and its implementation, see the vignette: vignette("ssarima","smooth")

The model, implemented in this function, is discussed in Svetunkov & Boylan (2019).

The basic ARIMA(p,d,q) used in the function has the following form:

(1 - B)^d (1 - a_1 B - a_2 B^2 - ... - a_p B^p) y_[t] = (1 + b_1 B + b_2 B^2 + ... + b_q B^q) \epsilon_[t] + c

where y_[t] is the actual values, \epsilon_[t] is the error term, a_i, b_j are the parameters for AR and MA respectively and c is the constant. In case of non-zero differences c acts as drift.

This model is then transformed into ARIMA in the Single Source of Error State space form (proposed in Snyder, 1985):

y_{t} = o_{t} (w' v_{t-l} + x_t a_{t-1} + \epsilon_{t})

v_{t} = F v_{t-l} + g \epsilon_{t}

a_{t} = F_{X} a_{t-1} + g_{X} \epsilon_{t} / x_{t}

Where o_{t} is the Bernoulli distributed random variable (in case of normal data equal to 1), v_{t} is the state vector (defined based on orders) and l is the vector of lags, x_t is the vector of exogenous parameters. w is the measurement vector, F is the transition matrix, g is the persistence vector, a_t is the vector of parameters for exogenous variables, F_{X} is the transitionX matrix and g_{X} is the persistenceX matrix.

Due to the flexibility of the model, multiple seasonalities can be used. For example, something crazy like this can be constructed: SARIMA(1,1,1)(0,1,1)[24](2,0,1)[24*7](0,0,1)[24*30], but the estimation may take some finite time... If you plan estimating a model with more than one seasonality, it is recommended to consider doing it using msarima.

The model selection for SSARIMA is done by the auto.ssarima function.

For some more information about the model and its implementation, see the vignette: vignette("ssarima","smooth")

Value

Object of class "adam" is returned with similar elements to the adam function.

Object of class "smooth" is returned. See ssarima for details.

Object of class "smooth" is returned. It contains the list of the following values:

• model - the name of the estimated model.

• timeElapsed - time elapsed for the construction of the model.

• states - the matrix of the fuzzy components of ssarima, where rows correspond to time and cols to states.

• transition - matrix F.

• persistence - the persistence vector. This is the place, where smoothing parameters live.

• measurement - measurement vector of the model.

• AR - the matrix of coefficients of AR terms.

• I - the matrix of coefficients of I terms.

• MA - the matrix of coefficients of MA terms.

• constant - the value of the constant term.

• initialType - Type of the initial values used.

• initial - the initial values of the state vector (extracted from states).

• nParam - table with the number of estimated / provided parameters. If a previous model was reused, then its initials are reused and the number of provided parameters will take this into account.

• fitted - the fitted values.

• forecast - the point forecast.

• lower - the lower bound of prediction interval. When interval="none" then NA is returned.

• upper - the higher bound of prediction interval. When interval="none" then NA is returned.

• residuals - the residuals of the estimated model.

• errors - The matrix of 1 to h steps ahead errors. Only returned when the multistep losses are used and semiparametric interval is needed.

• s2 - variance of the residuals (taking degrees of freedom into account).

• interval - type of interval asked by user.

• level - confidence level for interval.

• cumulative - whether the produced forecast was cumulative or not.

• y - the original data.

• holdout - the holdout part of the original data.

• xreg - provided vector or matrix of exogenous variables. If regressors="s", then this value will contain only selected exogenous variables.

• initialX - initial values for parameters of exogenous variables.

• ICs - values of information criteria of the model. Includes AIC, AICc, BIC and BICc.

• logLik - log-likelihood of the function.

• lossValue - Cost function value.

• loss - Type of loss function used in the estimation.

• FI - Fisher Information. Equal to NULL if FI=FALSE or when FI is not provided at all.

• accuracy - vector of accuracy measures for the holdout sample. In case of non-intermittent data includes: MPE, MAPE, SMAPE, MASE, sMAE, RelMAE, sMSE and Bias coefficient (based on complex numbers). In case of intermittent data the set of errors will be: sMSE, sPIS, sCE (scaled cumulative error) and Bias coefficient. This is available only when holdout=TRUE.

• B - the vector of all the estimated parameters.

Author(s)

Ivan Svetunkov, ivan@svetunkov.com

References

• Svetunkov I. (2023) Smooth forecasting with the smooth package in R. arXiv:2301.01790. doi:10.48550/arXiv.2301.01790.

• Svetunkov I. (2015 - Inf) "smooth" package for R - series of posts about the underlying models and how to use them: https://openforecast.org/category/r-en/smooth/.

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., Boylan, J.E., 2023a. iETS: State Space Model for Intermittent Demand Forecastings. International Journal of Production Economics. 109013. doi:10.1016/j.ijpe.2023.109013

• Teunter R., Syntetos A., Babai Z. (2011). Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214, 606-615.

• Croston, J. (1972) Forecasting and stock control for intermittent demands. Operational Research Quarterly, 23(3), 289-303.

• Svetunkov, I., & Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. International Journal of Production Research, 0(0), 1–10. doi:10.1080/00207543.2019.1600764

• Svetunkov, I., 2023. Smooth Forecasting with the Smooth Package in R. arXiv. doi:10.48550/arXiv.2301.01790

• Snyder, R. D., 1985. Recursive Estimation of Dynamic Linear Models. Journal of the Royal Statistical Society, Series B (Methodological) 47 (2), 272-276.

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag. doi:10.1007/978-3-540-71918-2.

• Svetunkov, I., & Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. International Journal of Production Research, 0(0), 1–10. doi:10.1080/00207543.2019.1600764

See Also

auto.ssarima, auto.msarima, adam, es, ces

es, ces, sim.es, gum, ssarima

auto.ssarima, orders, msarima, auto.msarima, sim.ssarima, adam

Examples

# ARIMA(1,1,1) fitted to some data
ourModel <- ssarima(rnorm(118,100,3),orders=list(ar=c(1),i=c(1),ma=c(1)),lags=c(1))

# Model with the same lags and orders, applied to a different data
ssarima(rnorm(118,100,3),orders=orders(ourModel),lags=lags(ourModel))

# The same model applied to a different data
ssarima(rnorm(118,100,3),model=ourModel)

# Example of SARIMA(2,0,0)(1,0,0)[4]
ssarima(rnorm(118,100,3),orders=list(ar=c(2,1)),lags=c(1,4))

# SARIMA(1,1,1)(0,0,1)[4] with different initialisations
ssarima(rnorm(118,100,3),orders=list(ar=c(1),i=c(1),ma=c(1,1)),
        lags=c(1,4),h=18,holdout=TRUE,initial="backcasting")


set.seed(41)
x <- rnorm(118,100,3)

# The best ARIMA for the data
ourModel <- auto.ssarima(x,orders=list(ar=c(2,1),i=c(1,1),ma=c(2,1)),lags=c(1,12),
                                   h=18,holdout=TRUE)

# The other one using optimised states
auto.ssarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
                       initial="two",h=18,holdout=TRUE)

summary(ourModel)
forecast(ourModel)
plot(forecast(ourModel))


# ARIMA(1,1,1) fitted to some data
ourModel <- ssarima_old(rnorm(118,100,3),orders=list(ar=c(1),i=c(1),ma=c(1)),lags=c(1),h=18,
                             holdout=TRUE)

# Model with the same lags and orders, applied to a different data
ssarima_old(rnorm(118,100,3),orders=orders(ourModel),lags=lags(ourModel),h=18,holdout=TRUE)

# The same model applied to a different data
ssarima_old(rnorm(118,100,3),model=ourModel,h=18,holdout=TRUE)

# SARIMA(0,1,1) with exogenous variables
ssarima_old(rnorm(118,100,3),orders=list(i=1,ma=1),h=18,holdout=TRUE,xreg=c(1:118))

summary(ourModel)
forecast(ourModel)
plot(forecast(ourModel))