#' ADAM is Dynamic Adaptive Model
#'
#' Function constructs an advanced Single Source of Error model, based on ETS
#' taxonomy and ARIMA elements
#'
#' Function estimates ADAM in a form of the Single Source of Error state space
#' model of the following type:
#'
#' \deqn{y_{t} = o_t (w(v_{t-l}) + h(x_t, a_{t-1}) + r(v_{t-l}) \epsilon_{t})}
#'
#' \deqn{v_{t} = f(v_{t-l}, a_{t-1}) + g(v_{t-l}, a_{t-1}, x_{t}) \epsilon_{t}}
#'
#' Where \eqn{o_{t}} is the Bernoulli distributed random variable (in case of
#' normal data it equals to 1 for all observations), \eqn{v_{t}} is the state
#' vector and \eqn{l} is the vector of lags, \eqn{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 \eqn{a_t} is the vector of parameters for exogenous variables.
#' Finally, \eqn{\epsilon_{t}} is the error term.
#'
#' The implemented model allows introducing several seasonal states and supports
#' intermittent data via the \code{occurrence} variable.
#'
#' The error term \eqn{\epsilon_t} can follow different distributions, which
#' are regulated via the \code{distribution} parameter. This includes:
#' \enumerate{
#' \item \code{default} - Normal distribution is used for the Additive error models,
#' Inverse Gaussian is used for the Multiplicative error models.
#' \item \link[stats]{dnorm} - Normal distribution,
#' \item \link[greybox]{dlaplace} - Laplace distribution,
#' \item \link[greybox]{ds} - S distribution,
#' \item \link[gnorm]{dgnorm} - Generalised Normal distribution,
#' \item \link[stats]{dlogis} - Logistic Distribution,
#' \item \link[stats]{dt} - T distribution,
#' \item \link[greybox]{dalaplace} - Asymmetric Laplace distribution,
#' \item \link[stats]{dlnorm} - Log normal distribution,
#' \item dllaplace - Log Laplace distribution,
#' \item dls - Log S distribution,
#' \item dlgnorm - Log Generalised Normal distribution,
# \item \link[greybox]{dbcnorm} - Box-Cox normal distribution,
#' \item \link[statmod]{dinvgauss} - Inverse Gaussian distribution,
#' }
#'
#' For some more information about the model and its implementation, see the
#' vignette: \code{vignette("adam","smooth")}.
#'
#' The function \code{auto.adam()} tries out models with the specified
#' distributions and returns the one with the most suitable one.
#'
#' @template ssAuthor
#' @template ssKeywords
#'
#' @template ssGeneralRef
#' @template ssIntermittentRef
#' @template ssETSRef
#' @template ssIntervalsRef
#'
#' @param y Vector, containing data needed to be forecasted. If a matrix is
#' provided, then the first column is used as a response variable, while the rest
#' of the matrix is used as a set of explanatory variables.
#' @param model The type of ETS model. The first letter stands for the type of
#' the error term ("A" or "M"), the second (and sometimes the third as well) is for
#' the trend ("N", "A", "Ad", "M" or "Md"), and the last one is for the type of
#' seasonality ("N", "A" or "M"). In case of several lags, the seasonal components
#' are assumed to be the same. The model is then printed out as
#' ETS(M,Ad,M)[m1,m2,...], where m1, m2, ... are the lags specified by the
#' \code{lags} parameter.
#' There are several options for the \code{model} besides the conventional ones,
#' which rely on information criteria:
#' \enumerate{
#' \item \code{model="ZZZ"} means that the model will be selected based on the
#' chosen information criteria type. The Branch and Bound is used in the process.
#' \item \code{model="XXX"} means that only additive components are tested, using
#' Branch and Bound.
#' \item \code{model="YYY"} implies selecting between multiplicative components.
#' \item \code{model="CCC"} trigers the combination of forecasts of models using
#' information criteria weights (Kolassa, 2011).
#' \item combinations between these four and the classical components are also
#' accepted. For example, \code{model="CAY"} will combine models with additive
#' trend and either none or multiplicative seasonality.
#' \item \code{model="PPP"} will produce the selection between pure additive and
#' pure multiplicative models. "P" stands for "Pure". This cannot be mixed with
#' other types of components.
#' \item \code{model="FFF"} will select between all the 30 types of models. "F"
#' stands for "Full". This cannot be mixed with other types of components.
#' \item The parameter \code{model} can also be a vector of names of models for a
#' finer tuning (pool of models). For example, \code{model=c("ANN","AAA")} will
#' estimate only two models and select the best of them.
#' }
#'
#' Also, \code{model} can accept a previously estimated adam and use all
#' its parameters.
#'
#' Keep in mind that model selection with "Z" components uses Branch and Bound
#' algorithm and may skip some models that could have slightly smaller
#' information criteria. If you want to do a exhaustive search, you would need
#' to list all the models to check as a vector.
#'
#' The default value is set to \code{"ZXZ"}, because the multiplicative trend is explosive
#' and dangerous. It should be used only for each separate time series, not for the
#' atomated predictions for big datasets.
#'
#' @param lags Defines lags for the corresponding components. All components
#' count, starting from level, so ETS(M,M,M) model for monthly data will have
#' \code{lags=c(1,1,12)}. However, the function will also accept \code{lags=c(12)},
#' assuming that the lags 1 were dropped.
#' @param orders The order of ARIMA to be included in the model. This should be passed
#' either as a vector (in which case the non-seasonal ARIMA is assumed) or as a list of
#' a type \code{orders=list(ar=c(p,P),i=c(d,D),ma=c(q,Q))}, in which case the \code{lags}
#' variable is used in order to determine the seasonality m. See \link[smooth]{msarima}
#' for details. Note that ARIMA here is mainly treated as an addition to the ETS model and
#' does not allow constant / drift.
#'
#' In case of \code{auto.adam()} function, \code{orders} accepts one more parameters:
#' \code{orders=list(select=FALSE)}. If \code{TRUE}, then the function will select the most
#' appropriate order using a mechanism similar to \code{auto.msarima()}. The values
#' \code{list(ar=...,i=...,ma=...)} specify the maximum orders to check in this case.
#' @param formula Formula to use in case of explanatory variables. If \code{NULL},
#' then all the variables are used as is. Only considered if \code{xreg} is not
#' \code{NULL} and \code{xregDo="use"}.
#' @param distribution what density function to assume for the error term. The full
#' name of the distribution should be provided, starting with the letter "d" -
#' "density". The names align with the names of distribution functions in R.
#' For example, see \link[stats]{dnorm}. For detailed explanation of available
#' distributions, see vignette in greybox package: \code{vignette("greybox","alm")}.
#' @param loss The type of Loss Function used in optimization. \code{loss} can
#' be:
#' \itemize{
#' \item \code{likelihood} - the model is estimated via the maximisation of the
#' likelihood of the function specified in \code{distribution};
#' \item \code{MSE} (Mean Squared Error),
#' \item \code{MAE} (Mean Absolute Error),
#' \item \code{HAM} (Half Absolute Moment),
#' \item \code{LASSO} - use LASSO to shrink the parameters of the model;
#' \item \code{RIDGE} - use RIDGE to shrink the parameters of the model;
#' \item \code{TMSE} - Trace Mean Squared Error,
#' \item \code{GTMSE} - Geometric Trace Mean Squared Error,
#' \item \code{MSEh} - optimisation using only h-steps ahead error,
#' \item \code{MSCE} - Mean Squared Cumulative Error.
#' }
#' In case of LASSO / RIDGE, the variables are not normalised prior to the estimation,
#' but the parameters are divided by the mean values of explanatory variables.
#'
#' Note that model selection and combination works properly only for the default
#' \code{loss="likelihood"}.
#'
#' Furthermore, just for fun the absolute and half analogues of multistep estimators
#' are available: \code{MAEh}, \code{TMAE}, \code{GTMAE}, \code{MACE},
#' \code{HAMh}, \code{THAM}, \code{GTHAM}, \code{CHAM}.
#'
#' Last but not least, user can provide their own function here as well, making sure
#' that it accepts parameters \code{actual}, \code{fitted} and \code{B}. Here is an
#' example:
#' \code{lossFunction <- function(actual, fitted, B) return(mean(abs(actual-fitted)))}
#' \code{loss=lossFunction}
#' @param h The forecast horizon. Mainly needed for the multistep loss functions.
#' @param holdout Logical. If \code{TRUE}, then the holdout of the size \code{h}
#' is taken from the data (can be used for the model testing purposes).
#' @param persistence Persistence vector \eqn{g}, containing smoothing
#' parameters. If \code{NULL}, then estimated. Can be also passed as a names list of
#' the type: \code{persistence=list(level=0.1, trend=0.05, seasonal=c(0.1,0.2),
#' xreg=c(0.1,0.2))}. Dropping some elements from the named list will make the function
#' estimate them. e.g. if you don't specify seasonal in the persistence for the ETS(M,N,M)
#' model, it will be estimated.
#' @param phi Value of damping parameter. If \code{NULL} then it is estimated.
#' Only applicable for damped-trend models.
#' @param initial Can be either character or a list, or a vector of initial states.
#' If it is character, then it can be \code{"optimal"}, meaning that the initial
#' states are optimised, or \code{"backcasting"}, meaning that the initials are
#' produced using backcasting procedure (advised for data with high frequency). In
#' case of the list, it is recommended to use the named one and to provide those
#' initial components that are available. For example:
#' \code{initial=list(level=1000,trend=10,seasonal=list(c(1,2),c(1,2,3,4)),
#' arima=1,xreg=100)}. If some of the components are needed by the model, but are
#' not provided in the list, they will be estimated. If the vector is provided,
#' then it is expected that the components will be provided one after another
#' without any gaps.
#' @param arma Either the named list or a vector with AR / MA parameters ordered lag-wise.
#' The number of elements should correspond to the specified orders e.g.
#' \code{orders=list(ar=c(1,1),ma=c(1,1)), lags=c(1,4), arma=list(ar=c(0.9,0.8),ma=c(-0.3,0.3))}
#' @param occurrence The type of model used in probability estimation. Can be
#' \code{"none"} - none,
#' \code{"fixed"} - constant probability,
#' \code{"general"} - the general Beta model with two parameters,
#' \code{"odds-ratio"} - the Odds-ratio model with b=1 in Beta distribution,
#' \code{"inverse-odds-ratio"} - the model with a=1 in Beta distribution,
#' \code{"direct"} - the TSB-like (Teunter et al., 2011) probability update
#' mechanism a+b=1,
#' \code{"auto"} - the automatically selected type of occurrence model.
#'
#' The type of model used in the occurrence is equal to the one provided in the
#' \code{model} parameter.
#'
#' Also, a model produced using \link[smooth]{oes} or \link[greybox]{alm} function
#' can be used here.
#' @param ic The information criterion to use in the model selection / combination
#' procedure.
#' @param bounds The type of bounds for the persistence to use in the model
#' estimation. Can be either \code{admissible} - guaranteeing the stability of the
#' model, \code{traditional} - restricting the values with (0, 1) or \code{none} - no
#' restrictions (potentially dangerous).
#' @param xreg The vector (either numeric or time series) or the matrix (or
#' data.frame / data.table) of exogenous variables that should be included in the
#' model. If matrix is included than columns should contain variables and rows -
#' observations.
#' Note that \code{xreg} should have number of observations equal to
#' the length of the response variable \code{y}. If it is not equal, then the
#' function will either trim or extrapolate the data.
#' @param xregDo The variable defines what to do with the provided xreg:
#' \code{"use"} means that all of the data should be used, while
#' \code{"select"} means that a selection using \code{ic} should be done,
#' \code{"adapt"} will trigger the mechanism of time varying parameters for the
#' explanatory variables.
#' @param silent Specifies, whether to provide the progress of the function or not.
#' If \code{TRUE}, then the function will print what it does and how much it has
#' already done.
#' @param ... Other non-documented parameters. For example, \code{FI=TRUE} will
#' make the function also produce Fisher Information matrix, which then can be
#' used to calculated variances of smoothing parameters and initial states of
#' the model. This is used in the \link[stats]{vcov} method.
#' Starting values of parameters can be passed via \code{B}, while the upper and lower
#' bounds should be passed in \code{ub} and \code{lb} respectively. In this case they
#' will be used for optimisation. These values should have the length equal
#' to the number of parameters to estimate in the following order:
#' \enumerate{
#' \item All smoothing parameters (for the states and then for the explanatory variables);
#' \item Damping parameter (if needed);
#' \item ARMA parameters;
#' \item All the initial values (for the states and then for the explanatory variables).
#' }
#' You can also pass parameters to the optimiser in order to fine tune its work:
#' \itemize{
#' \item \code{maxeval} - maximum number of evaluations to carry out (default is 20 per
#' estimated parameter and at least 1000 if pure ARIMA is considered);
#' \item \code{maxtime} - stop, when the optimisation time (in seconds) exceeds this;
#' \item \code{xtol_rel} - the relative precision of the optimiser (the default is 1E-6);
#' \item \code{xtol_abs} - the absolute precision of the optimiser (the default is 1E-8);
#' \item \code{ftol_rel} - the stopping criterion in case of the relative change in the loss
#' function (the default is 1E-4);
#' \item \code{ftol_abs} - the stopping criterion in case of the absolute change in the loss
#' function (the default is 0 - not used);
#' \item \code{algorithm} - the algorithm to use in optimisation
#' (by default, \code{"NLOPT_LN_SBPLX"} is used);
#' \item \code{print_level} - the level of output for the optimiser (0 by default).
#' If equal to 41, then the detailed results of the optimisation are returned.
#' }
#' You can read more about these parameters by running the function
#' \link[nloptr]{nloptr.print.options}.
#' Finally, the parameter \code{lambda} for LASSO / RIDGE, Asymmetric Laplace and df
#' of Student's distribution can be provided here as well.
#'
#' @return Object of class "adam" is returned. It contains the list of the
#' following values:
#' \itemize{
#' \item \code{model} - the name of the constructed model,
#' \item \code{timeElapsed} - the time elapsed for the estimation of the model,
#' \item \code{y} - the in-sample part of the data used for the training of the model,
#' \item \code{holdout} - the holdout part of the data, excluded for purposes of model evaluation,
#' \item \code{fitted} - the vector of fitted values,
#' \item \code{residuals} - the vector of residuals,
#' \item \code{forecast} - the point forecast for h steps ahead (by default NA is returned),
#' \item \code{states} - the matrix of states with observations in rows and states in columns,
#' \item \code{persisten} - the vector of smoothing parameters,
#' \item \code{phi} - the value of damping parameter,
#' \item \code{transition} - the transition matrix,
#' \item \code{measurement} - the measurement matrix with observations in rows and state elements
#' in columns,
#' \item \code{initial} - the named list of initial values, including level, trend, seasonal, ARIMA
#' and xreg components,
#' \item \code{initialEstimated} - the named vector, defining which of the initials were estimated in
#' the model,
#' \item \code{initialType} - the type of initialisation used ("optimal" / "backcasting" / "provided"),
#' \item \code{orders} - the orders of ARIMA used in the estimation,
#' \item \code{arma} - the list of AR / MA parameters used in the model,
#' \item \code{nParam} - the matrix of the estimated / provided parameters,
#' \item \code{occurrence} - the oes model used for the occurrence part of the model,
#' \item \code{xreg} - the matrix of explanatory variables after all expansions and transformations,
#' \item \code{formula} - the formula used for the explanatory variables expansion,
#' \item \code{loss} - the type of loss function used in the estimation,
#' \item \code{lossValue} - the value of that loss function,
#' \item \code{logLik} - the value of the log-likelihood,
#' \item \code{distribution} - the distribution function used in the calculation of the likelihood,
#' \item \code{scale} - the value of the scale parameter,
#' \item \code{lambda} - the value of the parameter used in LASSO / dalaplace / dt,
#' \item \code{B} - the vector of all estimated parameters,
#' \item \code{lags} - the vector of lags used in the model construction,
#' \item \code{lagsAll} - the vector of the internal lags used in the model.
#' }
#'
#' @seealso \code{\link[forecast]{ets}, \link[smooth]{es}}
#'
#' @examples
#'
#' # Model selection using a specified pool of models
#' ourModel <- adam(rnorm(100,100,10), model=c("ANN","ANA","AAA"), lags=c(5,10))
#'
#' summary(ourModel)
#' forecast(ourModel)
#' plot(forecast(ourModel))
#'
#' # Model combination using a specified pool
#' ourModel <- adam(rnorm(100,100,10), model=c("ANN","AAN","MNN","CCC"), lags=c(5,10))
#'
#' @importFrom forecast forecast na.interp
#' @importFrom greybox dlaplace dalaplace ds stepwise alm is.occurrence polyprod
#' @importFrom stats dnorm dlogis dt dlnorm frequency
#' @importFrom statmod dinvgauss
#' @importFrom gnorm dgnorm
#' @importFrom nloptr nloptr
#' @importFrom pracma hessian
#' @importFrom zoo zoo
#' @useDynLib mes
#' @rdname adam
#' @export adam
adam <- function(y, model="ZXZ", lags=c(1,frequency(y)), orders=list(ar=c(0),i=c(0),ma=c(0)), formula=NULL,
distribution=c("default","dnorm","dlaplace","ds","dgnorm","dlogis","dt","dalaplace",
"dlnorm","dllaplace","dls","dlgnorm","dinvgauss"),
loss=c("likelihood","MSE","MAE","HAM","LASSO","RIDGE","MSEh","TMSE","GTMSE","MSCE"),
h=0, holdout=FALSE,
persistence=NULL, phi=NULL, initial=c("optimal","backcasting"), arma=NULL,
occurrence=c("none","auto","fixed","general","odds-ratio","inverse-odds-ratio","direct"),
ic=c("AICc","AIC","BIC","BICc"), bounds=c("usual","admissible","none"),
xreg=NULL, xregDo=c("use","select","adapt"), silent=TRUE, ...){
# Copyright (C) 2019 - Inf Ivan Svetunkov
# Start measuring the time of calculations
startTime <- Sys.time();
ellipsis <- list(...);
# If a previous model is provided as a model, write down the variables
if(is.adam(model) || is.adam.sim(model)){
# If this is the simulated data, extract the parameters
# if(is.adam.sim(model) & !is.null(dim(model$data))){
# warning("The provided model has several submodels. Choosing a random one.",call.=FALSE);
# i <- round(runif(1,1:length(model$persistence)));
# persistence <- model$persistence[,i];
# initial <- model$initial[,i];
# initialSeason <- model$initialSeason[,i];
# if(any(model$iprob!=1)){
# occurrence <- "a";
# }
# }
# else{
persistence <- model$persistence;
initial <- model$initial;
initialEstimated <- model$initialEstimated;
distribution <- model$distribution;
loss <- model$loss;
persistence <- model$persistence;
phi <- model$phi;
if(model$initialType!="backcasting"){
initial <- model$initial;
}
else{
initial <- "b";
}
occurrence <- model$occurrence;
ic <- model$ic;
bounds <- model$bounds;
ellipsis$lambda <- model$lambda;
ellipsis$B <- model$B;
CFValue <- model$lossValue;
logLikADAMValue <- logLik(model);
if(is.null(xreg)){
xreg <- model$xreg;
xregDo <- model$xregDo;
}
else{
if(is.null(model$xreg)){
xreg <- NULL;
}
else{
if(ncol(xreg)!=ncol(model$xreg)){
xreg <- xreg[,colnames(model$xreg)];
}
}
xregDo <- model$xregDo;
}
formula <- formula(model);
# Parameters of the original ARIMA model
lags <- lags(model);
orders <- orders(model);
arma <- model$arma;
model <- modelType(model);
modelDo <- "use";
# if(any(unlist(gregexpr("C",model))!=-1)){
# initial <- "o";
# }
}
else if(inherits(model,"ets")){
# Extract smoothing parameters
i <- 1;
lags <- 1;
persistence <- coef(model)[i];
if(model$components[2]!="N"){
i <- i+1;
persistence <- c(persistence,coef(model)[i]);
if(model$components[3]!="N"){
i <- i+1;
persistence <- c(persistence,coef(model)[i]);
}
}
else{
if(model$components[3]!="N"){
i <- i+1;
persistence <- c(persistence,coef(model)[i]);
}
}
# Damping parameter
if(model$components[4]=="TRUE"){
i <- i+1;
phi <- coef(model)[i];
}
# Initials
i <- i+1;
initial <- coef(model)[i];
# Initial for the trend
if(model$components[2]!="N"){
i <- i+1;
lags <- c(lags,1);
initial <- c(initial,coef(model)[i]);
}
# Initials of seasonal component
if(model$components[3]!="N"){
if(model$components[2]!="N"){
initial <- c(initial,rev(model$states[1,-c(1:2)]));
}
else{
initial <- c(initial,rev(model$states[1,-c(1)]));
}
lags <- c(lags,model$m);
}
model <- modelType(model);
distribution <- "dnorm";
loss <- "likelihood";
modelDo <- "use"
}
else if(is.character(model)){
modelDo <- "";
# Everything is okay
}
else{
modelDo <- "";
warning("A model of an unknown class was provided. Switching to 'ZZZ'.",call.=FALSE);
model <- "ZZZ";
}
# paste0() is needed in order to get rid of potential issues with names
yName <- paste0(deparse(substitute(y)),collapse="");
#### Check the parameters of the function and create variables based on them ####
checkerReturn <- parametersChecker(y, model, lags, formula, orders, arma,
persistence, phi, initial,
distribution, loss, h, holdout, occurrence, ic, bounds,
xreg, xregDo, yName,
silent, modelDo, ParentEnvironment=environment(), ellipsis, fast=FALSE);
#### Return regression if it is pure ####
if(is.alm(checkerReturn)){
obsInSample <- nobs(checkerReturn);
nParam <- length(checkerReturn$coefficient);
if(!is.null(dim(y))){
xreg <- y[,-1,drop=FALSE];
}
modelReturned <- list(model="Regression");
modelReturned$timeElapsed <- Sys.time()-startTime;
modelReturned$y <- actuals(checkerReturn);
if(holdout){
if(is.null(dim(y))){
yHoldout <- y[obsInSample+c(1:h)];
}
else{
yHoldout <- y[obsInSample+c(1:h),1];
}
modelReturned$holdout <- yHoldout;
}
else{
modelReturned$holdout <- NULL;
}
# Extract indeces from the data
yIndex <- try(time(y),silent=TRUE);
# If we cannot extract time, do something
if(inherits(yIndex,"try-error")){
if(!is.data.frame(y) && !is.null(dim(y))){
yIndex <- as.POSIXct(rownames(y));
}
else if(is.data.frame(y)){
yIndex <- c(1:nrow(y));
}
else{
yIndex <- c(1:length(y));
}
}
# Prepare fitted, residuals and the forecasts
if(inherits(y,"zoo")){
modelReturned$y <- zoo(modelReturned$y, order.by=yIndex[1:obsInSample]);
modelReturned$fitted <- zoo(fitted(checkerReturn), order.by=yIndex[1:obsInSample]);
modelReturned$residuals <- zoo(residuals(checkerReturn), order.by=yIndex[1:obsInSample]);
if(h>0){
modelReturned$forecast <- zoo(forecast(checkerReturn,h=h,newdata=tail(xreg,h))$mean,
order.by=yIndex[obsInSample+1:h]);
}
else{
modelReturned$forecast <- zoo(NA, order.by=yIndex[obsInSample+1]);
}
modelReturned$states <- zoo(matrix(coef(checkerReturn), obsInSample+1, nParam, byrow=TRUE,
dimnames=list(NULL, names(coef(checkerReturn)))),
order.by=c(yIndex[1]-diff(yIndex[1:2]),yIndex[1:obsInSample]));
}
else{
yFrequency <- frequency(y);
modelReturned$y <- ts(modelReturned$y, start=yIndex[1], frequency=yFrequency);
modelReturned$fitted <- ts(fitted(checkerReturn), start=yIndex[1], frequency=yFrequency);
modelReturned$residuals <- ts(residuals(checkerReturn), start=yIndex[1], frequency=yFrequency);
if(h>0){
modelReturned$forecast <- ts(forecast(checkerReturn,h=h,newdata=tail(xreg,h))$mean,
start=yIndex[obsInSample+1], frequency=yFrequency);
}
else{
modelReturned$forecast <- ts(NA, start=yIndex[obsInSample]+diff(yIndex[1:2]), frequency=yFrequency);
}
modelReturned$states <- ts(matrix(coef(checkerReturn), obsInSample+1, nParam, byrow=TRUE,
dimnames=list(NULL, names(coef(checkerReturn)))),
start=yIndex-diff(yIndex[1:2]), frequency=yFrequency);
}
modelReturned$persistence <- rep(0, nParam);
names(modelReturned$persistence) <- paste0("delta",c(1:nParam));
modelReturned$phi <- 1;
modelReturned$transition <- diag(nParam);
modelReturned$measurement <- checkerReturn$data;
modelReturned$measurement[,1] <- 1;
colnames(modelReturned$measurement) <- colnames(modelReturned$states);
modelReturned$initial <- list(xreg=coef(checkerReturn));
modelReturned$initialType <- "optimal";
modelReturned$initialEstimated <- TRUE;
names(modelReturned$initialEstimated) <- "xreg";
modelReturned$orders <- list(ar=0,i=0,ma=0);
modelReturned$arma <- NULL;
# Number of estimated parameters
parametersNumber <- matrix(0,2,4,
dimnames=list(c("Estimated","Provided"),
c("nParamInternal","nParamXreg","nParamOccurrence","nParamAll")));
parametersNumber[1,2] <- nParam;
if(is.occurrence(checkerReturn$occurrence)){
parametersNumber[1,3] <- nParam;
}
parametersNumber[1,4] <- sum(parametersNumber[1,1:3]);
modelReturned$nParam <- parametersNumber;
modelReturned$occurrence <- checkerReturn$occurrence;
modelReturned$xreg <- checkerReturn$data[,-1,drop=FALSE];
modelReturned$formula <- formula(checkerReturn);
modelReturned$xregDo <- "use";
modelReturned$loss <- checkerReturn$loss;
modelReturned$lossValue <- checkerReturn$lossValue;
modelReturned$lossFunction <- checkerReturn$lossFunction;
modelReturned$logLik <- logLik(checkerReturn);
modelReturned$distribution <- checkerReturn$distribution;
modelReturned$scale <- checkerReturn$scale;
modelReturned$lambda <- checkerReturn$other;
modelReturned$B <- coef(checkerReturn);
modelReturned$lags <- 1;
modelReturned$lagsAll <- rep(1,nParam);
modelReturned$FI <- checkerReturn$FI;
if(holdout){
modelReturned$accuracy <- measures(modelReturned$holdout,modelReturned$forecast,modelReturned$y)
}
else{
modelReturned$accuracy <- NULL;
}
class(modelReturned) <- c("adam","smooth");
if(!silent){
plot(modelReturned,7)
}
return(modelReturned);
}
# Remove xreg if it was provided, just to preserve some memory
rm(xreg);
#### The function creates the technical variables (lags etc) based on the type of the model ####
architector <- function(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal,
xregNumber, obsInSample, initialType,
arimaModel, lagsModelARIMA, xregModel){
# If there is ETS
if(etsModel){
modelIsTrendy <- Ttype!="N";
if(modelIsTrendy){
# Make lags (1, 1)
lagsModel <- matrix(c(1,1),ncol=1);
componentsNamesETS <- c("level","trend");
}
else{
# Make lags (1, ...)
lagsModel <- matrix(c(1),ncol=1);
componentsNamesETS <- c("level");
}
modelIsSeasonal <- Stype!="N";
if(modelIsSeasonal){
# If the lags are for the non-seasonal model
lagsModel <- matrix(c(lagsModel,lagsModelSeasonal),ncol=1);
componentsNumberETSSeasonal <- length(lagsModelSeasonal);
if(componentsNumberETSSeasonal>1){
componentsNamesETS <- c(componentsNamesETS,paste0("seasonal",c(1:componentsNumberETSSeasonal)));
}
else{
componentsNamesETS <- c(componentsNamesETS,"seasonal");
}
}
else{
componentsNumberETSSeasonal <- 0;
}
lagsModelAll <- lagsModel;
componentsNumberETS <- length(lagsModel);
}
else{
modelIsTrendy <- modelIsSeasonal <- FALSE;
componentsNumberETS <- componentsNumberETSSeasonal <- 0;
componentsNamesETS <- NULL;
lagsModelAll <- lagsModel <- NULL;
}
# If there is ARIMA
if(arimaModel){
lagsModelAll <- matrix(c(lagsModel,lagsModelARIMA), ncol=1);
}
# If there are xreg
if(xregModel){
lagsModelAll <- matrix(c(lagsModelAll,rep(1,xregNumber)), ncol=1);
}
lagsModelMax <- max(lagsModelAll);
# Define the number of cols that should be in the matvt
obsStates <- obsInSample + lagsModelMax*switch(initialType,
"backcasting"=2,
1);
return(list(lagsModel=lagsModel,lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax,
componentsNumberETS=componentsNumberETS, componentsNumberETSSeasonal=componentsNumberETSSeasonal,
componentsNumberETSNonSeasonal=componentsNumberETS-componentsNumberETSSeasonal,
componentsNamesETS=componentsNamesETS, obsStates=obsStates, modelIsTrendy=modelIsTrendy,
modelIsSeasonal=modelIsSeasonal));
}
#### The function creates the necessary matrices based on the model and provided parameters ####
# This is needed in order to initialise the estimation
creator <- function(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
lags, lagsModel, lagsModelARIMA, lagsModelAll, lagsModelMax,
obsStates, obsInSample, obsAll, componentsNumberETS, componentsNumberETSSeasonal,
componentsNamesETS, otLogical, yInSample,
# Persistence
persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate, persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi,
# Initials
initialType, initialEstimate,
initialLevel, initialLevelEstimate, initialTrend, initialTrendEstimate,
initialSeasonal, initialSeasonalEstimate,
initialArima, initialArimaEstimate, initialArimaNumber,
initialXregEstimate, initialXregProvided,
# ARIMA elements
arimaModel, arRequired, iRequired, maRequired, armaParameters,
arOrders, iOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA,
# Explanatory variables
xregModel, xregModelInitials, xregData, xregNumber, xregNames){
# Matrix of states. Time in columns, components in rows
matVt <- matrix(NA, componentsNumberETS+componentsNumberARIMA+xregNumber, obsStates,
dimnames=list(c(componentsNamesETS,componentsNamesARIMA,xregNames),NULL));
# Measurement rowvector
matWt <- matrix(1, obsAll, componentsNumberETS+componentsNumberARIMA+xregNumber,
dimnames=list(NULL,c(componentsNamesETS,componentsNamesARIMA,xregNames)));
# If xreg are provided, then fill in the respective values in Wt vector
if(xregModel){
matWt[,componentsNumberETS+componentsNumberARIMA+1:xregNumber] <- xregData;
}
# Transition matrix
matF <- diag(componentsNumberETS+componentsNumberARIMA+xregNumber);
# Persistence vector
vecG <- matrix(0, componentsNumberETS+componentsNumberARIMA+xregNumber, 1,
dimnames=list(c(componentsNamesETS,componentsNamesARIMA,xregNames),NULL));
j <- 0;
# ETS model
if(etsModel){
j <- j+1;
rownames(vecG)[j] <- "alpha";
if(!persistenceLevelEstimate){
vecG[j,] <- persistenceLevel;
}
if(modelIsTrendy){
j <- j+1;
rownames(vecG)[j] <- "beta";
if(!persistenceTrendEstimate){
vecG[j,] <- persistenceTrend;
}
}
if(modelIsSeasonal){
if(!all(persistenceSeasonalEstimate)){
vecG[j+which(!persistenceSeasonalEstimate),] <- persistenceSeasonal;
}
if(componentsNumberETSSeasonal>1){
rownames(vecG)[j+c(1:componentsNumberETSSeasonal)] <- paste0("gamma",c(1:componentsNumberETSSeasonal));
}
else{
rownames(vecG)[j+1] <- "gamma";
}
j <- j+componentsNumberETSSeasonal;
}
}
# ARIMA model, names for persistence
if(arimaModel){
# Remove diagonal from the ARIMA part of the matrix
matF[j+1:componentsNumberARIMA,j+1:componentsNumberARIMA] <- 0;
if(componentsNumberARIMA>1){
rownames(vecG)[j+1:componentsNumberARIMA] <- paste0("psi",c(1:componentsNumberARIMA));
}
else{
rownames(vecG)[j+1:componentsNumberARIMA] <- "psi";
}
j <- j+componentsNumberARIMA;
}
# Regression
if(xregModel){
if(persistenceXregProvided && !persistenceXregEstimate){
vecG[j+1:xregNumber,] <- persistenceXreg;
}
rownames(vecG)[j+1:xregNumber] <- paste0("delta",c(1:xregNumber));
}
# Damping parameter value
if(etsModel && modelIsTrendy){
matF[1,2] <- phi;
matF[2,2] <- phi;
matWt[,2] <- phi;
}
# If the arma parameters were provided, fill in the persistence
if(arimaModel && (!arEstimate && !maEstimate)){
# Call polynomial
arimaPolynomials <- polynomialiser(NULL, arOrders, iOrders, maOrders,
arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags);
# Fill in the transition matrix
if(nrow(nonZeroARI)>0){
matF[componentsNumberETS+nonZeroARI[,2],componentsNumberETS+nonZeroARI[,2]] <-
-arimaPolynomials$ariPolynomial[nonZeroARI[,1]];
}
# Fill in the persistence vector
if(nrow(nonZeroARI)>0){
vecG[componentsNumberETS+nonZeroARI[,2]] <- -arimaPolynomials$ariPolynomial[nonZeroARI[,1]];
}
if(nrow(nonZeroMA)>0){
vecG[componentsNumberETS+nonZeroMA[,2]] <- vecG[componentsNumberETS+nonZeroMA[,2]] +
arimaPolynomials$maPolynomial[nonZeroMA[,1]];
}
}
else{
arimaPolynomials <- NULL;
}
# ETS model, initial state
# If something needs to be estimated...
if(etsModel){
if(initialEstimate){
# For the seasonal models
if(modelIsSeasonal){
if(obsNonzero>=lagsModelMax*2){
# If either Etype or Stype are multiplicative, do multiplicative decomposition
decompositionType <- c("additive","multiplicative")[any(c(Etype,Stype)=="M")+1];
yDecomposition <- msdecompose(yInSample, lags[lags!=1], type=decompositionType);
j <- 1;
# level
if(initialLevelEstimate){
matVt[j,1:lagsModelMax] <- mean(yInSample[1:lagsModelMax]);
if(xregModel){
if(Etype=="A"){
matVt[j,1:lagsModelMax] <- matVt[j,1:lagsModelMax] -
as.vector(xregModelInitials[[1]]$initialXreg %*% xregData[1,]);
}
else{
matVt[j,1:lagsModelMax] <- matVt[j,1:lagsModelMax] /
as.vector(exp(xregModelInitials[[2]]$initialXreg %*% xregData[1,]));
}
}
}
else{
matVt[j,1:lagsModelMax] <- initialLevel;
}
j <- j+1;
# If trend is needed
if(modelIsTrendy){
if(initialTrendEstimate){
if(Ttype=="A" && Stype=="M"){
if(initialLevelEstimate){
# level fix
matVt[j-1,1:lagsModelMax] <- exp(mean(log(yInSample[otLogical][1:lagsModelMax])));
}
# trend
matVt[j,1:lagsModelMax] <- prod(yDecomposition$initial)-yDecomposition$initial[1];
}
else if(Ttype=="M" && Stype=="A"){
if(initialLevelEstimate){
# level fix
matVt[j-1,1:lagsModelMax] <- exp(mean(log(yInSample[otLogical][1:lagsModelMax])));
}
# trend
matVt[j,1:lagsModelMax] <- sum(yDecomposition$initial)/yDecomposition$initial[1];
}
else{
# trend
matVt[j,1:lagsModelMax] <- yDecomposition$initial[2];
}
# This is a failsafe for multiplicative trend models, so that the thing does not explode
if(Ttype=="M" && matVt[j,1:lagsModelMax]>1.1){
matVt[j,1:lagsModelMax] <- 1;
}
}
else{
matVt[j,1:lagsModelMax] <- initialTrend;
}
j <- j+1;
}
#### Seasonal components
# For pure models use stuff as is
if(all(c(Etype,Stype)=="A") || all(c(Etype,Stype)=="M") ||
(Etype=="A" & Stype=="M")){
for(i in 1:componentsNumberETSSeasonal){
if(initialSeasonalEstimate[i]){
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <- yDecomposition$seasonal[[i]];
# Renormalise the initial seasons
if(Stype=="A"){
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] -
mean(matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]]);
}
else{
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] /
exp(mean(log(matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]])));
}
}
else{
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <- initialSeasonal[[i]];
}
}
}
# For mixed models use a different set of initials
else if(Etype=="M" && Stype=="A"){
for(i in 1:componentsNumberETSSeasonal){
if(initialSeasonalEstimate[i]){
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+
1:lagsModel[i+j-1]] <- log(yDecomposition$seasonal[[i]])*min(yInSample[otLogical]);
# Renormalise the initial seasons
if(Stype=="A"){
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] -
mean(matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]]);
}
else{
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] /
exp(mean(log(matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]])));
}
}
else{
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <- initialSeasonal[[i]];
}
}
}
}
else{
# If either Etype or Stype are multiplicative, do multiplicative decomposition
j <- 1;
# level
if(initialLevelEstimate){
matVt[j,1:lagsModelMax] <- mean(yInSample[1:lagsModelMax]);
if(xregModel){
if(Etype=="A"){
matVt[j,1:lagsModelMax] <- matVt[j,1:lagsModelMax] -
as.vector(xregModelInitials[[1]]$initialXreg %*% xregData[1,]);
}
else{
matVt[j,1:lagsModelMax] <- matVt[j,1:lagsModelMax] /
as.vector(exp(xregModelInitials[[2]]$initialXreg %*% xregData[1,]));
}
}
}
else{
matVt[j,1:lagsModelMax] <- initialLevel;
}
j <- j+1;
if(modelIsTrendy){
if(initialTrendEstimate){
if(Ttype=="A"){
# trend
matVt[j,1:lagsModelMax] <- yInSample[2]-yInSample[1];
}
else if(Ttype=="M"){
if(initialLevelEstimate){
# level fix
matVt[j-1,1:lagsModelMax] <- exp(mean(log(yInSample[otLogical][1:lagsModelMax])));
}
# trend
matVt[j,1:lagsModelMax] <- yInSample[2]/yInSample[1];
}
# This is a failsafe for multiplicative trend models, so that the thing does not explode
if(Ttype=="M" && matVt[j,1:lagsModelMax]>1.1){
matVt[j,1:lagsModelMax] <- 1;
}
}
else{
matVt[j,1:lagsModelMax] <- initialTrend;
}
j <- j+1;
}
#### Seasonal components
# For pure models use stuff as is
if(Stype=="A"){
for(i in 1:componentsNumberETSSeasonal){
if(initialSeasonalEstimate[i]){
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
yInSample[1:lagsModel[i+j-1]]-matVt[1,1];
# Renormalise the initial seasons
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] -
mean(matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]]);
}
else{
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <- initialSeasonal[[i]];
}
}
}
# For mixed models use a different set of initials
else{
for(i in 1:componentsNumberETSSeasonal){
if(initialSeasonalEstimate[i]){
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
yInSample[1:lagsModel[i+j-1]]/matVt[1,1];
# Renormalise the initial seasons
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <-
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] /
exp(mean(log(matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]])));
}
else{
matVt[i+j-1,(lagsModelMax-lagsModel[i+j-1])+1:lagsModel[i+j-1]] <- initialSeasonal[[i]];
}
}
}
}
}
# Non-seasonal models
else{
# level
if(initialLevelEstimate){
matVt[1,lagsModelMax] <- mean(yInSample[1:max(lagsModelMax,ceiling(obsInSample*0.2))]);
if(xregModel){
if(Etype=="A"){
matVt[1,lagsModelMax] <- matVt[1,lagsModelMax] -
as.vector(xregModelInitials[[1]]$initialXreg %*% xregData[1,]);
}
else{
matVt[1,lagsModelMax] <- matVt[1,lagsModelMax] /
as.vector(exp(xregModelInitials[[2]]$initialXreg %*% xregData[1,]));
}
}
}
else{
matVt[1,lagsModelMax] <- initialLevel;
}
if(modelIsTrendy){
if(initialTrendEstimate){
matVt[2,lagsModelMax] <- switch(Ttype,
"A" = mean(diff(yInSample[1:max(lagsModelMax+1,ceiling(obsInSample*0.2))]),na.rm=TRUE),
"M" = exp(mean(diff(log(yInSample[otLogical])),na.rm=TRUE)));
}
else{
matVt[2,lagsModelMax] <- initialTrend;
}
}
}
if(initialLevelEstimate && Etype=="M" && matVt[1,lagsModelMax]==0){
matVt[1,lagsModelMax] <- mean(yInSample);
}
}
# Else, insert the provided ones... make sure that this is not a backcasting
else if(!initialEstimate && initialType=="provided"){
j <- 1;
matVt[j,1:lagsModelMax] <- initialLevel;
if(modelIsTrendy){
j <- j+1;
matVt[j,1:lagsModelMax] <- initialTrend;
}
if(modelIsSeasonal){
for(i in 1:componentsNumberETSSeasonal){
matVt[j+i,(lagsModelMax-lagsModel[j+i])+1:lagsModel[j+i]] <- initialSeasonal[[i]];
}
}
j <- j+componentsNumberETSSeasonal;
}
}
# If ARIMA orders are specified, prepare initials
if(arimaModel){
if(initialArimaEstimate){
matVt[componentsNumberETS+1:componentsNumberARIMA,
1:lagsModelARIMA[componentsNumberARIMA]+(lagsModelMax-lagsModelARIMA[componentsNumberARIMA])] <-
switch(Etype, "A"=0, "M"=1);
# If this is just ARIMA with optimisation, refine the initials
if(!etsModel && initialType!="backcasting"){
# This is needed in order to make the initial components more realistic
# matVt[1:componentsNumberARIMA,
# 1:lagsModelARIMA[componentsNumberARIMA]+(lagsModelMax-lagsModelARIMA[componentsNumberARIMA])] <-
# yInSample[lagsModelMax:1];
arimaPolynomials <- polynomialiser(rep(0.1,sum(c(arOrders,maOrders))), arOrders, iOrders, maOrders,
arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags);
if(nrow(nonZeroARI)>0 && nrow(nonZeroARI)>=nrow(nonZeroMA)){
matVt[componentsNumberETS+nonZeroARI[,2],
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*%
t(matVt[componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber]) /
tail(arimaPolynomials$ariPolynomial,1),
"M"=exp(arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*%
t(log(matVt[componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber])) /
tail(arimaPolynomials$ariPolynomial,1)));
}
else{
matVt[componentsNumberETS+nonZeroMA[,2],
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*%
t(matVt[componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber]) /
tail(arimaPolynomials$maPolynomial,1),
"M"=exp(arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*%
t(log(matVt[componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber])) /
tail(arimaPolynomials$maPolynomial,1)));
}
}
}
else{
# Fill in the matrix with 0 / 1, just in case if the state will not be updated anymore
matVt[componentsNumberETS+1:componentsNumberARIMA,
1:lagsModelARIMA[componentsNumberARIMA]+(lagsModelMax-lagsModelARIMA[componentsNumberARIMA])] <-
switch(Etype, "A"=0, "M"=1);
# Insert the provided initials
matVt[componentsNumberETS+componentsNumberARIMA, 1:lagsModelARIMA[componentsNumberARIMA]+
(lagsModelMax-lagsModelARIMA[componentsNumberARIMA])] <-
initialArima[1:initialArimaNumber];
# If only AR is needed, but provided or if both are needed, but provided
if(((arRequired && !arEstimate) && !maRequired) ||
((arRequired && !arEstimate) && (maRequired && !maEstimate)) ||
(iRequired && !arEstimate && !maEstimate)){
matVt[componentsNumberETS+nonZeroARI[,2],lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*% t(initialArima[1:initialArimaNumber]) /
tail(arimaPolynomials$ariPolynomial,1),
"M"=exp(arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*% t(log(initialArima[1:initialArimaNumber])) /
tail(arimaPolynomials$ariPolynomial,1)));
}
# If only MA is needed, but provided
else if(((maRequired && !maEstimate) && !arRequired)){
matVt[componentsNumberETS+nonZeroMA[,2],lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*% t(initialArima[1:initialArimaNumber]) /
tail(arimaPolynomials$maPolynomial,1),
"M"=exp(arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*% t(log(initialArima[1:initialArimaNumber])) /
tail(arimaPolynomials$maPolynomial,1)));
}
}
}
# Fill in the initials for xreg
if(xregModel){
if(Etype=="A" || initialXregProvided || is.null(xregModelInitials[[2]])){
matVt[componentsNumberETS+componentsNumberARIMA+1:xregNumber,1:lagsModelMax] <- xregModelInitials[[1]]$initialXreg;
}
else{
matVt[componentsNumberETS+componentsNumberARIMA+1:xregNumber,1:lagsModelMax] <- xregModelInitials[[2]]$initialXreg;
}
}
return(list(matVt=matVt, matWt=matWt, matF=matF, vecG=vecG, arimaPolynomials=arimaPolynomials));
}
#### ARI and MA polynomials function ####
if(arimaModel){
polynomialiser <- function(B, arOrders, iOrders, maOrders,
arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags){
# Number of parameters that we have
nParamAR <- sum(arOrders);
nParamMA <- sum(maOrders);
# Matrices with parameters
arParameters <- matrix(0, max(arOrders * lags) + 1, length(arOrders));
iParameters <- matrix(0, max(iOrders * lags) + 1, length(iOrders));
maParameters <- matrix(0, max(maOrders * lags) + 1, length(maOrders));
# The first element is always 1
arParameters[1,] <- iParameters[1,] <- maParameters[1,] <- 1;
# nParam is used for B
nParam <- 1;
# armanParam is used for the provided arma parameters
armanParam <- 1;
# Fill in the matrices with the provided parameters
for(i in 1:length(lags)){
if(arOrders[i]*lags[i]!=0){
if(arEstimate){
arParameters[1+(1:arOrders[i])*lags[i],i] <- -B[nParam+c(1:arOrders[i])-1];
nParam <- nParam + arOrders[i];
}
else if(!arEstimate && arRequired){
arParameters[1+(1:arOrders[i])*lags[i],i] <- -armaParameters[armanParam+c(1:arOrders[i])-1];
armanParam <- armanParam + arOrders[i];
}
}
if(iOrders[i]*lags[i] != 0){
iParameters[1+lags[i],i] <- -1;
}
if(maOrders[i]*lags[i]!=0){
if(maEstimate){
maParameters[1+(1:maOrders[i])*lags[i],i] <- B[nParam+c(1:maOrders[i])-1];
nParam <- nParam + maOrders[i];
}
else if(!maEstimate && maRequired){
maParameters[1+(1:maOrders[i])*lags[i],i] <- armaParameters[armanParam+c(1:maOrders[i])-1];
armanParam <- armanParam + maOrders[i];
}
}
}
# Vectors of polynomials for the ARIMA
arPolynomial <- vector("numeric", sum(arOrders * lags) + 1);
iPolynomial <- vector("numeric", sum(iOrders * lags) + 1);
maPolynomial <- vector("numeric", sum(maOrders * lags) + 1);
ariPolynomial <- vector("numeric", sum(arOrders * lags) + sum(iOrders * lags) + 1);
# Fill in the first polynomials
arPolynomial[0:(arOrders[1]*lags[1])+1] <- arParameters[0:(arOrders[1]*lags[1])+1,1];
iPolynomial[0:(iOrders[1]*lags[1])+1] <- iParameters[0:(iOrders[1]*lags[1])+1,1];
maPolynomial[0:(maOrders[1]*lags[1])+1] <- maParameters[0:(maOrders[1]*lags[1])+1,1];
index1 <- 0
index2 <- 0;
# Fill in all the other polynomials
for(i in 1:length(lags)){
if(i!=1){
if(arOrders[i]>0){
index1[] <- tail(which(arPolynomial!=0),1);
index2[] <- tail(which(arParameters[,i]!=0),1);
arPolynomial[1:(index1+index2-1)] <- polyprod(arPolynomial[1:index1], arParameters[1:index2,i]);
}
if(maOrders[i]>0){
index1[] <- tail(which(maPolynomial!=0),1);
index2[] <- tail(which(maParameters[,i]!=0),1);
maPolynomial[1:(index1+index2-1)] <- polyprod(maPolynomial[1:index1], maParameters[1:index2,i]);
}
if(iOrders[i]>0){
index1[] <- tail(which(iPolynomial!=0),1);
index2[] <- tail(which(iParameters[,i]!=0),1);
iPolynomial[1:(index1+index2-1)] <- polyprod(iPolynomial[1:index1], iParameters[1:index2,i]);
}
}
# This part takes the power of (1-B)^D
if(iOrders[i]>1){
for(j in 2:iOrders[i]){
index1[] <- tail(which(iPolynomial!=0),1);
index2[] <- tail(which(iParameters[,i]!=0),1);
iPolynomial[1:(index1+index2-1)] = polyprod(iPolynomial[1:index1], iParameters[1:index2,i]);
}
}
}
# ARI polynomials
ariPolynomial[] <- polyprod(arPolynomial, iPolynomial);
return(list(arPolynomial=arPolynomial,iPolynomial=iPolynomial,
ariPolynomial=ariPolynomial,maPolynomial=maPolynomial));
}
}
#### The function fills in the existing matrices with values of A ####
# This is needed in order to do the estimation and the fit
filler <- function(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETS, componentsNumberETSNonSeasonal,
componentsNumberETSSeasonal, componentsNumberARIMA,
lags, lagsModel, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arEstimate, maEstimate, arOrders, iOrders, maOrders,
arRequired, maRequired, armaParameters,
nonZeroARI, nonZeroMA, arimaPolynomials,
xregModel, xregNumber){
j <- 0;
# Fill in persistence
if(persistenceEstimate){
# Persistence of ETS
if(etsModel){
i <- 1;
# alpha
if(persistenceLevelEstimate){
j[] <- j+1;
vecG[i] <- B[j];
}
# beta
if(modelIsTrendy){
i[] <- 2;
if(persistenceTrendEstimate){
j[] <- j+1;
vecG[i] <- B[j];
}
}
# gamma1, gamma2, ...
if(modelIsSeasonal){
if(any(persistenceSeasonalEstimate)){
vecG[i+which(persistenceSeasonalEstimate)] <- B[j+c(1:sum(persistenceSeasonalEstimate))];
j[] <- j+sum(persistenceSeasonalEstimate);
}
i[] <- componentsNumberETS;
}
}
# Persistence of xreg
if(xregModel && persistenceXregEstimate){
vecG[j+componentsNumberARIMA+1:xregNumber] <- B[j+1:xregNumber];
j[] <- j+xregNumber;
}
}
# Damping parameter
if(etsModel && phiEstimate){
j[] <- j+1;
matWt[,2] <- B[j];
matF[1:2,2] <- B[j];
}
# ARMA parameters. This goes before xreg in persistence
if(arimaModel){
# Call the function returning ARI and MA polynomials
arimaPolynomials <- polynomialiser(B[j+1:sum(c(arOrders*arEstimate,maOrders*maEstimate))], arOrders, iOrders, maOrders,
arRequired, maRequired, arEstimate, maEstimate, armaParameters, lags);
# Fill in the transition matrix
if(nrow(nonZeroARI)>0){
matF[componentsNumberETS+nonZeroARI[,2],componentsNumberETS+1:componentsNumberARIMA] <-
-arimaPolynomials$ariPolynomial[nonZeroARI[,1]];
}
# Fill in the persistence vector
if(nrow(nonZeroARI)>0){
vecG[componentsNumberETS+nonZeroARI[,2]] <- -arimaPolynomials$ariPolynomial[nonZeroARI[,1]];
}
if(nrow(nonZeroMA)>0){
vecG[componentsNumberETS+nonZeroMA[,2]] <- vecG[componentsNumberETS+nonZeroMA[,2]] +
arimaPolynomials$maPolynomial[nonZeroMA[,1]];
}
j[] <- j+sum(c(arOrders*arEstimate,maOrders*maEstimate));
}
# Initials of ETS if something needs to be estimated
if(etsModel && (initialType!="backcasting") && initialEstimate){
i <- 1;
if(initialLevelEstimate){
j[] <- j+1;
matVt[i,1:lagsModelMax] <- B[j];
}
i[] <- i+1;
if(modelIsTrendy && initialTrendEstimate){
j[] <- j+1;
matVt[i,1:lagsModelMax] <- B[j];
i[] <- i+1;
}
if(modelIsSeasonal && any(initialSeasonalEstimate)){
for(k in 1:componentsNumberETSSeasonal){
if(initialSeasonalEstimate[k]){
matVt[componentsNumberETSNonSeasonal+k,
(lagsModelMax-lagsModel[componentsNumberETSNonSeasonal+k])+
2:lagsModel[componentsNumberETSNonSeasonal+k]-1] <-
B[j+2:(lagsModel[componentsNumberETSNonSeasonal+k])-1];
matVt[componentsNumberETSNonSeasonal+k,
(lagsModelMax-lagsModel[componentsNumberETSNonSeasonal+k])+
lagsModel[componentsNumberETSNonSeasonal+k]] <-
switch(Stype,
"A"=-sum(B[j+2:(lagsModel[componentsNumberETSNonSeasonal+k])-1]),
"M"=1/prod(B[j+2:(lagsModel[componentsNumberETSNonSeasonal+k])-1]));
j[] <- j+lagsModel[componentsNumberETSNonSeasonal+k]-1;
}
}
}
}
# Initials of ARIMA
if(arimaModel){
if((initialType!="backcasting") && initialArimaEstimate){
if(nrow(nonZeroARI)>0 && nrow(nonZeroARI)>=nrow(nonZeroMA)){
matVt[componentsNumberETS+componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <- B[j+1:initialArimaNumber];
matVt[componentsNumberETS+nonZeroARI[,2],
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*% t(B[j+1:initialArimaNumber]) /
tail(arimaPolynomials$ariPolynomial,1),
"M"=exp(arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*% t(log(B[j+1:initialArimaNumber])) /
tail(arimaPolynomials$ariPolynomial,1)));
}
else{
matVt[componentsNumberETS+componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <- B[j+1:initialArimaNumber];
matVt[componentsNumberETS+nonZeroMA[,2],
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*% t(B[j+1:initialArimaNumber]) /
tail(arimaPolynomials$maPolynomial,1),
"M"=exp(arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*% t(log(B[j+1:initialArimaNumber])) /
tail(arimaPolynomials$maPolynomial,1)));
}
j[] <- j+componentsNumberARIMA;
}
# This is needed in order to propagate initials of ARIMA to all components
else if(any(c(arEstimate,maEstimate))){
if(nrow(nonZeroARI)>0 && nrow(nonZeroARI)>=nrow(nonZeroMA)){
matVt[componentsNumberETS+nonZeroARI[,2],
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"= arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*%
t(matVt[componentsNumberETS+componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber]) /
tail(arimaPolynomials$ariPolynomial,1),
"M"=exp(arimaPolynomials$ariPolynomial[nonZeroARI[,1]] %*%
t(log(matVt[componentsNumberETS+componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber])) /
tail(arimaPolynomials$ariPolynomial,1)));
}
else{
matVt[componentsNumberETS+nonZeroMA[,2],
lagsModelMax-initialArimaNumber+1:initialArimaNumber] <-
switch(Etype,
"A"=arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*%
t(matVt[componentsNumberETS+componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber]) /
tail(arimaPolynomials$maPolynomial,1),
"M"=exp(arimaPolynomials$maPolynomial[nonZeroMA[,1]] %*%
t(log(matVt[componentsNumberETS+componentsNumberARIMA,
lagsModelMax-initialArimaNumber+1:initialArimaNumber])) /
tail(arimaPolynomials$maPolynomial,1)));
}
}
}
# Initials of the xreg
if(xregModel && (initialType!="backcasting") && initialEstimate && initialXregEstimate){
matVt[componentsNumberETS+componentsNumberARIMA+1:xregNumber,1:lagsModelMax] <- B[j+1:xregNumber];
}
return(list(matVt=matVt, matWt=matWt, matF=matF, vecG=vecG, arimaPolynomials=arimaPolynomials));
}
#### The function initialises the vector B for ETS ####
initialiser <- function(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETSNonSeasonal, componentsNumberETSSeasonal, componentsNumberETS,
lags, lagsModel, lagsModelSeasonal, lagsModelARIMA, lagsModelMax,
matVt,
# persistence values
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
# initials
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
# ARIMA elements
arimaModel, arRequired, maRequired, arEstimate, maEstimate, arOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA, initialArimaNumber,
# Explanatory variables
xregModel, xregNumber, lambdaEstimate){
# The vector of logicals for persistence elements
persistenceEstimateVector <- c(persistenceLevelEstimate,modelIsTrendy&persistenceTrendEstimate,
modelIsSeasonal&persistenceSeasonalEstimate);
# The order:
# Persistence of states and for xreg, phi, AR and MA parameters, initials, initialsARIMA, initials for xreg
B <- Bl <- Bu <- vector("numeric",
# Values of the persistence vector + phi
etsModel*(persistenceLevelEstimate + modelIsTrendy*persistenceTrendEstimate +
modelIsSeasonal*sum(persistenceSeasonalEstimate) + phiEstimate) +
xregModel*persistenceXregEstimate*xregNumber +
# AR and MA values
arimaModel*(arEstimate*sum(arOrders)+maEstimate*sum(maOrders)) +
# initials of ETS
etsModel*(initialType!="backcasting")*
(initialLevelEstimate +
(modelIsTrendy*initialTrendEstimate) +
(modelIsSeasonal*sum(initialSeasonalEstimate*(lagsModelSeasonal-1)))) +
# initials of ARIMA
(initialType!="backcasting")*arimaModel*initialArimaNumber*initialArimaEstimate +
# initials of xreg
(initialType!="backcasting")*xregModel*initialXregEstimate*xregNumber+
lambdaEstimate);
j <- 0;
if(etsModel){
# Fill in persistence
if(persistenceEstimate){
if(any(c(Etype,Ttype,Stype)=="M")){
# A special type of model which is not safe: AAM, MAA, MAM
if((Etype=="A" && Ttype=="A" && Stype=="M") || (Etype=="A" && Ttype=="M" && Stype=="A") ||
((initialType=="backcasting") &&
((Etype=="M" && Ttype=="A" && Stype=="A") || (Etype=="M" && Ttype=="A" && Stype=="M")))){
B[1:sum(persistenceEstimateVector)] <-
c(0.01,0,rep(0,componentsNumberETSSeasonal))[which(persistenceEstimateVector)];
}
# MMA is the worst. Set everything to zero and see if anything can be done...
else if((Etype=="M" && Ttype=="M" && Stype=="A")){
B[1:sum(persistenceEstimateVector)] <-
c(0,0,rep(0,componentsNumberETSSeasonal))[which(persistenceEstimateVector)];
}
else if(Etype=="M" && Ttype=="A"){
if(initialType=="backcasting"){
B[1:sum(persistenceEstimateVector)] <-
c(0.1,0,rep(0.11,componentsNumberETSSeasonal))[which(persistenceEstimateVector)];
}
else{
B[1:sum(persistenceEstimateVector)] <-
c(0.1,0.05,rep(0.11,componentsNumberETSSeasonal))[which(persistenceEstimateVector)];
}
}
else{
B[1:sum(persistenceEstimateVector)] <-
c(0.1,0.05,rep(0.11,componentsNumberETSSeasonal))[which(persistenceEstimateVector)];
}
}
else{
B[1:sum(persistenceEstimateVector)] <-
c(0.1,0.05,rep(0.11,componentsNumberETSSeasonal))[which(persistenceEstimateVector)];
}
Bl[1:sum(persistenceEstimateVector)] <- rep(-5, sum(persistenceEstimateVector));
Bu[1:sum(persistenceEstimateVector)] <- rep(5, sum(persistenceEstimateVector));
# Names for B
if(persistenceLevelEstimate){
j[] <- j+1
names(B)[j] <- "alpha";
}
if(modelIsTrendy && persistenceTrendEstimate){
j[] <- j+1
names(B)[j] <- "beta";
}
if(modelIsSeasonal && any(persistenceSeasonalEstimate)){
if(componentsNumberETSSeasonal>1){
names(B)[j+c(1:sum(persistenceSeasonalEstimate))] <-
paste0("gamma",c(1:componentsNumberETSSeasonal));
}
else{
names(B)[j+1] <- "gamma";
}
j[] <- j+sum(persistenceSeasonalEstimate);
}
}
}
# Persistence if xreg is provided
if(xregModel && persistenceXregEstimate){
B[j+1:xregNumber] <- rep(0.01,xregNumber);
Bl[j+1:xregNumber] <- rep(-5, xregNumber);
Bu[j+1:xregNumber] <- rep(5, xregNumber);
names(B)[j+1:xregNumber] <- paste0("delta",c(1:xregNumber));
j[] <- j+xregNumber;
}
# Damping parameter
if(etsModel && phiEstimate){
j[] <- j+1;
B[j] <- 0.95;
names(B)[j] <- "phi";
Bl[j] <- 0;
Bu[j] <- 1;
}
# ARIMA parameters (AR / MA)
if(arimaModel){
# These are filled in lags-wise
if(any(c(arEstimate,maEstimate))){
for(i in 1:length(lags)){
if(arRequired && arEstimate && arOrders[i]>0){
B[j+c(1:arOrders[i])] <- rep(0.1,arOrders[i]);
Bl[j+c(1:arOrders[i])] <- -5;
Bu[j+c(1:arOrders[i])] <- 5;
names(B)[j+1:arOrders[i]] <- paste0("phi",1:arOrders[i],"[",lags[i],"]");
j[] <- j + arOrders[i];
}
if(maRequired && maEstimate && maOrders[i]>0){
B[j+c(1:maOrders[i])] <- rep(0.1,maOrders[i]);
Bl[j+c(1:maOrders[i])] <- -5;
Bu[j+c(1:maOrders[i])] <- 5;
names(B)[j+1:maOrders[i]] <- paste0("theta",1:maOrders[i],"[",lags[i],"]");
j[] <- j + maOrders[i];
}
}
}
}
# Initials
if(etsModel && initialType!="backcasting" && initialEstimate){
if(initialLevelEstimate){
j[] <- j+1;
B[j] <- matVt[1,lagsModelMax];
names(B)[j] <- "level";
if(Etype=="A"){
Bl[j] <- -Inf;
Bu[j] <- Inf;
}
else{
Bl[j] <- 0;
Bu[j] <- Inf;
}
}
if(modelIsTrendy && initialTrendEstimate){
j[] <- j+1;
B[j] <- matVt[2,lagsModelMax];
names(B)[j] <- "trend";
if(Ttype=="A"){
Bl[j] <- -Inf;
Bu[j] <- Inf;
}
else{
Bl[j] <- 0;
Bu[j] <- Inf;
}
}
if(modelIsSeasonal && any(initialSeasonalEstimate)){
if(componentsNumberETSSeasonal>1){
for(k in 1:componentsNumberETSSeasonal){
if(initialSeasonalEstimate[k]){
# -1 is needed in order to remove the redundant seasonal element (normalisation)
B[j+2:lagsModel[componentsNumberETSNonSeasonal+k]-1] <-
matVt[componentsNumberETSNonSeasonal+k,
(lagsModelMax-lagsModel[componentsNumberETSNonSeasonal+k])+
2:lagsModel[componentsNumberETSNonSeasonal+k]-1];
names(B)[j+2:(lagsModel[componentsNumberETSNonSeasonal+k])-1] <-
paste0("seasonal",k,"_",2:lagsModel[componentsNumberETSNonSeasonal+k]-1);
if(Stype=="A"){
Bl[j+2:lagsModel[componentsNumberETSNonSeasonal+k]-1] <- -Inf;
Bu[j+2:lagsModel[componentsNumberETSNonSeasonal+k]-1] <- Inf;
}
else{
Bl[j+2:lagsModel[componentsNumberETSNonSeasonal+k]-1] <- 0;
Bu[j+2:lagsModel[componentsNumberETSNonSeasonal+k]-1] <- Inf;
}
j[] <- j+(lagsModelSeasonal[k]-1);
}
}
}
else{
# -1 is needed in order to remove the redundant seasonal element (normalisation)
B[j+2:(lagsModel[componentsNumberETS])-1] <- matVt[componentsNumberETS,(lagsModelMax-lagsModel[componentsNumberETS])+
2:lagsModel[componentsNumberETS]-1];
names(B)[j+2:(lagsModel[componentsNumberETS])-1] <- paste0("seasonal_",2:lagsModel[componentsNumberETS]-1);
if(Stype=="A"){
Bl[j+2:(lagsModel[componentsNumberETS])-1] <- -Inf;
Bu[j+2:(lagsModel[componentsNumberETS])-1] <- Inf;
}
else{
Bl[j+2:(lagsModel[componentsNumberETS])-1] <- 0;
Bu[j+2:(lagsModel[componentsNumberETS])-1] <- Inf;
}
j[] <- j+(lagsModel[componentsNumberETS]-1);
}
}
}
# ARIMA initials
if(initialType!="backcasting" && arimaModel && initialArimaEstimate){
B[j+1:initialArimaNumber] <- tail(matVt[componentsNumberETS+componentsNumberARIMA,1:lagsModelMax],initialArimaNumber);
names(B)[j+1:initialArimaNumber] <- paste0("ARIMAState",1:initialArimaNumber);
j[] <- j+initialArimaNumber;
}
# Initials of the xreg
if(initialType!="backcasting" && initialXregEstimate){
B[j+1:xregNumber] <- matVt[componentsNumberETS+componentsNumberARIMA+1:xregNumber,lagsModelMax];
names(B)[j+1:xregNumber] <- rownames(matVt)[componentsNumberETS+componentsNumberARIMA+1:xregNumber];
if(Etype=="A"){
Bl[j+1:xregNumber] <- -Inf;
Bu[j+1:xregNumber] <- Inf;
}
else{
Bl[j+1:xregNumber] <- -Inf;
Bu[j+1:xregNumber] <- Inf;
}
j[] <- j+xregNumber;
}
# Add lambda if it is needed
if(lambdaEstimate){
j[] <- j+1;
B[j] <- lambda;
names(B)[j] <- "lambda";
Bl[j] <- 1e-10;
Bu[j] <- Inf;
}
return(list(B=B,Bl=Bl,Bu=Bu));
}
##### Function returns scale parameter for the provided parameters #####
scaler <- function(distribution, Etype, errors, yFitted, obsInSample, lambda){
# as.complex() is needed in order to make the optimiser work in exotic cases
scale <- switch(distribution,
"dnorm"=sqrt(sum(errors^2)/obsInSample),
"dlaplace"=sum(abs(errors))/obsInSample,
"ds"=sum(sqrt(abs(errors))) / (obsInSample*2),
"dgnorm"=(lambda*sum(abs(errors)^lambda)/obsInSample)^{1/lambda},
"dlogis"=sqrt(sum(errors^2)/obsInSample * 3 / pi^2),
"dt"=sqrt(sum(errors^2)/obsInSample),
"dalaplace"=sum(errors*(lambda-(errors<=0)*1))/obsInSample,
"dlnorm"=switch(Etype,
"A"=Re(sqrt(sum(log(as.complex(1+errors/yFitted))^2)/obsInSample)),
"M"=sqrt(sum(log(1+errors)^2)/obsInSample)),
"dllaplace"=switch(Etype,
"A"=Re(sum(abs(log(as.complex(1+errors/yFitted))))/obsInSample),
"M"=sum(abs(log(1+errors))/obsInSample)),
"dls"=switch(Etype,
"A"=Re(sum(sqrt(abs(log(as.complex(1+errors/yFitted))))/obsInSample)),
"M"=sum(sqrt(abs(log(1+errors)))/obsInSample)),
"dlgnorm"=switch(Etype,
"A"=Re((lambda*sum(abs(log(as.complex(1+errors/yFitted)))^lambda)/obsInSample)^{1/lambda}),
"M"=(lambda*sum(abs(log(as.complex(1+errors)))^lambda)/obsInSample)^{1/lambda}),
"dinvgauss"=switch(Etype,
"A"=sum((errors/yFitted)^2/(1+errors/yFitted))/obsInSample,
"M"=sum((errors)^2/(1+errors))/obsInSample),
# "M"=mean((errors)^2/(1+errors))),
);
return(scale);
}
##### Cost Function for ETS #####
CF <- function(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, yInSample,
ot, otLogical, occurrenceModel, obsInSample,
componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal,
componentsNumberARIMA,
lags, lagsModel, lagsModelAll, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, nonZeroARI, nonZeroMA, arEstimate, maEstimate, arimaPolynomials,
arOrders, iOrders, maOrders, arRequired, maRequired, armaParameters,
xregModel, xregNumber,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate,
arPolynomialMatrix, maPolynomialMatrix){
# Fill in the matrices
adamElements <- filler(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETS, componentsNumberETSNonSeasonal,
componentsNumberETSSeasonal, componentsNumberARIMA,
lags, lagsModel, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arEstimate, maEstimate, arOrders, iOrders, maOrders,
arRequired, maRequired, armaParameters,
nonZeroARI, nonZeroMA, arimaPolynomials,
xregModel, xregNumber);
# If we estimate lambda, take it from the B vector
if(lambdaEstimate){
# Take absolute value, just to be on safe side. We don't need negatives anyway.
lambda[] <- abs(B[length(B)]);
# Lambda is restricted by 0.25 if it is optimised.
if(any(distribution==c("dgnorm","dlgnorm")) && lambda<0.25){
return(1E+10/lambda);
}
}
# Check the bounds, classical restrictions
if(bounds=="usual"){
# Stationarity and invertibility conditions for ARIMA
if(arimaModel && any(c(arEstimate,maEstimate))){
# Calculate the polynomial roots for AR
if(arEstimate){
arPolynomialMatrix[,1] <- -adamElements$arimaPolynomials$arPolynomial[-1];
arPolyroots <- abs(eigen(arPolynomialMatrix, symmetric=TRUE, only.values=TRUE)$values);
if(any(arPolyroots>1)){
return(1E+100*max(arPolyroots));
}
}
# Calculate the polynomial roots of MA
if(maEstimate){
maPolynomialMatrix[,1] <- adamElements$arimaPolynomials$maPolynomial[-1];
maPolyroots <- abs(eigen(maPolynomialMatrix, symmetric=TRUE ,only.values=TRUE)$values);
if(any(maPolyroots>1)){
return(1E+100*max(abs(maPolyroots)));
}
}
}
# Smoothing parameters & phi restrictions in case of ETS
if(etsModel){
if(any(adamElements$vecG[1:componentsNumberETS]>1) || any(adamElements$vecG[1:componentsNumberETS]<0)){
return(1E+300);
}
if(modelIsTrendy){
if((adamElements$vecG[2]>adamElements$vecG[1])){
return(1E+300);
}
if(modelIsSeasonal && any(adamElements$vecG[componentsNumberETSNonSeasonal+c(1:componentsNumberETSSeasonal)]>
(1-adamElements$vecG[1]))){
return(1E+300);
}
}
else{
if(modelIsSeasonal && any(adamElements$vecG[componentsNumberETSNonSeasonal+c(1:componentsNumberETSSeasonal)]>
(1-adamElements$vecG[1]))){
return(1E+300);
}
}
# This is the restriction on the damping parameter
if(phiEstimate && (adamElements$matF[2,2]>1 || adamElements$matF[2,2]<0)){
return(1E+300);
}
}
# Smoothing parameters for the explanatory variables (0, 1) region
if(xregModel && xregDo=="adapt"){
if(any(adamElements$vecG[componentsNumberETS+componentsNumberARIMA+1:xregNumber]>1) ||
any(adamElements$vecG[componentsNumberETS+componentsNumberARIMA+1:xregNumber]<0)){
return(1E+100*max(abs(adamElements$vecG[componentsNumberETS+componentsNumberARIMA+1:xregNumber]-0.5)));
}
}
}
else if(bounds=="admissible"){
# Stationarity condition of ARIMA
if(arimaModel){
# Calculate the polynomial roots for AR
if(arEstimate){
arPolynomialMatrix[,1] <- -adamElements$arimaPolynomials$arPolynomial[-1];
arPolyroots <- abs(eigen(arPolynomialMatrix, symmetric=TRUE, only.values=TRUE)$values);
if(any(arPolyroots>1)){
return(1E+100*max(arPolyroots));
}
}
}
# Stability / invertibility condition for ETS / ARIMA.
if(etsModel || arimaModel){
if(xregModel){
# We check the condition on average
eigenValues <- abs(eigen((adamElements$matF -
matrix(adamElements$vecG, ncol(adamElements$matWt), obsInSample) %*%
adamElements$matWt[1:obsInSample,,drop=FALSE] / obsInSample),
symmetric=TRUE, only.values=TRUE)$values);
}
else{
eigenValues <- abs(eigen(adamElements$matF -
adamElements$vecG %*% adamElements$matWt[obsInSample,,drop=FALSE],
symmetric=TRUE, only.values=TRUE)$values);
}
if(any(eigenValues>1+1E-50)){
return(1E+100*max(eigenValues));
}
}
}
# Produce fitted values and errors
adamFitted <- adamFitterWrap(adamElements$matVt, adamElements$matWt, adamElements$matF, adamElements$vecG,
lagsModelAll, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, yInSample, ot, initialType=="backcasting");
if(!multisteps){
if(loss=="likelihood"){
# Scale for different functions
scale <- scaler(distribution, Etype, adamFitted$errors[otLogical],
adamFitted$yFitted[otLogical], obsInSample, lambda);
# Calculate the likelihood
## as.complex() is needed for failsafe in case of exotic models
CFValue <- -sum(switch(distribution,
"dnorm"=switch(Etype,
"A"=dnorm(x=yInSample[otLogical], mean=adamFitted$yFitted[otLogical],
sd=scale, log=TRUE),
"M"=dnorm(x=yInSample[otLogical], mean=adamFitted$yFitted[otLogical],
sd=scale*adamFitted$yFitted[otLogical], log=TRUE)),
"dlaplace"=switch(Etype,
"A"=dlaplace(q=yInSample[otLogical], mu=adamFitted$yFitted[otLogical],
scale=scale, log=TRUE),
"M"=dlaplace(q=yInSample[otLogical], mu=adamFitted$yFitted[otLogical],
scale=scale*adamFitted$yFitted[otLogical], log=TRUE)),
"ds"=switch(Etype,
"A"=ds(q=yInSample[otLogical],mu=adamFitted$yFitted[otLogical],
scale=scale, log=TRUE),
"M"=ds(q=yInSample[otLogical],mu=adamFitted$yFitted[otLogical],
scale=scale*sqrt(adamFitted$yFitted[otLogical]), log=TRUE)),
"dgnorm"=switch(Etype,
"A"=dgnorm(x=yInSample[otLogical],mu=adamFitted$yFitted[otLogical],
alpha=scale, beta=lambda, log=TRUE),
# Suppres Warnings is needed, because the check is done for scalar alpha
"M"=suppressWarnings(dgnorm(x=yInSample[otLogical],
mu=adamFitted$yFitted[otLogical],
alpha=scale*adamFitted$yFitted[otLogical],
beta=lambda, log=TRUE))),
"dlogis"=switch(Etype,
"A"=dlogis(x=yInSample[otLogical], location=adamFitted$yFitted[otLogical],
scale=scale, log=TRUE),
"M"=dlogis(x=yInSample[otLogical], location=adamFitted$yFitted[otLogical],
scale=scale*adamFitted$yFitted[otLogical], log=TRUE)),
"dt"=switch(Etype,
"A"=dt(adamFitted$errors[otLogical], df=abs(lambda), log=TRUE),
"M"=dt(adamFitted$errors[otLogical]*adamFitted$yFitted[otLogical],
df=abs(lambda), log=TRUE)),
"dalaplace"=switch(Etype,
"A"=dalaplace(q=yInSample[otLogical], mu=adamFitted$yFitted[otLogical],
scale=scale, alpha=lambda, log=TRUE),
"M"=dalaplace(q=yInSample[otLogical], mu=adamFitted$yFitted[otLogical],
scale=scale*adamFitted$yFitted[otLogical], alpha=lambda, log=TRUE)),
"dlnorm"=dlnorm(x=yInSample[otLogical], meanlog=Re(log(as.complex(adamFitted$yFitted[otLogical]))),
sdlog=scale, log=TRUE),
"dllaplace"=dlaplace(q=log(yInSample[otLogical]), mu=Re(log(as.complex(adamFitted$yFitted[otLogical]))),
scale=scale, log=TRUE) -log(yInSample[otLogical]),
"dls"=ds(q=log(yInSample[otLogical]), mu=Re(log(as.complex(adamFitted$yFitted[otLogical]))),
scale=scale, log=TRUE) -log(yInSample[otLogical]),
"dlgnorm"=dgnorm(x=log(yInSample[otLogical]),mu=Re(log(as.complex(adamFitted$yFitted[otLogical]))),
alpha=scale, beta=lambda, log=TRUE) -log(yInSample[otLogical]),
# "dinvgauss"=dinvgauss(x=1+adamFitted$errors, mean=1,
# dispersion=scale, log=TRUE)));
# "dinvgauss"=dinvgauss(x=yInSampleNew, mean=adamFitted$yFitted,
# dispersion=scale/adamFitted$yFitted, log=TRUE)));
"dinvgauss"=dinvgauss(x=yInSample[otLogical], mean=adamFitted$yFitted[otLogical],
dispersion=scale/adamFitted$yFitted[otLogical], log=TRUE)));
# Differential entropy for the logLik of occurrence model
if(occurrenceModel || any(!otLogical)){
CFValue <- CFValue + switch(distribution,
"dnorm" =,
"dlnorm" = obsZero*(log(sqrt(2*pi)*scale)+0.5),
"dlogis" = obsZero*2,
"dlaplace" =,
"dllaplace" =,
"dalaplace" = obsZero*(1 + log(2*scale)),
"ds" =,
"dls" = obsZero*(2 + 2*log(2*scale)),
"dgnorm" =,
"dlgnorm" = obsZero*(1/lambda-log(lambda/(2*scale*gamma(1/lambda)))),
"dt" = obsZero*((scale+1)/2 *
(digamma((scale+1)/2)-digamma(scale/2)) +
log(sqrt(scale) * beta(scale/2,0.5))),
# "dinvgauss" = obsZero*(0.5*(log(pi/2)+1+suppressWarnings(log(scale)))));
# "dinvgauss" =0);
"dinvgauss" = 0.5*(obsZero*(log(pi/2)+1+suppressWarnings(log(scale)))-
sum(log(adamFitted$yFitted[!otLogical]))));
}
}
else if(loss=="MSE"){
CFValue <- sum(adamFitted$errors^2)/obsInSample;
}
else if(loss=="MAE"){
CFValue <- sum(abs(adamFitted$errors))/obsInSample;
}
else if(loss=="HAM"){
CFValue <- sum(sqrt(abs(adamFitted$errors)))/obsInSample;
}
else if(any(loss==c("LASSO","RIDGE"))){
### All of this is needed in order to normalise level, trend, seasonal and xreg parameters
# Define, how many elements to skip (we don't normalise smoothing parameters)
if(persistenceXregEstimate){
persistenceToSkip <- componentsNumberETS+componentsNumberARIMA+xregNumber;
}
else{
persistenceToSkip <- componentsNumberETS+componentsNumberARIMA;
}
j <- 1;
if(phiEstimate){
j[] <- 2;
}
if(initialType=="optimal"){
# Standardise the level
B[persistenceToSkip+j] <- (B[persistenceToSkip+j] - mean(yInSample[1:lagsModelMax])) / mean(yInSample[1:lagsModelMax]);
# Change B values for the trend, so that it shrinks properly
if(Ttype=="M"){
j[] <- j+1;
B[persistenceToSkip+j] <- log(B[persistenceToSkip+j]);
}
else if(Ttype=="A"){
j[] <- j+1;
B[persistenceToSkip+j] <- B[persistenceToSkip+j]/mean(yInSample);
}
# Change B values for seasonality, so that it shrinks properly
if(Stype=="M"){
B[persistenceToSkip+j+1:(sum(lagsModel)-j)] <- log(B[persistenceToSkip+j+1:(sum(lagsModel)-j)]);
}
else if(Stype=="A"){
B[persistenceToSkip+j+1:(sum(lagsModel)-j)] <- B[persistenceToSkip+j+1:(sum(lagsModel)-j)]/mean(yInSample);
}
# Normalise parameters of xreg if they are additive. Otherwise leave - they will be small and close to zero
if(xregNumber>0 && Etype=="A"){
denominator <- tail(colMeans(matWt),xregNumber);
# If it is lower than 1, then we are probably dealing with (0, 1). No need to normalise
denominator[abs(denominator)<1] <- 1;
B[persistenceToSkip+sum(lagsModel)+c(1:xregNumber)] <- tail(B,xregNumber) / denominator;
}
}
CFValue <- (switch(Etype,
"A"=(1-lambda)* sqrt(sum(adamFitted$errors^2))/obsInSample,
"M"=(1-lambda)* sqrt(sum(log(1+adamFitted$errors)^2))/obsInSample) +
switch(loss,
"LASSO"=lambda * sum(abs(B)),
"RIDGE"=lambda * sqrt(sum((B)^2))));
}
else if(loss=="custom"){
CFValue <- lossFunction(actual=yInSample,fitted=adamFitted$yFitted,B=B);
}
}
else{
# Call for the Rcpp function to produce a matrix of multistep errors
adamErrors <- adamErrorerWrap(adamFitted$matVt, adamElements$matWt, adamElements$matF,
lagsModelAll, Etype, Ttype, Stype,
componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, h,
yInSample, ot);
# This is a fix for the multistep in case of Etype=="M", assuming logN
if(Etype=="M"){
adamErrors[] <- log(1+adamErrors);
}
# Not done yet: "aMSEh","aTMSE","aGTMSE","aMSCE","aGPL"
CFValue <- switch(loss,
"MSEh"=sum(adamErrors[,h]^2)/(obsInSample-h),
"TMSE"=sum(colSums(adamErrors^2)/(obsInSample-h)),
"GTMSE"=sum(log(colSums(adamErrors^2)/(obsInSample-h))),
"MSCE"=sum(rowSums(adamErrors)^2)/(obsInSample-h),
"MAEh"=sum(abs(adamErrors[,h]))/(obsInSample-h),
"TMAE"=sum(colSums(abs(adamErrors))/(obsInSample-h)),
"GTMAE"=sum(log(colSums(abs(adamErrors))/(obsInSample-h))),
"MACE"=sum(abs(rowSums(adamErrors)))/(obsInSample-h),
"HAMh"=sum(sqrt(abs(adamErrors[,h])))/(obsInSample-h),
"THAM"=sum(colSums(sqrt(abs(adamErrors)))/(obsInSample-h)),
"GTHAM"=sum(log(colSums(sqrt(abs(adamErrors)))/(obsInSample-h))),
"CHAM"=sum(sqrt(abs(rowSums(adamErrors))))/(obsInSample-h),
"GPL"=log(det(t(adamErrors) %*% adamErrors/(obsInSample-h))),
0);
}
if(is.na(CFValue) || is.nan(CFValue)){
CFValue[] <- 1e+300;
}
return(CFValue);
}
#### The function returns log-likelihood of the model ####
logLikADAM <- function(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, yInSample,
ot, otLogical, occurrenceModel, pFitted, obsInSample,
componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal,
componentsNumberARIMA,
lags, lagsModel, lagsModelAll, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, nonZeroARI, nonZeroMA, arEstimate, maEstimate, arimaPolynomials,
arOrders, iOrders, maOrders, arRequired, maRequired, armaParameters,
xregModel, xregNumber,
bounds, loss, lossFunction, distribution, horizon, multisteps, lambda, lambdaEstimate,
arPolynomialMatrix, maPolynomialMatrix){
if(!multisteps){
if(any(loss==c("LASSO","RIDGE"))){
return(0);
}
else{
distributionNew <- switch(loss,
"MSE"=switch(Etype,"A"="dnorm","M"="dlnorm"),
"MAE"=switch(Etype,"A"="dlaplace","M"="dllaplace"),
"HAM"=switch(Etype,"A"="ds","M"="dls"),
distribution);
lossNew <- switch(loss,
"MSE"=,"MAE"=,"HAM"="likelihood",
loss)
# bounds="none" switches off the checks of parameters.
logLikReturn <- -CF(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, yInSample,
ot, otLogical, occurrenceModel, obsInSample,
componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal,
componentsNumberARIMA,
lags, lagsModel, lagsModelAll, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, nonZeroARI, nonZeroMA, arEstimate, maEstimate, arimaPolynomials,
arOrders, iOrders, maOrders, arRequired, maRequired, armaParameters,
xregModel, xregNumber,
bounds="none", lossNew, lossFunction, distributionNew,
horizon, multisteps, lambda, lambdaEstimate,
arPolynomialMatrix, maPolynomialMatrix);
# If this is an occurrence model, add the probabilities
if(occurrenceModel){
if(is.infinite(logLikReturn)){
logLikReturn[] <- 0;
}
if(any(c(1-pFitted[!otLogical]==0,pFitted[otLogical]==0))){
# return(-Inf);
ptNew <- pFitted[(pFitted!=0) & (pFitted!=1)];
otNew <- ot[(pFitted!=0) & (pFitted!=1)];
# Just return the original likelihood if the probability is weird
if(length(ptNew)==0){
return(logLikReturn);
}
else{
return(logLikReturn + sum(log(ptNew[otNew==1])) + sum(log(1-ptNew[otNew==0])));
}
}
else{
return(logLikReturn + sum(log(pFitted[otLogical])) + sum(log(1-pFitted[!otLogical])));
}
}
else{
return(logLikReturn);
}
}
}
else{
# Use the predictive likelihoods from the GPL paper:
# - Normal for MSEh, MSCE, GPL and their analytical counterparts
# - Laplace for MAEh and MACE,
# - S for HAMh and CHAM
# bounds="none" switches off the checks of parameters.
logLikReturn <- CF(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, yInSample,
ot, otLogical, occurrenceModel, obsInSample,
componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal,
componentsNumberARIMA,
lags, lagsModel, lagsModelAll, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, nonZeroARI, nonZeroMA, arEstimate, maEstimate, arimaPolynomials,
arOrders, iOrders, maOrders, arRequired, maRequired, armaParameters,
xregModel, xregNumber,
bounds="none", loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate,
arPolynomialMatrix, maPolynomialMatrix);
logLikReturn[] <- -switch(loss,
"MSEh"=, "aMSEh"=, "TMSE"=, "aTMSE"=, "MSCE"=, "aMSCE"=
(obsInSample-h)/2*(log(2*pi)+1+log(logLikReturn)),
"GTMSE"=, "aGTMSE"=
(obsInSample-h)/2*(log(2*pi)+1+logLikReturn),
"MAEh"=, "TMAE"=, "GTMAE"=, "MACE"=
(obsInSample-h)*(log(2)+1+log(logLikReturn)),
"HAMh"=, "THAM"=, "GTHAM"=, "CHAM"=
(obsInSample-h)*(log(4)+2+2*log(logLikReturn)),
#### Divide GPL by 8 in order to make it comparable with the univariate ones
"GPL"=, "aGPL"=
(obsInSample-h)/2*(h*log(2*pi)+h+logLikReturn)/h);
# This is not well motivated at the moment, but should make likelihood comparable, taking T instead of T-h
logLikReturn[] <- logLikReturn / (obsInSample-h) * obsInSample;
# In case of multiplicative model, we assume log- distribution
if(Etype=="M"){
logLikReturn[] <- logLikReturn - sum(log(yInSample));
}
return(logLikReturn);
}
}
#### The function estimates the ETS model and returns B, logLik, nParam and CF(B) ####
estimator <- function(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate){
# Create the basic variables
adamArchitect <- architector(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal,
xregNumber, obsInSample, initialType,
arimaModel, lagsModelARIMA, xregModel);
list2env(adamArchitect, environment());
# Create the matrices for the specific ETS model
adamCreated <- creator(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
lags, lagsModel, lagsModelARIMA, lagsModelAll, lagsModelMax,
obsStates, obsInSample, obsAll, componentsNumberETS, componentsNumberETSSeasonal,
componentsNamesETS, otLogical, yInSample,
persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate, persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi,
initialType, initialEstimate,
initialLevel, initialLevelEstimate, initialTrend, initialTrendEstimate,
initialSeasonal, initialSeasonalEstimate,
initialArima, initialArimaEstimate, initialArimaNumber,
initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
arOrders, iOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames);
# If B is not provided, initialise it
if(is.null(B)){
BValues <- initialiser(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETSNonSeasonal, componentsNumberETSSeasonal, componentsNumberETS,
lags, lagsModel, lagsModelSeasonal, lagsModelARIMA, lagsModelMax,
adamCreated$matVt,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arRequired, maRequired, arEstimate, maEstimate, arOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA, initialArimaNumber,
xregModel, xregNumber, lambdaEstimate);
}
# print(BValues$B);
# Preheat the initial state of ARIMA. Do this only for optimal initials and if B is not provided
if(arimaModel && initialType=="optimal" && initialArimaEstimate && is.null(B)){
adamElements <- filler(BValues$B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETS, componentsNumberETSNonSeasonal,
componentsNumberETSSeasonal, componentsNumberARIMA,
lags, lagsModel, lagsModelMax,
adamCreated$matVt, adamCreated$matWt, adamCreated$matF, adamCreated$vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arEstimate, maEstimate, arOrders, iOrders, maOrders,
arRequired, maRequired, armaParameters,
nonZeroARI, nonZeroMA, adamCreated$arimaPolynomials,
xregModel, xregNumber);
# Do initial fit to get the state values from the backcasting
adamFitted <- adamFitterWrap(cbind(adamElements$matVt,adamElements$matVt[,lagsModelMax+1:lagsModelMax,drop=FALSE]),
adamElements$matWt, adamElements$matF, adamElements$vecG,
lagsModelAll, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, yInSample, ot, TRUE);
adamElements$matVt[,1:lagsModelMax] <- adamFitted$matVt[,1:lagsModelMax];
# Produce new initials
BValuesNew <- initialiser(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETSNonSeasonal, componentsNumberETSSeasonal, componentsNumberETS,
lags, lagsModel, lagsModelSeasonal, lagsModelARIMA, lagsModelMax,
adamElements$matVt,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arRequired, maRequired, arEstimate, maEstimate, arOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA, initialArimaNumber,
xregModel, xregNumber, lambdaEstimate);
B <- BValuesNew$B;
# Failsafe, just in case if the initial values contain NA / NaN
if(any(is.na(B))){
B[is.na(B)] <- BValues$B[is.na(B)];
}
if(any(is.nan(B))){
B[is.nan(B)] <- BValues$B[is.nan(B)];
}
}
# Create the vector of initials for the optimisation
if(is.null(B)){
B <- BValues$B
# lb <- BValues$Bl;
# ub <- BValues$Bu;
}
# Matrices needed for the polynomials calculation -> stationarity / stability checks
if(arimaModel){
# AR polynomials
arPolynomialMatrix <- matrix(0, arOrders %*% lags, arOrders %*% lags);
if(nrow(arPolynomialMatrix)>1){
arPolynomialMatrix[2:nrow(arPolynomialMatrix)-1,2:nrow(arPolynomialMatrix)] <- diag(nrow(arPolynomialMatrix)-1);
}
# MA polynomials
maPolynomialMatrix <- matrix(0, maOrders %*% lags, maOrders %*% lags);
if(nrow(maPolynomialMatrix)>1){
maPolynomialMatrix[2:nrow(maPolynomialMatrix)-1,2:nrow(maPolynomialMatrix)] <- diag(nrow(maPolynomialMatrix)-1);
}
}
else{
maPolynomialMatrix <- arPolynomialMatrix <- NULL;
}
# If the distribution is default, change it according to the error term
if(distribution=="default"){
distributionNew <- switch(Etype,
"A"=switch(loss,
"MAEh"=, "MACE"=, "MAE"="dlaplace",
"HAMh"=, "CHAM"=, "HAM"="ds",
"MSEh"=, "MSCE"=, "MSE"=, "GPL"=, "likelihood"=, "dnorm"),
"M"=switch(loss,
"MAEh"=, "MACE"=, "MAE"="dllaplace",
"HAMh"=, "CHAM"=, "HAM"="dls",
"MSEh"=, "MSCE"=, "MSE"=, "GPL"="dlnorm",
"likelihood"=, "dinvgauss"));
}
else{
distributionNew <- distribution;
}
# print(B)
# print(Etype)
# print(Ttype)
# print(Stype)
# stop()
print_level_hidden <- print_level;
if(print_level==41){
print_level[] <- 0;
}
# Parameters are chosen to speed up the optimisation process and have decent accuracy
res <- suppressWarnings(nloptr(B, CF, lb=lb, ub=ub,
opts=list(algorithm=algorithm, xtol_rel=xtol_rel, xtol_abs=xtol_abs,
ftol_rel=ftol_rel, ftol_abs=ftol_abs,
maxeval=maxeval, maxtime=maxtime, print_level=print_level),
etsModel=etsModel, Etype=Etype, Ttype=Ttype, Stype=Stype, modelIsTrendy=modelIsTrendy,
modelIsSeasonal=modelIsSeasonal, yInSample=yInSample,
ot=ot, otLogical=otLogical, occurrenceModel=occurrenceModel, obsInSample=obsInSample,
componentsNumberETS=componentsNumberETS, componentsNumberETSSeasonal=componentsNumberETSSeasonal,
componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal,
componentsNumberARIMA=componentsNumberARIMA,
lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax,
matVt=adamCreated$matVt, matWt=adamCreated$matWt, matF=adamCreated$matF, vecG=adamCreated$vecG,
persistenceEstimate=persistenceEstimate, persistenceLevelEstimate=persistenceLevelEstimate,
persistenceTrendEstimate=persistenceTrendEstimate,
persistenceSeasonalEstimate=persistenceSeasonalEstimate,
persistenceXregEstimate=persistenceXregEstimate,
phiEstimate=phiEstimate, initialType=initialType,
initialEstimate=initialEstimate, initialLevelEstimate=initialLevelEstimate,
initialTrendEstimate=initialTrendEstimate, initialSeasonalEstimate=initialSeasonalEstimate,
initialArimaEstimate=initialArimaEstimate, initialXregEstimate=initialXregEstimate,
arimaModel=arimaModel, nonZeroARI=nonZeroARI, nonZeroMA=nonZeroMA,
arimaPolynomials=adamCreated$arimaPolynomials,
arEstimate=arEstimate, maEstimate=maEstimate,
arOrders=arOrders, iOrders=iOrders, maOrders=maOrders,
arRequired=arRequired, maRequired=maRequired, armaParameters=armaParameters,
xregModel=xregModel, xregNumber=xregNumber,
bounds=bounds, loss=loss, lossFunction=lossFunction, distribution=distributionNew,
horizon=horizon, multisteps=multisteps, lambda=lambda, lambdaEstimate=lambdaEstimate,
arPolynomialMatrix=arPolynomialMatrix, maPolynomialMatrix=maPolynomialMatrix));
if(is.infinite(res$objective) || res$objective==1e+300){
# If the optimisation didn't work, give it another try with zero initials for smoothing parameters
if(etsModel){
B[1:sum(persistenceLevelEstimate,persistenceTrendEstimate,persistenceSeasonalEstimate)] <- 0;
}
if(arimaModel){
B[sum(persistenceLevelEstimate,persistenceTrendEstimate,persistenceSeasonalEstimate,
persistenceXregEstimate*xregNumber)+c(1:sum(arOrders*arEstimate,maOrders*maEstimate))] <- 0.01;
}
# print(B)
res <- suppressWarnings(nloptr(B, CF, lb=lb, ub=ub,
opts=list(algorithm="NLOPT_LN_SBPLX", xtol_rel=xtol_rel,
ftol_rel=ftol_rel, ftol_abs=ftol_abs,
maxeval=maxeval, maxtime=maxtime, print_level=print_level),
etsModel=etsModel, Etype=Etype, Ttype=Ttype, Stype=Stype, modelIsTrendy=modelIsTrendy,
modelIsSeasonal=modelIsSeasonal, yInSample=yInSample,
ot=ot, otLogical=otLogical, occurrenceModel=occurrenceModel, obsInSample=obsInSample,
componentsNumberETS=componentsNumberETS, componentsNumberETSSeasonal=componentsNumberETSSeasonal,
componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal,
componentsNumberARIMA=componentsNumberARIMA,
lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax,
matVt=adamCreated$matVt, matWt=adamCreated$matWt, matF=adamCreated$matF, vecG=adamCreated$vecG,
persistenceEstimate=persistenceEstimate,
persistenceLevelEstimate=persistenceLevelEstimate,
persistenceTrendEstimate=persistenceTrendEstimate,
persistenceSeasonalEstimate=persistenceSeasonalEstimate,
persistenceXregEstimate=persistenceXregEstimate,
phiEstimate=phiEstimate, initialType=initialType,
initialEstimate=initialEstimate, initialLevelEstimate=initialLevelEstimate,
initialTrendEstimate=initialTrendEstimate, initialSeasonalEstimate=initialSeasonalEstimate,
initialArimaEstimate=initialArimaEstimate, initialXregEstimate=initialXregEstimate,
arimaModel=arimaModel, nonZeroARI=nonZeroARI, nonZeroMA=nonZeroMA,
arimaPolynomials=adamCreated$arimaPolynomials,
arEstimate=arEstimate, maEstimate=maEstimate,
arOrders=arOrders, iOrders=iOrders, maOrders=maOrders,
arRequired=arRequired, maRequired=maRequired, armaParameters=armaParameters,
xregModel=xregModel, xregNumber=xregNumber,
bounds=bounds, loss=loss, lossFunction=lossFunction, distribution=distributionNew,
horizon=horizon, multisteps=multisteps, lambda=lambda, lambdaEstimate=lambdaEstimate,
arPolynomialMatrix=arPolynomialMatrix, maPolynomialMatrix=maPolynomialMatrix));
}
if(print_level_hidden>0){
print(res);
}
##### !!! Check the obtained parameters and the loss value and remove redundant parameters !!! #####
# Cases to consider:
# 1. Some smoothing parameters are zero or one;
# 2. The cost function value is -Inf (due to no variability in the sample);
# Prepare the values to return
B[] <- res$solution;
CFValue <- res$objective;
# In case of likelihood, we typically have one more parameter to estimate - scale
nParamEstimated <- length(B) + (loss=="likelihood");
# Return a proper logLik class
logLikADAMValue <- structure(logLikADAM(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, yInSample,
ot, otLogical, occurrenceModel, pFitted, obsInSample,
componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal,
componentsNumberARIMA,
lags, lagsModel, lagsModelAll, lagsModelMax,
adamCreated$matVt, adamCreated$matWt, adamCreated$matF, adamCreated$vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, nonZeroARI, nonZeroMA, arEstimate, maEstimate,
adamCreated$arimaPolynomials,
arOrders, iOrders, maOrders, arRequired, maRequired, armaParameters,
xregModel, xregNumber,
bounds, loss, lossFunction, distributionNew, horizon, multisteps,
lambda, lambdaEstimate, arPolynomialMatrix, maPolynomialMatrix),
nobs=obsInSample,df=nParamEstimated,class="logLik");
#### If we do variables selection, do it here, then reestimate the model. ####
if(xregDo=="select"){
# Fill in the matrices
adamElements <- filler(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETS, componentsNumberETSNonSeasonal,
componentsNumberETSSeasonal, componentsNumberARIMA,
lags, lagsModel, lagsModelMax,
adamCreated$matVt, adamCreated$matWt, adamCreated$matF, adamCreated$vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arEstimate, maEstimate, arOrders, iOrders, maOrders,
arRequired, maRequired, armaParameters,
nonZeroARI, nonZeroMA, adamCreated$arimaPolynomials,
xregModel, xregNumber);
# Fit the model to the data
adamFitted <- adamFitterWrap(adamElements$matVt, adamElements$matWt, adamElements$matF, adamElements$vecG,
lagsModelAll, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, yInSample, ot, initialType=="backcasting");
# Extract the errors corrrectly
errors <- switch(distribution,
"dlnorm"=, "dllaplace"=, "dls"=, "dlgnorm"=, "dinvgauss"=switch(Etype,
"A"=1+adamFitted$errors/adamFitted$yFitted,
"M"=adamFitted$errors),
"dnorm"=, "dlaplace"=, "ds"=, "dgnorm"=, "dlogis"=, "dt"=, "dalaplace"=,adamFitted$errors);
# Extract the errors and amend them to correspond to the distribution
errors[] <- errors + switch(Etype,"A"=0,"M"=1);
if(any(distribution==c("dlnorm","dllaplace","dls","dlgnorm")) && Etype=="A" && any(errors<0)){
errors[errors<0] <- 1e-100;
}
# Call the xregSelector providing the original matrix with the data
xregIndex <- switch(Etype,"A"=1,"M"=2);
xregModelInitials[[xregIndex]] <- xregSelector(errors=errors, xregData=xregDataOriginal, ic=ic,
df=length(B)+1, distribution=distributionNew, occurrence=oesModel);
xregNumber <- length(xregModelInitials[[xregIndex]]$initialXreg);
xregNames <- names(xregModelInitials[[xregIndex]]$initialXreg);
# Make sure that the accidental "`" don't reappear...
xregNames[] <- gsub("\`","",xregNames,ignore.case=TRUE);
# If there are some variables, then do the proper reestimation and return the new values
if(xregNumber>0){
xregModel[] <- TRUE;
initialXregEstimate[] <- initialXregEstimateOriginal;
persistenceXregEstimate[] <- persistenceXregEstimateOriginal;
xregData <- xregDataOriginal[,xregNames,drop=FALSE];
# Redefine loss for ALM
lossNew <- switch(loss,
"MSEh"=,"TMSE"=,"GTMSE"=,"MSCE"="MSE",
"MAEh"=,"TMAE"=,"GTMAE"=,"MACE"="MAE",
"HAMh"=,"THAM"=,"GTHAM"=,"CHAM"="HAM",
loss);
if(lossNew=="custom"){
lossNew <- lossFunction;
}
# Estimate alm again in order to get proper initials
almModel <- do.call(alm,list(formula=as.formula("y~."),
data=matrix(c(yInSample,xregData[1:obsInSample,,drop=FALSE]),
obsInSample,xregNumber+1,
dimnames=list(NULL,c("y",xregNames))),
distribution=distributionNew, loss=lossNew, occurrence=oesModel));
xregModelInitials[[xregIndex]]$initialXreg <- coef(almModel)[-1];
return(estimator(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo="use",
ot, otLogical, occurrenceModel, pFitted,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate));
}
}
return(list(B=B, CFValue=CFValue, nParamEstimated=nParamEstimated, logLikADAMValue=logLikADAMValue,
xregModel=xregModel, xregData=xregData, xregNumber=xregNumber, xregNames=xregNames, xregModelInitials=xregModelInitials,
initialXregEstimate=initialXregEstimate, persistenceXregEstimate=persistenceXregEstimate,
arimaPolynomials=adamCreated$arimaPolynomials));
}
#### The function creates a pool of models and selects the best of them ####
selector <- function(model, modelsPool, allowMultiplicative,
etsModel, Etype, Ttype, Stype, damped, lags,
lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted, ICFunction,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate){
# Check if the pool was provided. In case of "no", form the big and the small ones
if(is.null(modelsPool)){
# The variable saying that the pool was not provided.
if(!silent){
cat("Forming the pool of models based on... ");
}
# Define the whole pool of errors
if(!allowMultiplicative){
poolErrors <- c("A");
poolTrends <- c("N","A","Ad");
poolSeasonals <- c("N","A");
}
else{
poolErrors <- c("A","M");
poolTrends <- c("N","A","Ad","M","Md");
poolSeasonals <- c("N","A","M");
}
# Some preparation variables
# If Etype is not Z, then check on additive errors
if(Etype!="Z"){
poolErrors <- poolErrorsSmall <- Etype;
}
else{
poolErrorsSmall <- "A";
}
# If Ttype is not Z, then create a pool with specified type
if(Ttype!="Z"){
if(Ttype=="X"){
poolTrendsSmall <- c("N","A");
poolTrends <- c("N","A","Ad");
checkTrend <- TRUE;
}
else if(Ttype=="Y"){
poolTrendsSmall <- c("N","M");
poolTrends <- c("N","M","Md");
checkTrend <- TRUE;
}
else{
if(damped){
poolTrends <- poolTrendsSmall <- paste0(Ttype,"d");
}
else{
poolTrends <- poolTrendsSmall <- Ttype;
}
checkTrend <- FALSE;
}
}
else{
poolTrendsSmall <- c("N","A");
checkTrend <- TRUE;
}
# If Stype is not Z, then crete specific pools
if(Stype!="Z"){
if(Stype=="X"){
poolSeasonals <- poolSeasonalsSmall <- c("N","A");
checkSeasonal <- TRUE;
}
else if(Stype=="Y"){
poolSeasonalsSmall <- c("N","M");
poolSeasonals <- c("N","M");
checkSeasonal <- TRUE;
}
else{
poolSeasonalsSmall <- Stype;
poolSeasonals <- Stype;
checkSeasonal <- FALSE;
}
}
else{
poolSeasonalsSmall <- c("N","A","M");
checkSeasonal <- TRUE;
}
# If ZZZ, then the vector is: "ANN" "ANA" "ANM" "AAN" "AAA" "AAM"
# Otherwise id depends on the provided restrictions
poolSmall <- paste0(rep(poolErrorsSmall,length(poolTrendsSmall)*length(poolSeasonalsSmall)),
rep(poolTrendsSmall,each=length(poolSeasonalsSmall)),
rep(poolSeasonalsSmall,length(poolTrendsSmall)));
# Align error and seasonality, if the error was not forced to be additive
if(any(substr(poolSmall,3,3)=="M") && all(Etype!=c("A","X"))){
multiplicativeSeason <- (substr(poolSmall,3,3)=="M");
poolSmall[multiplicativeSeason] <- paste0("M",substr(poolSmall[multiplicativeSeason],2,3));
}
modelsTested <- NULL;
modelCurrent <- NA;
# Counter + checks for the components
j <- 1;
i <- 0;
check <- TRUE;
besti <- bestj <- 1;
results <- vector("list",length(poolSmall));
#### Branch and bound is here ####
while(check){
i <- i + 1;
modelCurrent[] <- poolSmall[j];
if(!silent){
cat(paste0(modelCurrent,", "));
}
Etype[] <- substring(modelCurrent,1,1);
Ttype[] <- substring(modelCurrent,2,2);
if(nchar(modelCurrent)==4){
phi[] <- 0.95;
phiEstimate[] <- TRUE;
Stype[] <- substring(modelCurrent,4,4);
}
else{
phi[] <- 1;
phiEstimate[] <- FALSE;
Stype[] <- substring(modelCurrent,3,3);
}
results[[i]] <- estimator(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate);
results[[i]]$IC <- ICFunction(results[[i]]$logLikADAMValue);
results[[i]]$Etype <- Etype;
results[[i]]$Ttype <- Ttype;
results[[i]]$Stype <- Stype;
results[[i]]$phiEstimate <- phiEstimate;
if(phiEstimate){
results[[i]]$phi <- results[[i]]$B[names(results[[i]]$B)=="phi"];
}
else{
results[[i]]$phi <- 1;
}
results[[i]]$model <- modelCurrent;
modelsTested <- c(modelsTested,modelCurrent);
if(j>1){
# If the first is better than the second, then choose first
if(results[[besti]]$IC <= results[[i]]$IC){
# If Ttype is the same, then we check seasonality
if(substring(modelCurrent,2,2)==substring(poolSmall[bestj],2,2)){
poolSeasonals <- results[[besti]]$Stype;
checkSeasonal <- FALSE;
# j[] <- j+1;
j <- which(poolSmall!=poolSmall[bestj] &
substring(poolSmall,nchar(poolSmall),nchar(poolSmall))==poolSeasonals);
}
# Otherwise we checked trend
else{
poolTrends <- results[[bestj]]$Ttype;
checkTrend[] <- FALSE;
}
}
else{
if(substring(modelCurrent,2,2) == substring(poolSmall[besti],2,2)){
poolSeasonals <- poolSeasonals[poolSeasonals!=results[[besti]]$Stype];
if(length(poolSeasonals)>1){
# Select another seasonal model, that is not from the previous iteration and not the current one
bestj[] <- j;
besti[] <- i;
# j[] <- 3;
j <- 3;
}
else{
bestj[] <- j;
besti[] <- i;
j <- which(substring(poolSmall,nchar(poolSmall),nchar(poolSmall))==poolSeasonals &
substring(poolSmall,2,2)!=substring(modelCurrent,2,2));
checkSeasonal[] <- FALSE;
}
}
else{
poolTrends <- poolTrends[poolTrends!=results[[bestj]]$Ttype];
besti[] <- i;
bestj[] <- j;
checkTrend[] <- FALSE;
}
}
if(all(!c(checkTrend,checkSeasonal))){
check[] <- FALSE;
}
}
else{
j <- 2;
}
if(j>=length(poolSmall)){
check[] <- FALSE;
}
}
# Prepare a bigger pool based on the small one
modelsPool <- unique(c(modelsTested,
paste0(rep(poolErrors,each=length(poolTrends)*length(poolSeasonals)),
poolTrends,
rep(poolSeasonals,each=length(poolTrends)))));
j <- length(modelsTested);
}
else{
j <- 0;
results <- vector("list",length(modelsPool));
}
modelsNumber <- length(modelsPool);
#### Run the full pool of models ####
if(!silent){
cat("Estimation progress: ");
}
# Start loop of models
while(j < modelsNumber){
j <- j + 1;
if(!silent){
if(j==1){
cat("\b");
}
cat(paste0(rep("\b",nchar(round((j-1)/modelsNumber,2)*100)+1),collapse=""));
cat(paste0(round(j/modelsNumber,2)*100,"%"));
}
modelCurrent <- modelsPool[j];
# print(modelCurrent)
Etype <- substring(modelCurrent,1,1);
Ttype <- substring(modelCurrent,2,2);
if(nchar(modelCurrent)==4){
phi[] <- 0.95;
Stype <- substring(modelCurrent,4,4);
phiEstimate <- TRUE;
}
else{
phi[] <- 1;
Stype <- substring(modelCurrent,3,3);
phiEstimate <- FALSE;
}
results[[j]] <- estimator(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate);
results[[j]]$IC <- ICFunction(results[[j]]$logLikADAMValue);
results[[j]]$Etype <- Etype;
results[[j]]$Ttype <- Ttype;
results[[j]]$Stype <- Stype;
results[[j]]$phiEstimate <- phiEstimate;
if(phiEstimate){
results[[j]]$phi <- results[[j]]$B[names(results[[j]]$B)=="phi"];
}
else{
results[[j]]$phi <- 1;
}
results[[j]]$model <- modelCurrent;
}
if(!silent){
cat("... Done! \n");
}
# Extract ICs and find the best
icSelection <- vector("numeric",modelsNumber);
for(i in 1:modelsNumber){
icSelection[i] <- results[[i]]$IC;
}
names(icSelection) <- modelsPool;
icSelection[is.nan(icSelection)] <- 1E100;
return(list(results=results,icSelection=icSelection));
}
##### Function uses residuals in order to determine the needed xreg #####
xregSelector <- function(errors, xregData, ic, df, distribution, occurrence){
stepwiseModel <- suppressWarnings(stepwise(cbind(as.data.frame(errors),xregData[1:obsInSample,,drop=FALSE]),
ic=ic, df=df, distribution=distribution, occurrence=occurrence, silent=TRUE));
return(list(initialXreg=coef(stepwiseModel)[-1],other=stepwiseModel$other,formula=formula(stepwiseModel)));
}
##### Function prepares all the matrices and vectors for return #####
preparator <- function(B, etsModel, Etype, Ttype, Stype,
lagsModel, lagsModelMax, lagsModelAll,
componentsNumberETS, componentsNumberETSSeasonal,
xregNumber, distribution, loss,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, lambdaEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
matVt, matWt, matF, vecG,
occurrenceModel, ot, oesModel,
parametersNumber, CFValue,
arimaModel, arRequired, maRequired,
arEstimate, maEstimate, arOrders, iOrders, maOrders,
nonZeroARI, nonZeroMA,
arimaPolynomials, armaParameters){
if(modelDo!="use"){
# Fill in the matrices
adamElements <- filler(B,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETS, componentsNumberETSNonSeasonal,
componentsNumberETSSeasonal, componentsNumberARIMA,
lags, lagsModel, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, arEstimate, maEstimate, arOrders, iOrders, maOrders,
arRequired, maRequired, armaParameters,
nonZeroARI, nonZeroMA, arimaPolynomials,
xregModel, xregNumber);
list2env(adamElements, environment());
}
# Write down lambda. It is always positive, so take abs
if(lambdaEstimate){
lambda[] <- abs(tail(B,1));
}
# Write down phi
if(phiEstimate){
phi[] <- B[names(B)=="phi"];
}
# Fit the model to the data
adamFitted <- adamFitterWrap(matVt, matWt, matF, vecG,
lagsModelAll, Etype, Ttype, Stype, componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, yInSample, ot, initialType=="backcasting");
if(any(yClasses=="ts")){
yFitted <- ts(rep(NA,obsInSample), start=yStart, frequency=yFrequency);
errors <- ts(rep(NA,obsInSample), start=yStart, frequency=yFrequency);
}
else{
yFitted <- zoo(rep(NA,obsInSample), order.by=yInSampleIndex);
errors <- zoo(rep(NA,obsInSample), order.by=yInSampleIndex);
}
errors[] <- adamFitted$errors;
yFitted[] <- adamFitted$yFitted;
# Check what was returned in the end
if(any(is.nan(yFitted)) || any(is.na(yFitted))){
warning("Something went wrong in the estimation of the model ETS(",
Etype,Ttype,ifelse(damped,"d",""),Stype,") and NaNs were produced. ",
"If this is a mixed model, consider using the pure ones instead.",
call.=FALSE);
}
if(occurrenceModel){
yFitted[] <- yFitted * pFitted;
}
matVt[] <- adamFitted$matVt;
if(initialType=="backcasting"){
matVt <- matVt[,1:(obsInSample+lagsModelMax), drop=FALSE];
}
# Produce forecasts if the horizon is non-zero
if(horizon>0){
if(any(yClasses=="ts")){
yForecast <- ts(rep(NA, horizon), start=yForecastStart, frequency=yFrequency);
}
else{
yForecast <- zoo(rep(NA, horizon), order.by=yForecastIndex);
}
yForecast[] <- adamForecasterWrap(matVt[,obsInSample+(1:lagsModelMax),drop=FALSE], tail(matWt,horizon), matF,
lagsModelAll, Etype, Ttype, Stype,
componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber,
horizon);
#### Make safety checks
# If there are NaN values
if(any(is.nan(yForecast))){
yForecast[is.nan(yForecast)] <- 0;
}
# Amend forecasts, multiplying by probability
if(occurrenceModel && !occurrenceModelProvided){
yForecast[] <- yForecast * c(forecast(oesModel, h=h)$mean);
}
else if(occurrenceModel && occurrenceModelProvided){
yForecast[] <- yForecast * pForecast;
}
}
else{
if(any(yClasses=="ts")){
yForecast <- ts(NA, start=yForecastStart, frequency=yFrequency);
}
else{
yForecast <- zoo(rep(NA, horizon), order.by=yForecastIndex);
}
}
# If the distribution is default, change it according to the error term
if(distribution=="default"){
distribution[] <- switch(Etype,
"A"=switch(loss,
"MAEh"=, "MACE"=, "MAE"="dlaplace",
"HAMh"=, "CHAM"=, "HAM"="ds",
"MSEh"=, "MSCE"=, "GPL"=, "MSE"=,
"aMSEh"=, "aMSCE"=, "aGPL"=, "likelihood"=, "dnorm"),
"M"="dinvgauss");
if(multisteps && Etype=="M"){
distribution[] <- switch(loss,
"MAEh"=, "MACE"=, "MAE"="dllaplace",
"HAMh"=, "CHAM"=, "HAM"="dls",
"MSEh"=, "MSCE"=, "GPL"=, "MSE"=,
"aMSEh"=, "aMSCE"=, "aGPL"=, "dlnorm");
}
}
#### Initial values to return ####
initialValue <- vector("list", etsModel*(1+modelIsTrendy+modelIsSeasonal)+arimaModel+xregModel);
initialValueETS <- vector("list", etsModel*length(lagsModel));
initialValueNames <- vector("character", etsModel*(1+modelIsTrendy+modelIsSeasonal)+arimaModel+xregModel);
# The vector that defines what was estimated in the model
initialEstimated <- vector("logical", etsModel*(1+modelIsTrendy+modelIsSeasonal*componentsNumberETSSeasonal)+
arimaModel+xregModel);
# Write down the initials of ETS
j <- 0;
if(etsModel){
# Write down level, trend and seasonal
for(i in 1:length(lagsModel)){
initialValueETS[[i]] <- tail(matVt[i,1:lagsModelMax],lagsModel[i]);
}
j[] <- j+1;
# Write down level in the final list
initialEstimated[j] <- initialLevelEstimate;
initialValue[[j]] <- initialValueETS[[j]];
initialValueNames[j] <- c("level");
names(initialEstimated)[j] <- initialValueNames[j];
if(modelIsTrendy){
j[] <- 2;
initialEstimated[j] <- initialTrendEstimate;
# Write down trend in the final list
initialValue[[j]] <- initialValueETS[[j]];
# Remove the trend from ETS list
initialValueETS[[j]] <- NULL;
initialValueNames[j] <- c("trend");
names(initialEstimated)[j] <- initialValueNames[j];
}
# Write down the initial seasonals
if(modelIsSeasonal){
initialEstimated[j+c(1:componentsNumberETSSeasonal)] <- initialSeasonalEstimate;
# Remove the level from ETS list
initialValueETS[[1]] <- NULL;
j[] <- j+1;
if(length(initialSeasonalEstimate)>1){
initialValue[[j]] <- initialValueETS;
initialValueNames[[j]] <- "seasonal";
names(initialEstimated)[j+0:(componentsNumberETSSeasonal-1)] <-
paste0(initialValueNames[j],c(1:componentsNumberETSSeasonal));
}
else{
initialValue[[j]] <- initialValueETS[[1]];
initialValueNames[[j]] <- "seasonal";
names(initialEstimated)[j] <- initialValueNames[j];
}
}
}
# Write down the ARIMA initials
if(arimaModel){
j[] <- j+1;
initialEstimated[j] <- initialArimaEstimate;
if(initialArimaEstimate && initialType=="optimal"){
initialValue[[j]] <- B[substr(names(B),1,10)=="ARIMAState"];
}
else if(initialArimaEstimate && initialType=="backcasting"){
initialValue[[j]] <- head(matVt[componentsNumberETS+componentsNumberARIMA,],initialArimaNumber);
}
else{
initialValue[[j]] <- initialArima;
}
initialValueNames[j] <- "arima";
names(initialEstimated)[j] <- initialValueNames[j];
}
# Write down the xreg initials
if(xregModel){
j[] <- j+1;
initialEstimated[j] <- initialXregEstimate;
initialValue[[j]] <- matVt[componentsNumberETS+componentsNumberARIMA+1:xregNumber,lagsModelMax];
initialValueNames[j] <- "xreg";
names(initialEstimated)[j] <- initialValueNames[j];
}
names(initialValue) <- initialValueNames;
#### Persistence to return ####
persistence <- as.vector(vecG);
names(persistence) <- rownames(vecG);
# Remove xreg persistence from the returned vector
if(xregModel && xregDo!="adapt"){
persistence <- persistence[substr(names(persistence),1,5)!="delta"];
# We've selected the variables, so there's nothing to select anymore
xregDo <- "use";
}
else if(!xregModel){
xregDo <- NULL;
}
if(arimaModel){
armaParametersList <- vector("list",arRequired+maRequired);
j[] <- 1;
if(arRequired && arEstimate){
# Avoid damping parameter phi
armaParametersList[[j]] <- B[nchar(names(B))>3 & substr(names(B),1,3)=="phi"];
names(armaParametersList)[j] <- "ar";
j[] <- j+1;
}
# If this was provided
else if(arRequired && !arEstimate){
# Avoid damping parameter phi
armaParametersList[[j]] <- armaParameters[substr(names(armaParameters),1,3)=="phi"];
names(armaParametersList)[j] <- "ar";
j[] <- j+1;
}
if(maRequired && maEstimate){
armaParametersList[[j]] <- B[substr(names(B),1,5)=="theta"];
names(armaParametersList)[j] <- "ma";
}
else if(maRequired && !maEstimate){
armaParametersList[[j]] <- armaParameters[substr(names(armaParameters),1,5)=="theta"];
names(armaParametersList)[j] <- "ma";
}
}
else{
armaParametersList <- NULL;
}
# which() is needed in order to overcome weird behaviour of zoo
scale <- scaler(distribution, Etype, errors[which(otLogical)], yFitted[which(otLogical)], obsInSample, lambda);
# Amend the class of state matrix
if(any(yClasses=="ts")){
matVt <- ts(t(matVt), start=(time(y)[1]-deltat(y)*lagsModelMax), frequency=yFrequency);
}
else{
yStatesIndex <- yInSampleIndex[1] - lagsModelMax*diff(tail(yInSampleIndex,2)) +
c(1:lagsModelMax-1)*diff(tail(yInSampleIndex,2));
yStatesIndex <- c(yStatesIndex, yInSampleIndex);
matVt <- zoo(t(matVt), order.by=yStatesIndex);
}
parametersNumber[2,4] <- sum(parametersNumber[2,1:3]);
return(list(model=NA, timeElapsed=NA,
y=NA, holdout=NA, fitted=yFitted, residuals=errors,
forecast=yForecast, states=matVt,
persistence=persistence, phi=phi, transition=matF,
measurement=matWt, initial=initialValue, initialType=initialType,
initialEstimated=initialEstimated, orders=orders, arma=armaParametersList,
nParam=parametersNumber, occurrence=oesModel, xreg=xregData,
formula=formula, xregDo=xregDo,
loss=loss, lossValue=CFValue, logLik=logLikADAMValue, distribution=distribution,
scale=scale, lambda=lambda, B=B, lags=lags, lagsAll=lagsModelAll, FI=FI));
}
#### Deal with occurrence model ####
if(occurrenceModel && !occurrenceModelProvided){
modelForOES <- model;
if(model=="NNN"){
modelForOES[] <- "MNN";
}
oesModel <- suppressWarnings(oes(ot, model=modelForOES, occurrence=occurrence, ic=ic, h=horizon,
holdout=FALSE, bounds="usual", xreg=xregData, xregDo=xregDo, silent=TRUE));
pFitted[] <- fitted(oesModel);
parametersNumber[1,3] <- nparam(oesModel);
# print(oesModel)
# This should not happen, but just in case...
if(oesModel$occurrence=="n"){
occurrence <- "n";
otLogical <- rep(TRUE,obsInSample);
occurrenceModel <- FALSE;
ot <- matrix(otLogical*1,ncol=1);
obsNonzero <- sum(ot);
obsZero <- obsInSample - obsNonzero;
Etype[] <- switch(Etype,
"M"="A",
"Y"=,
"Z"="X",
Etype);
Ttype[] <- switch(Ttype,
"M"="A",
"Y"=,
"Z"="X",
Ttype);
Stype[] <- switch(Stype,
"M"="A",
"Y"=,
"Z"="X",
Stype);
}
}
else if(occurrenceModel && occurrenceModelProvided){
parametersNumber[2,3] <- nparam(oesModel);
}
##### Prepare stuff for the variables selection if xregDo="select" #####
if(xregDo=="select"){
# First, record the original parameters
xregExistOriginal <- xregModel;
initialXregsProvidedOriginal <- initialXregProvided;
initialXregEstimateOriginal <- initialXregEstimate;
persistenceXregOriginal <- persistenceXreg;
persistenceXregProvidedOriginal <- persistenceXregProvided;
persistenceXregEstimateOriginal <- persistenceXregEstimate;
xregModelOriginal <- xregModelInitials;
xregDataOriginal <- xregData;
xregNumberOriginal <- xregNumber;
xregNamesOriginal <- xregNames;
# Set the parameters to zero and do simple ETS
xregModel[] <- FALSE;
initialXregProvided <- FALSE;
initialXregEstimate[] <- FALSE;
persistenceXreg <- 0;
persistenceXregProvided <- FALSE;
persistenceXregEstimate[] <- FALSE;
xregModelInitials[[1]] <- NULL;
xregModelInitials[[2]] <- NULL;
xregData <- NULL;
xregNumber[] <- 0;
xregNames <- NULL;
}
##### Estimate the specified model #####
if(modelDo=="estimate"){
# Estimate the parameters of the demand sizes model
adamEstimated <- estimator(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate);
list2env(adamEstimated, environment());
#### This part is needed in order for the filler to do its job later on
# Create the basic variables based on the estimated model
adamArchitect <- architector(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal,
xregNumber, obsInSample, initialType,
arimaModel, lagsModelARIMA, xregModel);
list2env(adamArchitect, environment());
# Create the matrices for the specific ETS model
adamCreated <- creator(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
lags, lagsModel, lagsModelARIMA, lagsModelAll, lagsModelMax,
obsStates, obsInSample, obsAll, componentsNumberETS, componentsNumberETSSeasonal,
componentsNamesETS, otLogical, yInSample,
persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate, persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi,
initialType, initialEstimate,
initialLevel, initialLevelEstimate, initialTrend, initialTrendEstimate,
initialSeasonal, initialSeasonalEstimate,
initialArima, initialArimaEstimate, initialArimaNumber,
initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
arOrders, iOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames);
list2env(adamCreated, environment());
icSelection <- ICFunction(adamEstimated$logLikADAMValue);
####!!! If the occurrence is auto, then compare this with the model with no occurrence !!!####
parametersNumber[1,1] <- nParamEstimated;
if(xregModel){
parametersNumber[1,2] <- xregNumber*initialXregEstimate + xregNumber*persistenceXregEstimate;
}
parametersNumber[1,4] <- sum(parametersNumber[1,1:3]);
parametersNumber[2,4] <- sum(parametersNumber[2,1:3]);
}
#### Selection of the best model ####
else if(modelDo=="select"){
adamSelected <- selector(model, modelsPool, allowMultiplicative,
etsModel, Etype, Ttype, Stype, damped, lags,
lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted, ICFunction,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate);
icSelection <- adamSelected$icSelection;
# Take the parameters of the best model
list2env(adamSelected$results[[which.min(icSelection)[1]]], environment());
#### This part is needed in order for the filler to do its job later on
# Create the basic variables based on the estimated model
adamArchitect <- architector(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal,
xregNumber, obsInSample, initialType,
arimaModel, lagsModelARIMA, xregModel);
list2env(adamArchitect, environment());
# Create the matrices for the specific ETS model
adamCreated <- creator(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
lags, lagsModel, lagsModelARIMA, lagsModelAll, lagsModelMax,
obsStates, obsInSample, obsAll, componentsNumberETS, componentsNumberETSSeasonal,
componentsNamesETS, otLogical, yInSample,
persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate, persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi,
initialType, initialEstimate,
initialLevel, initialLevelEstimate, initialTrend, initialTrendEstimate,
initialSeasonal, initialSeasonalEstimate,
initialArima, initialArimaEstimate, initialArimaNumber,
initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
arOrders, iOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames);
list2env(adamCreated, environment());
parametersNumber[1,1] <- nParamEstimated;
if(xregModel){
parametersNumber[1,2] <- xregNumber*initialXregEstimate + xregNumber*persistenceXregEstimate;
}
parametersNumber[1,4] <- sum(parametersNumber[1,1:3]);
parametersNumber[2,4] <- sum(parametersNumber[2,1:3]);
}
#### Combination of models ####
else if(modelDo=="combine"){
modelOriginal <- model;
# If the pool is not provided, then create one
if(is.null(modelsPool)){
# Define the whole pool of errors
if(!allowMultiplicative){
poolErrors <- c("A");
poolTrends <- c("N","A","Ad");
poolSeasonals <- c("N","A");
}
else{
poolErrors <- c("A","M");
poolTrends <- c("N","A","Ad","M","Md");
poolSeasonals <- c("N","A","M");
}
# Some preparation variables
# If Etype is not Z, then check on additive errors
if(Etype!="Z"){
poolErrors <- switch(Etype,
"N"="N",
"A"=,
"X"="A",
"M"=,
"Y"="M");
}
# If Ttype is not Z, then create a pool with specified type
if(Ttype!="Z"){
poolTrends <- switch(Ttype,
"N"="N",
"A"=ifelse(damped,"Ad","A"),
"M"=ifelse(damped,"Md","M"),
"X"=c("N","A","Ad"),
"Y"=c("N","M","Md"));
}
# If Stype is not Z, then crete specific pools
if(Stype!="Z"){
poolSeasonals <- switch(Stype,
"N"="N",
"A"="A",
"X"=c("N","A"),
"M"="M",
"Y"=c("N","M"));
}
modelsPool <- paste0(rep(poolErrors,length(poolTrends)*length(poolSeasonals)),
rep(poolTrends,each=length(poolSeasonals)),
rep(poolSeasonals,length(poolTrends)));
}
adamSelected <- selector(model, modelsPool, allowMultiplicative,
etsModel, Etype, Ttype, Stype, damped, lags,
lagsModelSeasonal, lagsModelARIMA,
obsStates, obsInSample,
yInSample, persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate,
persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi, phiEstimate,
initialType, initialLevel, initialTrend, initialSeasonal,
initialArima, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames, xregDo,
ot, otLogical, occurrenceModel, pFitted, ICFunction,
bounds, loss, lossFunction, distribution,
horizon, multisteps, lambda, lambdaEstimate);
icSelection <- adamSelected$icSelection;
icBest <- min(icSelection);
adamSelected$icWeights <- (exp(-0.5*(icSelection-icBest)) /
sum(exp(-0.5*(icSelection-icBest))));
# This is a failsafe mechanism, just to make sure that the ridiculous models don't impact forecasts
adamSelected$icWeights[adamSelected$icWeights<1e-5] <- 0
adamSelected$icWeights <- adamSelected$icWeights/sum(adamSelected$icWeights);
# adamArchitect <- vector("list",10)
for(i in 1:length(adamSelected$results)){
# Take the parameters of the best model
list2env(adamSelected$results[[i]], environment());
#### This part is needed in order for the filler to do its job later on
# Create the basic variables based on the estimated model
adamArchitect <- architector(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal,
xregNumber, obsInSample, initialType,
arimaModel, lagsModelARIMA, xregModel);
list2env(adamArchitect, environment());
adamSelected$results[[i]]$modelIsTrendy <- adamArchitect$modelIsTrendy;
adamSelected$results[[i]]$modelIsSeasonal <- adamArchitect$modelIsSeasonal;
adamSelected$results[[i]]$lagsModel <- adamArchitect$lagsModel;
adamSelected$results[[i]]$lagsModelAll <- adamArchitect$lagsModelAll;
adamSelected$results[[i]]$lagsModelMax <- adamArchitect$lagsModelMax;
adamSelected$results[[i]]$componentsNumberETS <- adamArchitect$componentsNumberETS;
adamSelected$results[[i]]$componentsNumberETSSeasonal <- adamArchitect$componentsNumberETSSeasonal;
adamSelected$results[[i]]$componentsNumberETSNonSeasonal <- adamArchitect$componentsNumberETSNonSeasonal;
adamSelected$results[[i]]$componentsNamesETS <- adamArchitect$componentsNamesETS;
# Create the matrices for the specific ETS model
adamCreated <- creator(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
lags, lagsModel, lagsModelARIMA, lagsModelAll, lagsModelMax,
obsStates, obsInSample, obsAll, componentsNumberETS, componentsNumberETSSeasonal,
componentsNamesETS, otLogical, yInSample,
persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate, persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi,
initialType, initialEstimate,
initialLevel, initialLevelEstimate, initialTrend, initialTrendEstimate,
initialSeasonal, initialSeasonalEstimate,
initialArima, initialArimaEstimate, initialArimaNumber,
initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
arOrders, iOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames);
adamSelected$results[[i]]$matVt <- adamCreated$matVt;
adamSelected$results[[i]]$matWt <- adamCreated$matWt;
adamSelected$results[[i]]$matF <- adamCreated$matF;
adamSelected$results[[i]]$vecG <- adamCreated$vecG;
adamSelected$results[[i]]$arimaPolynomials <- adamCreated$arimaPolynomials;
parametersNumber[1,1] <- adamSelected$results[[i]]$nParamEstimated;
if(xregModel){
parametersNumber[1,2] <- xregNumber*initialXregEstimate + xregNumber*persistenceXregEstimate;
}
parametersNumber[1,4] <- sum(parametersNumber[1,1:3]);
parametersNumber[2,4] <- sum(parametersNumber[2,1:3]);
adamSelected$results[[i]]$parametersNumber <- parametersNumber;
}
}
#### Use the provided model ####
else if(modelDo=="use"){
# If the distribution is default, change it according to the error term
if(distribution=="default"){
distributionNew <- switch(Etype,
"A"=switch(loss,
"MAEh"=, "MACE"=, "MAE"="dlaplace",
"HAMh"=, "CHAM"=, "HAM"="ds",
"MSEh"=, "MSCE"=, "MSE"=, "GPL"=, "likelihood"=, "dnorm"),
"M"=switch(loss,
"MAEh"=, "MACE"=, "MAE"="dllaplace",
"HAMh"=, "CHAM"=, "HAM"="dls",
"MSEh"=, "MSCE"=, "MSE"=, "GPL"="dlnorm",
"likelihood"=, "dinvgauss"));
}
else{
distributionNew <- distribution;
}
# Create the basic variables
adamArchitect <- architector(etsModel, Etype, Ttype, Stype, lags, lagsModelSeasonal,
xregNumber, obsInSample, initialType,
arimaModel, lagsModelARIMA, xregModel);
list2env(adamArchitect, environment());
# Create the matrices for the specific ETS model
adamCreated <- creator(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
lags, lagsModel, lagsModelARIMA, lagsModelAll, lagsModelMax,
obsStates, obsInSample, obsAll, componentsNumberETS, componentsNumberETSSeasonal,
componentsNamesETS, otLogical, yInSample,
persistence, persistenceEstimate,
persistenceLevel, persistenceLevelEstimate, persistenceTrend, persistenceTrendEstimate,
persistenceSeasonal, persistenceSeasonalEstimate,
persistenceXreg, persistenceXregEstimate, persistenceXregProvided,
phi,
initialType, initialEstimate,
initialLevel, initialLevelEstimate, initialTrend, initialTrendEstimate,
initialSeasonal, initialSeasonalEstimate,
initialArima, initialArimaEstimate, initialArimaNumber,
initialXregEstimate, initialXregProvided,
arimaModel, arRequired, iRequired, maRequired, armaParameters,
arOrders, iOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA,
xregModel, xregModelInitials, xregData, xregNumber, xregNames);
list2env(adamCreated, environment());
CFValue <- CF(B=0, etsModel=etsModel, Etype=Etype, Ttype=Ttype, Stype=Stype, modelIsTrendy=modelIsTrendy,
modelIsSeasonal=modelIsSeasonal, yInSample=yInSample,
ot=ot, otLogical=otLogical, occurrenceModel=occurrenceModel, obsInSample=obsInSample,
componentsNumberETS=componentsNumberETS, componentsNumberETSSeasonal=componentsNumberETSSeasonal,
componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal,
componentsNumberARIMA=componentsNumberARIMA,
lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax,
matVt=matVt, matWt=matWt, matF=matF, vecG=vecG,
persistenceEstimate=persistenceEstimate,
persistenceLevelEstimate=persistenceLevelEstimate,
persistenceTrendEstimate=persistenceTrendEstimate,
persistenceSeasonalEstimate=persistenceSeasonalEstimate,
persistenceXregEstimate=persistenceXregEstimate,
phiEstimate=phiEstimate, initialType=initialType,
initialEstimate=initialEstimate, initialLevelEstimate=initialLevelEstimate,
initialTrendEstimate=initialTrendEstimate, initialSeasonalEstimate=initialSeasonalEstimate,
initialArimaEstimate=initialArimaEstimate, initialXregEstimate=initialXregEstimate,
arimaModel=arimaModel, nonZeroARI=nonZeroARI, nonZeroMA=nonZeroMA,
arimaPolynomials=arimaPolynomials,
arEstimate=arEstimate, maEstimate=maEstimate,
arOrders=arOrders, iOrders=iOrders, maOrders=maOrders,
arRequired=arRequired, maRequired=maRequired, armaParameters=armaParameters,
xregModel=xregModel, xregNumber=xregNumber,
bounds=bounds, loss=loss, lossFunction=lossFunction, distribution=distributionNew,
horizon=horizon, multisteps=multisteps, lambda=lambda, lambdaEstimate=lambdaEstimate,
arPolynomialMatrix=NULL, maPolynomialMatrix=NULL);
parametersNumber[1,1] <- parametersNumber[1,4] <- 1;
logLikADAMValue <- structure(logLikADAM(B=0,
etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal, yInSample,
ot, otLogical, occurrenceModel, pFitted, obsInSample,
componentsNumberETS, componentsNumberETSSeasonal, componentsNumberETSNonSeasonal,
componentsNumberARIMA,
lags, lagsModel, lagsModelAll, lagsModelMax,
matVt, matWt, matF, vecG,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
arimaModel, nonZeroARI, nonZeroMA, arEstimate, maEstimate, arimaPolynomials,
arOrders, iOrders, maOrders, arRequired, maRequired, armaParameters,
xregModel, xregNumber,
bounds, loss, lossFunction, distributionNew, horizon,
multisteps, lambda, lambdaEstimate, arPolynomialMatrix=NULL,
maPolynomialMatrix=NULL)
,nobs=obsInSample,df=parametersNumber[1,4],class="logLik")
icSelection <- ICFunction(logLikADAMValue);
# If Fisher Information is required, do that analytically
if(FI){
# If B is not provided, then use the standard thing
if(is.null(B)){
BValues <- initialiser(etsModel, Etype, Ttype, Stype, modelIsTrendy, modelIsSeasonal,
componentsNumberETSNonSeasonal, componentsNumberETSSeasonal, componentsNumberETS,
lags, lagsModel, lagsModelSeasonal, lagsModelARIMA, lagsModelMax,
matVt,
TRUE, TRUE, modelIsTrendy, rep(modelIsSeasonal,componentsNumberETSSeasonal), FALSE,
damped, "optimal", TRUE,
TRUE, TRUE, rep(modelIsSeasonal,componentsNumberETSSeasonal),
arimaModel, xregModel,
arimaModel, arRequired, maRequired, arRequired, maRequired, arOrders, maOrders,
componentsNumberARIMA, componentsNamesARIMA, initialArimaNumber,
xregModel, xregNumber, FALSE);
# Create the vector of initials for the optimisation
B <- BValues$B;
}
# Define parameters just for FI calculation
if(initialType=="provided"){
initialLevelEstimateFI <- any(names(B)=="level");
initialTrendEstimateFI <- any(names(B)=="trend");
if(any(substr(names(B),1,8)=="seasonal")){
initialSeasonalEstimateFI <- vector("logical", componentsNumberETSSeasonal);
seasonalNames <- names(B)[substr(names(B),1,8)=="seasonal"];
# If there is only one seasonality
if(substr(seasonalNames,1,9)=="seasonal_"){
initialSeasonalEstimateFI[] <- TRUE;
}
# If there is several
else{
initialSeasonalEstimateFI[unique(as.numeric(substr(seasonalNames,9,9)))] <- TRUE;
}
}
else{
initialSeasonalEstimateFI <- FALSE;
}
if(arimaModel){
initialArimaEstimateFI <- any(substr(names(B),1,10)=="ARIMAState");
}
else{
initialArimaEstimateFI <- FALSE;
}
if(xregModel){
initialXregEstimateFI <- any(colnames(xregData) %in% names(B));
}
else{
initialXregEstimateFI <- FALSE;
}
initialTypeFI <- "optimal";
initialEstimateFI <- any(c(initialLevelEstimateFI,initialTrendEstimateFI,initialSeasonalEstimateFI,
initialArimaEstimateFI, initialXregEstimateFI));
}
else{
initialTypeFI <- initialType;
initialEstimateFI <- FALSE;
}
# If smoothing parmaeters were estimated, then alpha should be in the list
persistenceLevelEstimateFI <- any(names(B)=="alpha");
persistenceTrendEstimateFI <- any(names(B)=="beta");
if(any(substr(names(B),1,5)=="gamma")){
gammas <- (substr(names(B),1,5)=="gamma");
if(sum(gammas)==1){
persistenceSeasonalEstimateFI <- TRUE;
}
else{
persistenceSeasonalEstimateFI <- vector("logical",componentsNumberETSSeasonal);
persistenceSeasonalEstimateFI[as.numeric(substr(names(B),6,6)[gammas])] <- TRUE;
}
}
else{
persistenceSeasonalEstimateFI <- FALSE;
}
persistenceXregEstimateFI <- any(substr(names(B),1,5)=="delta");
persistenceEstimateFI <- any(c(persistenceLevelEstimateFI,persistenceTrendEstimateFI,
persistenceSeasonalEstimateFI,persistenceXregEstimateFI));
phiEstimateFI <- any(names(B)=="phi");
lambdaEstimateFI <- any(names(B)=="lambda");
if(arimaModel){
maEstimateFI <- maRequired;
arEstimateFI <- arRequired;
maPolynomialMatrix <- arPolynomialMatrix <- NULL;
}
FI <- -hessian(logLikADAM, B, etsModel=etsModel, Etype=Etype, Ttype=Ttype, Stype=Stype, modelIsTrendy=modelIsTrendy,
modelIsSeasonal=modelIsSeasonal, yInSample=yInSample,
ot=ot, otLogical=otLogical, occurrenceModel=occurrenceModel, pFitted=pFitted, obsInSample=obsInSample,
componentsNumberETS=componentsNumberETS, componentsNumberETSSeasonal=componentsNumberETSSeasonal,
componentsNumberETSNonSeasonal=componentsNumberETSNonSeasonal,
componentsNumberARIMA=componentsNumberARIMA,
lags=lags, lagsModel=lagsModel, lagsModelAll=lagsModelAll, lagsModelMax=lagsModelMax,
matVt=matVt, matWt=matWt, matF=matF, vecG=vecG,
persistenceEstimate=persistenceEstimateFI, persistenceLevelEstimate=persistenceLevelEstimateFI,
persistenceTrendEstimate=persistenceTrendEstimateFI,
persistenceSeasonalEstimate=persistenceSeasonalEstimateFI,
persistenceXregEstimate=persistenceXregEstimateFI,
phiEstimate=phiEstimateFI, initialType=initialTypeFI,
initialEstimate=initialEstimateFI, initialLevelEstimate=initialLevelEstimateFI,
initialTrendEstimate=initialTrendEstimateFI, initialSeasonalEstimate=initialSeasonalEstimateFI,
initialArimaEstimate=initialArimaEstimateFI, initialXregEstimate=initialXregEstimateFI,
arimaModel=arimaModel, nonZeroARI=nonZeroARI, nonZeroMA=nonZeroMA,
arEstimate=arEstimateFI, maEstimate=maEstimateFI, arimaPolynomials=arimaPolynomials,
arOrders=arOrders, iOrders=iOrders, maOrders=maOrders,
arRequired=arRequired, maRequired=maRequired, armaParameters=armaParameters,
xregModel=xregModel, xregNumber=xregNumber,
bounds=bounds, loss=loss, lossFunction=lossFunction, distribution=distribution,
horizon=horizon, multisteps=multisteps,
lambda=lambda, lambdaEstimate=lambdaEstimateFI,
arPolynomialMatrix=arPolynomialMatrix, maPolynomialMatrix=maPolynomialMatrix);
colnames(FI) <- names(B);
rownames(FI) <- names(B);
}
else{
FI <- NULL;
}
}
# Transform everything into appropriate classes
if(any(yClasses=="ts")){
yInSample <- ts(yInSample,start=yStart, frequency=yFrequency);
if(holdout){
yHoldout <- ts(yHoldout, start=yForecastStart, frequency=yFrequency);
}
}
else{
yInSample <- zoo(yInSample, order.by=yInSampleIndex);
if(holdout){
yHoldout <- zoo(yHoldout, order.by=yForecastIndex);
}
}
#### Prepare the return if we didn't combine anything ####
if(modelDo!="combine"){
modelReturned <- preparator(B, etsModel, Etype, Ttype, Stype,
lagsModel, lagsModelMax, lagsModelAll,
componentsNumberETS, componentsNumberETSSeasonal,
xregNumber, distribution, loss,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, lambdaEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
matVt, matWt, matF, vecG,
occurrenceModel, ot, oesModel,
parametersNumber, CFValue,
arimaModel, arRequired, maRequired,
arEstimate, maEstimate, arOrders, iOrders, maOrders,
nonZeroARI, nonZeroMA,
arimaPolynomials, armaParameters);
# Prepare the name of the model
modelName <- "";
if(etsModel){
if(model!="NNN"){
modelName[] <- "ETS";
if(xregModel){
modelName[] <- paste0(modelName,"X");
}
modelName[] <- paste0(modelName,"(",model,")");
if(componentsNumberETSSeasonal>1){
modelName[] <- paste0(modelName,"[",paste0(lags[lags!=1], collapse=", "),"]");
}
}
else{
modelName[] <- "Constant level";
}
}
if(arimaModel){
if(etsModel){
modelName[] <- paste0(modelName,"+");
}
# Either the lags are non-seasonal, or there are no orders for seasonal lags
if(all(lags==1) || (all(arOrders[lags>1]==0) && all(iOrders[lags>1]==0) && all(maOrders[lags>1]==0))){
modelName[] <- paste0(modelName,"ARIMA");
if(!etsModel && xregModel){
modelName[] <- paste0(modelName,"X");
}
modelName[] <- paste0(modelName,"(",arOrders[1],",",iOrders[1],",",maOrders[1],")");
}
else{
modelName[] <- paste0(modelName,"SARIMA");
if(!etsModel && xregModel){
modelName[] <- paste0(modelName,"X");
}
for(i in 1:length(arOrders)){
if(all(arOrders[i]==0) && all(iOrders[i]==0) && all(maOrders[i]==0)){
next;
}
modelName[] <- paste0(modelName,"(",arOrders[i],",");
modelName[] <- paste0(modelName,iOrders[i],",");
modelName[] <- paste0(modelName,maOrders[i],")[",lags[i],"]");
}
}
}
if(!etsModel && !arimaModel){
if(xregDo=="adapt"){
modelName[] <- paste0("Dynamic regression");
}
else{
modelName[] <- paste0("Regression");
}
}
if(all(occurrence!=c("n","none"))){
modelName[] <- paste0("i",modelName);
}
modelReturned$model <- modelName;
modelReturned$timeElapsed <- Sys.time()-startTime;
modelReturned$y <- yInSample;
modelReturned$holdout <- yHoldout;
if(any(yNAValues)){
modelReturned$y[yNAValues[1:obsInSample]] <- NA;
if(length(yNAValues)==obsAll){
modelReturned$holdout[yNAValues[-c(1:obsInSample)]] <- NA;
}
modelReturned$residuals[yNAValues[1:obsInSample]] <- NA;
}
class(modelReturned) <- c("adam","smooth");
}
#### Return the combined model ####
else{
modelReturned <- list(models=vector("list",length(adamSelected$results)));
yFittedCombined <- rep(0,obsInSample);
if(h>0){
yForecastCombined <- rep(0,h);
}
else{
yForecastCombined <- NA;
}
parametersNumberOverall <- parametersNumber;
for(i in 1:length(adamSelected$results)){
list2env(adamSelected$results[[i]], environment());
modelReturned$models[[i]] <- preparator(B, etsModel, Etype, Ttype, Stype,
lagsModel, lagsModelMax, lagsModelAll,
componentsNumberETS, componentsNumberETSSeasonal,
xregNumber, distribution, loss,
persistenceEstimate, persistenceLevelEstimate, persistenceTrendEstimate,
persistenceSeasonalEstimate, persistenceXregEstimate,
phiEstimate, lambdaEstimate,
initialType, initialEstimate,
initialLevelEstimate, initialTrendEstimate, initialSeasonalEstimate,
initialArimaEstimate, initialXregEstimate,
matVt, matWt, matF, vecG,
occurrenceModel, ot, oesModel,
parametersNumber, CFValue,
arimaModel, arRequired, maRequired,
arEstimate, maEstimate, arOrders, iOrders, maOrders,
nonZeroARI, nonZeroMA,
arimaPolynomials, armaParameters);
modelReturned$models[[i]]$fitted[is.na(modelReturned$models[[i]]$fitted)] <- 0;
yFittedCombined[] <- yFittedCombined + modelReturned$models[[i]]$fitted * adamSelected$icWeights[i];
if(h>0){
modelReturned$models[[i]]$forecast[is.na(modelReturned$models[[i]]$forecast)] <- 0;
yForecastCombined[] <- yForecastCombined + modelReturned$models[[i]]$forecast * adamSelected$icWeights[i];
}
# Prepare the name of the model
modelName <- "";
if(xregModel){
modelName[] <- "ETSX";
}
else{
modelName[] <- "ETS";
}
modelName[] <- paste0(modelName,"(",model,")");
if(all(occurrence!=c("n","none"))){
modelName[] <- paste0("i",modelName);
}
if(componentsNumberETSSeasonal>1){
modelName[] <- paste0(modelName,"[",paste0(lags[lags!=1], collapse=", "),"]");
}
if(arimaModel){
# Either the lags are non-seasonal, or there are no orders for seasonal lags
if(all(lags==1) || (all(arOrders[lags>1]==0) && all(iOrders[lags>1]==0) && all(maOrders[lags>1]==0))){
modelName[] <- paste0(modelName,"+ARIMA(",arOrders[1],",",iOrders[1],",",maOrders[1],")");
}
else{
modelName[] <- paste0(modelName,"+SARIMA");
for(i in 1:length(arOrders)){
if(all(arOrders[i]==0) && all(iOrders[i]==0) && all(maOrders[i]==0)){
next;
}
modelName[] <- paste0(modelName,"(",arOrders[i],",");
modelName[] <- paste0(modelName,iOrders[i],",");
modelName[] <- paste0(modelName,maOrders[i],")[",lags[i],"]");
}
}
}
if(all(occurrence!=c("n","none"))){
modelName[] <- paste0("i",modelName);
}
modelReturned$models[[i]]$model <- modelName;
modelReturned$models[[i]]$timeElapsed <- Sys.time()-startTime;
parametersNumberOverall[1,1] <- parametersNumber[1,1] + parametersNumber[1,1] * adamSelected$icWeights[i];
modelReturned$models[[i]]$y <- yInSample;
if(any(yNAValues)){
modelReturned$models[[i]]$y[yNAValues[1:obsInSample]] <- NA;
if(length(yNAValues)==obsAll){
modelReturned$models[[i]]$holdout[yNAValues[-c(1:obsInSample)]] <- NA;
}
modelReturned$models[[i]]$residuals[yNAValues[1:obsInSample]] <- NA;
}
class(modelReturned$models[[i]]) <- c("adam","smooth");
}
# Record the original name of the model.
model[] <- modelOriginal;
# Prepare the name of the model
modelName <- "";
if(xregModel){
modelName[] <- "ETSX";
}
else{
modelName[] <- "ETS";
}
modelName[] <- paste0(modelName,"(",model,")");
if(all(occurrence!=c("n","none"))){
modelName[] <- paste0("i",modelName);
}
if(componentsNumberETSSeasonal>1){
modelName[] <- paste0(modelName,"[",paste0(lags[lags!=1], collapse=", "),"]");
}
if(arimaModel){
# Either the lags are non-seasonal, or there are no orders for seasonal lags
if(all(lags==1) || (all(arOrders[lags>1]==0) && all(iOrders[lags>1]==0) && all(maOrders[lags>1]==0))){
modelName[] <- paste0(modelName,"+ARIMA(",arOrders[1],",",iOrders[1],",",maOrders[1],")");
}
else{
modelName[] <- paste0(modelName,"+SARIMA");
for(i in 1:length(arOrders)){
if(all(arOrders[i]==0) && all(iOrders[i]==0) && all(maOrders[i]==0)){
next;
}
modelName[] <- paste0(modelName,"(",arOrders[i],",");
modelName[] <- paste0(modelName,iOrders[i],",");
modelName[] <- paste0(modelName,maOrders[i],")[",lags[i],"]");
}
}
}
if(all(occurrence!=c("n","none"))){
modelName[] <- paste0("i",modelName);
}
modelReturned$model <- modelName;
modelReturned$timeElapsed <- Sys.time()-startTime;
modelReturned$holdout <- yHoldout;
modelReturned$y <- yInSample;
modelReturned$fitted <- ts(yFittedCombined,start=yStart, frequency=yFrequency);
modelReturned$residuals <- yInSample - yFittedCombined;
if(any(yNAValues)){
modelReturned$y[yNAValues[1:obsInSample]] <- NA;
if(length(yNAValues)==obsAll){
modelReturned$holdout[yNAValues[-c(1:obsInSample)]] <- NA;
}
modelReturned$residuals[yNAValues[1:obsInSample]] <- NA;
}
modelReturned$forecast <- ts(yForecastCombined,start=yForecastStart, frequency=yFrequency);
parametersNumberOverall[1,4] <- sum(parametersNumberOverall[1,1:3]);
parametersNumberOverall[2,4] <- sum(parametersNumberOverall[2,1:3]);
modelReturned$nParam <- parametersNumberOverall;
modelReturned$ICw <- adamSelected$icWeights;
# These two are needed just to make basic methods work
modelReturned$distribution <- distribution;
modelReturned$scale <- sqrt(mean(modelReturned$residuals^2,na.rm=TRUE));
class(modelReturned) <- c("adamCombined","adam","smooth");
}
modelReturned$ICs <- icSelection;
modelReturned$lossFunction <- lossFunction;
# Error measures if there is a holdout
if(holdout){
modelReturned$accuracy <- measures(yHoldout,modelReturned$forecast,yInSample);
}
if(!silent){
plot(modelReturned, 7);
}
return(modelReturned);
}
#### Technical methods ####
#' @export
lags.adam <- function(object, ...){
return(object$lags);
}
#' @rdname plot.smooth
#' @export
plot.adam <- function(x, which=c(1,2,4,6), level=0.95, legend=FALSE,
ask=prod(par("mfcol")) < length(which) && dev.interactive(),
lowess=TRUE, ...){
ellipsis <- list(...);
# Define, whether to wait for the hit of "Enter"
if(ask){
oask <- devAskNewPage(TRUE);
on.exit(devAskNewPage(oask));
}
# 1. Fitted vs Actuals values
plot1 <- function(x, ...){
ellipsis <- list(...);
# Get the actuals and the fitted values
ellipsis$y <- as.vector(actuals(x));
if(is.occurrence(x)){
if(any(x$distribution==c("plogis","pnorm")) || x$occurrence!="none"){
ellipsis$y <- (ellipsis$y!=0)*1;
}
}
ellipsis$x <- as.vector(fitted(x));
# If this is a mixture model, remove zeroes
if(is.occurrence(x$occurrence)){
ellipsis$x <- ellipsis$x[ellipsis$y!=0];
ellipsis$y <- ellipsis$y[ellipsis$y!=0];
}
# Remove NAs
if(any(is.na(ellipsis$x))){
ellipsis$y <- ellipsis$y[!is.na(ellipsis$x)];
ellipsis$x <- ellipsis$x[!is.na(ellipsis$x)];
}
if(any(is.na(ellipsis$y))){
ellipsis$x <- ellipsis$x[!is.na(ellipsis$y)];
ellipsis$y <- ellipsis$y[!is.na(ellipsis$y)];
}
# Title
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "Actuals vs Fitted";
}
# If type and ylab are not provided, set them...
if(!any(names(ellipsis)=="type")){
ellipsis$type <- "p";
}
if(!any(names(ellipsis)=="ylab")){
ellipsis$ylab <- "Actuals";
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Fitted";
}
# xlim and ylim
if(!any(names(ellipsis)=="xlim")){
ellipsis$xlim <- range(c(ellipsis$x,ellipsis$y));
}
if(!any(names(ellipsis)=="ylim")){
ellipsis$ylim <- range(c(ellipsis$x,ellipsis$y));
}
# Start plotting
do.call(plot,ellipsis);
abline(a=0,b=1,col="grey",lwd=2,lty=2)
if(lowess){
lines(lowess(ellipsis$x, ellipsis$y), col="red");
}
}
# 2 and 3: Standardised / studentised residuals vs Fitted
plot2 <- function(x, type="rstandard", ...){
ellipsis <- list(...);
ellipsis$x <- as.vector(fitted(x));
if(type=="rstandard"){
ellipsis$y <- as.vector(rstandard(x));
yName <- "Standardised";
}
else{
ellipsis$y <- as.vector(rstudent(x));
yName <- "Studentised";
}
if(is.occurrence(x$occurrence)){
ellipsis$x <- ellipsis$x[actuals(x$occurrence)!=0];
ellipsis$y <- ellipsis$y[actuals(x$occurrence)!=0];
}
# Remove NAs
if(any(is.na(ellipsis$x))){
ellipsis$x <- ellipsis$x[!is.na(ellipsis$x)];
ellipsis$y <- ellipsis$y[!is.na(ellipsis$y)];
}
# Main, labs etc
if(!any(names(ellipsis)=="main")){
if(any(x$distribution==c("dinvgauss","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$main <- paste0("log(",yName," Residuals) vs Fitted");
}
else{
ellipsis$main <- paste0(yName," Residuals vs Fitted");
}
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Fitted";
}
if(!any(names(ellipsis)=="ylab")){
ellipsis$ylab <- paste0(yName," Residuals");
}
if(legend){
if(ellipsis$x[length(ellipsis$x)]>mean(ellipsis$x)){
legendPosition <- "bottomright";
}
else{
legendPosition <- "topright";
}
}
zValues <- switch(x$distribution,
"dlaplace"=,
"dllaplace"=qlaplace(c((1-level)/2, (1+level)/2), 0, 1),
"dalaplace"=qalaplace(c((1-level)/2, (1+level)/2), 0, 1, x$lambda),
"dlogis"=qlogis(c((1-level)/2, (1+level)/2), 0, 1),
"dt"=qt(c((1-level)/2, (1+level)/2), nobs(x)-nparam(x)),
"ds"=,
"dls"=qs(c((1-level)/2, (1+level)/2), 0, 1),
"dgnorm"=,
"dlgnorm"=qgnorm(c((1-level)/2, (1+level)/2), 0, 1, x$lambda),
# In the next one, the scale is debiased, taking n-k into account
"dinvgauss"=qinvgauss(c((1-level)/2, (1+level)/2), mean=1,
dispersion=x$scale * nobs(x) / (nobs(x)-nparam(x))),
"dlnorm"=qlnorm(c((1-level)/2, (1+level)/2), 0, 1),
qnorm(c((1-level)/2, (1+level)/2), 0, 1));
# Analyse stuff in logarithms if the error is multiplicative
if(any(x$distribution==c("dinvgauss","dlnorm"))){
ellipsis$y[] <- log(ellipsis$y);
zValues <- log(zValues);
}
else if(any(x$distribution==c("dllaplace","dls","dlgnorm"))){
ellipsis$y[] <- log(ellipsis$y);
}
outliers <- which(ellipsis$y >zValues[2] | ellipsis$y <zValues[1]);
# cat(paste0(round(length(outliers)/length(ellipsis$y),3)*100,"% of values are outside the bounds\n"));
if(!any(names(ellipsis)=="ylim")){
ellipsis$ylim <- range(c(ellipsis$y,zValues), na.rm=TRUE)*1.2;
if(legend){
if(legendPosition=="bottomright"){
ellipsis$ylim[1] <- ellipsis$ylim[1] - 0.2*diff(ellipsis$ylim);
}
else{
ellipsis$ylim[2] <- ellipsis$ylim[2] + 0.2*diff(ellipsis$ylim);
}
}
}
xRange <- range(ellipsis$x, na.rm=TRUE);
xRange[1] <- xRange[1] - sd(ellipsis$x, na.rm=TRUE);
xRange[2] <- xRange[2] + sd(ellipsis$x, na.rm=TRUE);
do.call(plot,ellipsis);
abline(h=0, col="grey", lty=2);
polygon(c(xRange,rev(xRange)),c(zValues[1],zValues[1],zValues[2],zValues[2]),
col="lightgrey", border=NA, density=10);
abline(h=zValues, col="red", lty=2);
if(length(outliers)>0){
points(ellipsis$x[outliers], ellipsis$y[outliers], pch=16);
text(ellipsis$x[outliers], ellipsis$y[outliers], labels=outliers, pos=(ellipsis$y[outliers]>0)*2+1);
}
if(lowess){
lines(lowess(ellipsis$x[!is.na(ellipsis$y)], ellipsis$y[!is.na(ellipsis$y)]), col="red");
}
if(legend){
if(lowess){
legend(legendPosition,
legend=c(paste0(round(level,3)*100,"% bounds"),"outside the bounds","LOWESS line"),
col=c("red", "black","red"), lwd=c(1,NA,1), lty=c(2,1,1), pch=c(NA,16,NA));
}
else{
legend(legendPosition,
legend=c(paste0(round(level,3)*100,"% bounds"),"outside the bounds"),
col=c("red", "black"), lwd=c(1,NA), lty=c(2,1), pch=c(NA,16));
}
}
}
# 4 and 5. Fitted vs |Residuals| or Fitted vs Residuals^2
plot3 <- function(x, type="abs", ...){
ellipsis <- list(...);
ellipsis$x <- as.vector(fitted(x));
ellipsis$y <- as.vector(residuals(x));
if(any(x$distribution==c("dinvgauss","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$y[] <- log(ellipsis$y);
}
if(type=="abs"){
ellipsis$y[] <- abs(ellipsis$y);
}
else{
ellipsis$y[] <- as.vector(ellipsis$y)^2;
}
if(is.occurrence(x$occurrence)){
ellipsis$x <- ellipsis$x[ellipsis$y!=0];
ellipsis$y <- ellipsis$y[ellipsis$y!=0];
}
# Remove NAs
if(any(is.na(ellipsis$x))){
ellipsis$x <- ellipsis$x[!is.na(ellipsis$x)];
ellipsis$y <- ellipsis$y[!is.na(ellipsis$y)];
}
if(!any(names(ellipsis)=="main")){
if(type=="abs"){
if(any(x$distribution==c("dinvgauss","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$main <- "|log(Residuals)| vs Fitted";
}
else{
ellipsis$main <- "|Residuals| vs Fitted";
}
}
else{
if(any(x$distribution==c("dinvgauss","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$main <- "log(Residuals)^2 vs Fitted";
}
else{
ellipsis$main <- "Residuals^2 vs Fitted";
}
}
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Fitted";
}
if(!any(names(ellipsis)=="ylab")){
if(type=="abs"){
ellipsis$ylab <- "|Residuals|";
}
else{
ellipsis$ylab <- "Residuals^2";
}
}
do.call(plot,ellipsis);
abline(h=0, col="grey", lty=2);
if(lowess){
lines(lowess(ellipsis$x[!is.na(ellipsis$y)], ellipsis$y[!is.na(ellipsis$y)]), col="red");
}
}
# 6. Q-Q with the specified distribution
plot4 <- function(x, ...){
ellipsis <- list(...);
ellipsis$y <- as.vector(residuals(x));
if(is.occurrence(x$occurrence)){
ellipsis$y <- ellipsis$y[actuals(x$occurrence)!=0];
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Theoretical Quantile";
}
if(!any(names(ellipsis)=="ylab")){
ellipsis$ylab <- "Actual Quantile";
}
if(any(x$distribution=="dnorm")){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ plot of Normal distribution";
}
do.call(qqnorm, ellipsis);
qqline(ellipsis$y);
}
else if(any(x$distribution=="dlnorm")){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ plot of Log Normal distribution";
}
ellipsis$x <- qlnorm(ppoints(500), meanlog=0, sdlog=x$scale);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qlnorm(p, meanlog=0, sdlog=x$scale));
}
else if(x$distribution=="dlaplace"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Laplace distribution";
}
ellipsis$x <- qlaplace(ppoints(500), mu=0, scale=x$scale);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qlaplace(p, mu=0, scale=x$scale));
}
else if(x$distribution=="dllaplace"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Log Laplace distribution";
}
ellipsis$x <- exp(qlaplace(ppoints(500), mu=0, scale=x$scale));
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) exp(qlaplace(p, mu=0, scale=x$scale)));
}
else if(x$distribution=="ds"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of S distribution";
}
ellipsis$x <- qs(ppoints(500), mu=0, scale=x$scale);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qs(p, mu=0, scale=x$scale));
}
else if(x$distribution=="dls"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Log S distribution";
}
ellipsis$x <- exp(qs(ppoints(500), mu=0, scale=x$scale));
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) exp(qs(p, mu=0, scale=x$scale)));
}
else if(x$distribution=="dgnorm"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Generalised Normal distribution";
}
ellipsis$x <- qgnorm(ppoints(500), mu=0, alpha=x$scale, beta=x$lambda);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qgnorm(p, mu=0, alpha=x$scale, beta=x$lambda));
}
else if(x$distribution=="dlgnorm"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Log Generalised Normal distribution";
}
ellipsis$x <- exp(qgnorm(ppoints(500), mu=0, lpha=x$scale, beta=x$lambda));
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) exp(qgnorm(p, mu=0, lpha=x$scale, beta=x$lambda)));
}
else if(x$distribution=="dlogis"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Logistic distribution";
}
ellipsis$x <- qlogis(ppoints(500), location=0, scale=x$scale);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qlogis(p, location=0, scale=x$scale));
}
else if(x$distribution=="dt"){
# Standardise residuals
ellipsis$y[] <- ellipsis$y / sd(ellipsis$y);
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Student's distribution";
}
ellipsis$x <- qt(ppoints(500), df=x$scale);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qt(p, df=x$scale));
}
else if(x$distribution=="dalaplace"){
if(!any(names(ellipsis)=="main")){
ellipsis$main <- paste0("QQ-plot of Asymmetric Laplace with alpha=",round(x$lambda,3));
}
ellipsis$x <- qalaplace(ppoints(500), mu=0, scale=x$scale, alpha=x$lambda);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qalaplace(p, mu=0, scale=x$scale, alpha=x$lambda));
}
else if(x$distribution=="dinvgauss"){
# Transform residuals for something meaningful
# This is not 100% accurate, because the dispersion should change as well as mean...
if(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ-plot of Inverse Gaussian distribution";
}
ellipsis$x <- qinvgauss(ppoints(500), mean=1, dispersion=x$scale);
do.call(qqplot, ellipsis);
qqline(ellipsis$y, distribution=function(p) qinvgauss(p, mean=1, dispersion=x$scale));
}
}
# 7. Basic plot over time
plot5 <- function(x, ...){
ellipsis <- list(...);
ellipsis$actuals <- actuals(x);
if(!is.null(x$holdout)){
if(is.zoo(ellipsis$actuals)){
ellipsis$actuals <- zoo(c(as.vector(ellipsis$actuals),as.vector(x$holdout)),
order.by=c(time(ellipsis$actuals),time(x$holdout)));
}
else{
ellipsis$actuals <- ts(c(ellipsis$actuals,x$holdout),
start=start(ellipsis$actuals),
frequency=frequency(ellipsis$actuals));
}
}
if(is.null(ellipsis$main)){
ellipsis$main <- x$model;
}
ellipsis$forecast <- x$forecast;
ellipsis$fitted <- fitted(x);
ellipsis$legend <- FALSE;
ellipsis$parReset <- FALSE;
do.call(graphmaker, ellipsis);
}
# 8 and 9. Standardised / Studentised residuals vs time
plot6 <- function(x, type="rstandard", ...){
ellipsis <- list(...);
if(type=="rstandard"){
ellipsis$x <- rstandard(x);
yName <- "Standardised";
}
else{
ellipsis$x <- rstudent(x);
yName <- "Studentised";
}
if(is.occurrence(x$occurrence)){
ellipsis$x <- ellipsis$x[actuals(x$occurrence)!=0];
}
if(!any(names(ellipsis)=="main")){
ellipsis$main <- paste0(yName," Residuals vs Time");
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Time";
}
if(!any(names(ellipsis)=="ylab")){
ellipsis$ylab <- paste0(yName," Residuals");
}
# If type and ylab are not provided, set them...
if(!any(names(ellipsis)=="type")){
ellipsis$type <- "l";
}
zValues <- switch(x$distribution,
"dlnorm"=qlnorm(c((1-level)/2, (1+level)/2), 0, 1),
"dlaplace"=,
"dllaplace"=qlaplace(c((1-level)/2, (1+level)/2), 0, 1),
"ds"=,
"dls"=qs(c((1-level)/2, (1+level)/2), 0, 1),
"dgnorm"=,
"dlgnorm"=qgnorm(c((1-level)/2, (1+level)/2), 0, 1, x$lambda),
"dlogis"=qlogis(c((1-level)/2, (1+level)/2), 0, 1),
"dt"=qt(c((1-level)/2, (1+level)/2), nobs(x)-nparam(x)),
"dalaplace"=qalaplace(c((1-level)/2, (1+level)/2), 0, 1, x$lambda),
# In the next one, the scale is debiased, taking n-k into account
"dinvgauss"=qinvgauss(c((1-level)/2, (1+level)/2), mean=1,
dispersion=x$scale * nobs(x) / (nobs(x)-nparam(x))),
qnorm(c((1-level)/2, (1+level)/2), 0, 1));
# Analyse stuff in logarithms if the error is multiplicative
if(any(x$distribution==c("dinvgauss","dlnorm"))){
ellipsis$x[] <- log(ellipsis$x);
zValues <- log(zValues);
}
else if(any(x$distribution==c("dllaplace","dls","dlgnorm"))){
ellipsis$x[] <- log(ellipsis$x);
}
outliers <- which(ellipsis$x >zValues[2] | ellipsis$x <zValues[1]);
if(!any(names(ellipsis)=="ylim")){
ellipsis$ylim <- c(-max(abs(ellipsis$x)),max(abs(ellipsis$x)))*1.2;
}
if(legend){
legendPosition <- "topright";
ellipsis$ylim[2] <- ellipsis$ylim[2] + 0.2*diff(ellipsis$ylim);
ellipsis$ylim[1] <- ellipsis$ylim[1] - 0.2*diff(ellipsis$ylim);
}
# Start plotting
do.call(plot,ellipsis);
if(length(outliers)>0){
points(time(ellipsis$x)[outliers], ellipsis$x[outliers], pch=16);
text(time(ellipsis$x)[outliers], ellipsis$x[outliers], labels=outliers, pos=(ellipsis$x[outliers]>0)*2+1);
}
if(lowess){
lines(lowess(c(1:length(ellipsis$x)),ellipsis$x), col="red");
}
abline(h=0, col="grey", lty=2);
abline(h=zValues[1], col="red", lty=2);
abline(h=zValues[2], col="red", lty=2);
polygon(c(1:nobs(x), c(nobs(x):1)),
c(rep(zValues[1],nobs(x)), rep(zValues[2],nobs(x))),
col="lightgrey", border=NA, density=10);
if(legend){
legend(legendPosition,legend=c("Residuals",paste0(level*100,"% prediction interval")),
col=c("black","red"), lwd=rep(1,3), lty=c(1,1,2));
}
}
# 10 and 11. ACF and PACF
plot7 <- function(x, type="acf", ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="main")){
if(type=="acf"){
ellipsis$main <- "Autocorrelation Function of Residuals";
}
else{
ellipsis$main <- "Partial Autocorrelation Function of Residuals";
}
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Lags";
}
if(!any(names(ellipsis)=="ylab")){
if(type=="acf"){
ellipsis$ylab <- "ACF";
}
else{
ellipsis$ylab <- "PACF";
}
}
if(!any(names(ellipsis)=="ylim")){
ellipsis$ylim <- c(-1,1);
}
if(type=="acf"){
theValues <- acf(as.vector(residuals(x)), plot=FALSE, na.action=na.pass);
}
else{
theValues <- pacf(as.vector(residuals(x)), plot=FALSE, na.action=na.pass);
}
ellipsis$x <- theValues$acf[-1];
zValues <- qnorm(c((1-level)/2, (1+level)/2),0,sqrt(1/nobs(x)));
ellipsis$type <- "h"
do.call(plot,ellipsis);
abline(h=0, col="black", lty=1);
abline(h=zValues, col="red", lty=2);
if(any(ellipsis$x>zValues[2] | ellipsis$x<zValues[1])){
outliers <- which(ellipsis$x >zValues[2] | ellipsis$x <zValues[1]);
points(outliers, ellipsis$x[outliers], pch=16);
text(outliers, ellipsis$x[outliers], labels=outliers, pos=(ellipsis$x[outliers]>0)*2+1);
}
}
# 12. Plot of states
plot8 <- function(x, ...){
parDefault <- par(no.readonly = TRUE);
if(any(unlist(gregexpr("C",x$model))==-1)){
statesNames <- c("actuals",colnames(x$states),"residuals");
x$states <- cbind(actuals(x),x$states,residuals(x));
colnames(x$states) <- statesNames;
if(ncol(x$states)>10){
message("Too many states. Plotting them one by one on several graphs.");
if(is.null(ellipsis$main)){
ellipsisMain <- NULL;
}
else{
ellipsisMain <- ellipsis$main;
}
nPlots <- ceiling(ncol(x$states)/10);
for(i in 1:nPlots){
if(is.null(ellipsisMain)){
ellipsis$main <- paste0("States of ",x$model,", part ",i);
}
ellipsis$x <- x$states[,(1+(i-1)*10):min(i*10,ncol(x$states)),drop=FALSE];
do.call(plot, ellipsis);
}
}
else{
if(ncol(x$states)<=5){
ellipsis$nc <- 1;
}
if(is.null(ellipsis$main)){
ellipsis$main <- paste0("States of ",x$model);
}
ellipsis$x <- x$states;
do.call(plot, ellipsis);
}
}
else{
# If we did combinations, we cannot return anything
message("Combination of models was done. Sorry, but there is nothing to plot.");
}
par(parDefault);
}
# Do plots
if(any(which==1)){
plot1(x, ...);
}
if(any(which==2)){
plot2(x, ...);
}
if(any(which==3)){
plot2(x, "rstudent", ...);
}
if(any(which==4)){
plot3(x, ...);
}
if(any(which==5)){
plot3(x, type="squared", ...);
}
if(any(which==6)){
plot4(x, ...);
}
if(any(which==7)){
plot5(x, ...);
}
if(any(which==8)){
plot6(x, ...);
}
if(any(which==9)){
plot6(x, "rstudent", ...);
}
if(any(which==10)){
plot7(x, type="acf", ...);
}
if(any(which==11)){
plot7(x, type="pacf", ...);
}
if(any(which==12)){
plot8(x, ...);
}
}
#' @export
print.adam <- function(x, digits=4, ...){
etsModel <- any(unlist(gregexpr("ETS",x$model))!=-1);
arimaModel <- any(unlist(gregexpr("ARIMA",x$model))!=-1);
cat(paste0("Time elapsed: ",round(as.numeric(x$timeElapsed,units="secs"),2)," seconds"));
cat(paste0("\nModel estimated: ",x$model));
if(is.occurrence(x$occurrence)){
occurrence <- switch(x$occurrence$occurrence,
"f"=,
"fixed"="Fixed probability",
"o"=,
"odds-ratio"="Odds ratio",
"i"=,
"inverse-odds-ratio"="Inverse odds ratio",
"d"=,
"direct"="Direct",
"g"=,
"general"="General",
"p"=,
"provided"="Provided by user");
cat(paste0("\nOccurrence model type: ",occurrence));
}
distrib <- switch(x$distribution,
"dnorm" = "Normal",
"dlaplace" = "Laplace",
"ds" = "S",
"dgnorm" = paste0("Generalised Normal with shape=",round(x$lambda, digits)),
"dlogis" = "Logistic",
"dt" = paste0("Student t with df=",round(x$lambda, digits)),
"dalaplace" = paste0("Asymmetric Laplace with lambda=",round(x$lambda,digits)),
"dlnorm" = "Log Normal",
"dllaplace" = "Log Laplace",
"dls" = "Log S",
"dlgnorm" = paste0("Log Generalised Normal with shape=",round(x$lambda, digits)),
# "dbcnorm" = paste0("Box-Cox Normal with lambda=",round(x$other$lambda,2)),
"dinvgauss" = "Inverse Gaussian"
);
if(is.occurrence(x$occurrence)){
distrib <- paste0("Mixture of Bernoulli and ", distrib);
}
cat(paste0("\nDistribution assumed in the model: ", distrib));
cat(paste0("\nLoss function type: ",x$loss));
if(!is.null(x$lossValue)){
cat(paste0("; Loss function value: ",round(x$lossValue,digits)));
if(any(x$loss==c("LASSO","RIDGE"))){
cat(paste0("; lambda=",x$lambda));
}
}
if(etsModel){
if(!is.null(x$persistence)){
cat(paste0("\nPersistence vector g"));
if(!is.null(x$xreg)){
cat(" (excluding xreg):\n");
}
else{
cat(":\n");
}
persistence <- x$persistence[substr(names(x$persistence),1,5)!="delta"];
if(arimaModel){
persistence <- persistence[substr(names(persistence),1,3)!="psi"];
}
print(round(persistence,digits));
}
if(!is.null(x$phi)){
if(gregexpr("d",x$model)!=-1){
cat(paste0("Damping parameter: ", round(x$phi,digits)));
}
}
}
# If this is ARIMA model
if(!is.null(x$arma) && (!is.null(x$arma$ar) || !is.null(x$arma$ma))){
cat(paste0("\nARMA parameters of the model:\n"));
if(!is.null(x$arma$ar)){
cat("AR:\n")
print(round(x$arma$ar,digits));
}
if(!is.null(x$arma$ma)){
cat("MA:\n")
print(round(x$arma$ma,digits));
}
}
cat("\nSample size: "); cat(nobs(x));
cat("\nNumber of estimated parameters: "); cat(nparam(x));
cat("\nNumber of degrees of freedom: "); cat(nobs(x)-nparam(x));
if(x$nParam[2,4]>0){
cat("\nNumber of provided parameters: "); cat(x$nParam[2,4]);
}
if(x$loss=="likelihood" ||
(any(x$loss==c("MSE","MSEh","MSCE","GPL")) & any(x$distribution==c("dnorm","dlnorm"))) ||
(any(x$loss==c("aMSE","aMSEh","aMSCE","aGPL")) & any(x$distribution==c("dnorm","dlnorm"))) ||
(any(x$loss==c("MAE","MAEh","MACE")) & any(x$distribution==c("dlaplace","dllaplace"))) ||
(any(x$loss==c("HAM","HAMh","CHAM")) & any(x$distribution==c("ds","dls")))){
ICs <- c(AIC(x),AICc(x),BIC(x),BICc(x));
names(ICs) <- c("AIC","AICc","BIC","BICc");
cat("\nInformation criteria:\n");
print(round(ICs,digits));
}
else{
cat("\nInformation criteria are unavailable for the chosen loss & distribution.\n");
}
# If there are accuracy measures, print them out
if(!is.null(x$accuracy)){
cat("\nForecast errors:\n");
if(is.null(x$occurrence)){
cat(paste(paste0("ME: ",round(x$accuracy["ME"],3)),
paste0("MAE: ",round(x$accuracy["MAE"],3)),
paste0("RMSE: ",round(sqrt(x$accuracy["MSE"]),3),"\n")
# paste0("Bias: ",round(x$accuracy["cbias"],3)*100,"%"),
,sep="; "));
cat(paste(paste0("sCE: ",round(x$accuracy["sCE"],5)*100,"%"),
paste0("sMAE: ",round(x$accuracy["sMAE"],5)*100,"%"),
paste0("sMSE: ",round(x$accuracy["sMSE"],5)*100,"%\n")
,sep="; "));
cat(paste(paste0("MASE: ",round(x$accuracy["MASE"],3)),
paste0("RMSSE: ",round(x$accuracy["RMSSE"],3)),
paste0("rMAE: ",round(x$accuracy["rMAE"],3)),
paste0("rRMSE: ",round(x$accuracy["rRMSE"],3),"\n")
,sep="; "));
}
else{
cat(paste(paste0("Bias: ",round(x$accuracy["cbias"],5)*100,"%"),
paste0("sMSE: ",round(x$accuracy["sMSE"],5)*100,"%"),
paste0("rRMSE: ",round(x$accuracy["rRMSE"],3)),
paste0("sPIS: ",round(x$accuracy["sPIS"],5)*100,"%"),
paste0("sCE: ",round(x$accuracy["sCE"],5)*100,"%\n"),sep="; "));
}
}
}
#' @export
print.adamCombined <- function(x, digits=4, ...){
cat(paste0("Time elapsed: ",round(as.numeric(x$timeElapsed,units="secs"),2)," seconds"));
cat(paste0("\nModel estimated: ",x$model));
cat(paste0("\nLoss function type: ",x$models[[1]]$loss));
cat(paste0("\n\nNumber of models combined: ", length(x$ICw)));
cat("\nSample size: "); cat(nobs(x));
cat("\nAverage number of estimated parameters: "); cat(round(nparam(x),digits=digits));
cat("\nAverage number of degrees of freedom: "); cat(round(nobs(x)-nparam(x),digits=digits));
if(!is.null(x$accuracy)){
cat("\n\nForecast errors:\n");
if(is.null(x$occurrence)){
cat(paste(paste0("ME: ",round(x$accuracy["ME"],3)),
paste0("MAE: ",round(x$accuracy["MAE"],3)),
paste0("RMSE: ",round(sqrt(x$accuracy["MSE"]),3),"\n")
# paste0("Bias: ",round(x$accuracy["cbias"],3)*100,"%"),
,sep="; "));
cat(paste(paste0("sCE: ",round(x$accuracy["sCE"],5)*100,"%"),
paste0("sMAE: ",round(x$accuracy["sMAE"],5)*100,"%"),
paste0("sMSE: ",round(x$accuracy["sMSE"],5)*100,"%\n")
,sep="; "));
cat(paste(paste0("MASE: ",round(x$accuracy["MASE"],3)),
paste0("RMSSE: ",round(x$accuracy["RMSSE"],3)),
paste0("rMAE: ",round(x$accuracy["rMAE"],3)),
paste0("rRMSE: ",round(x$accuracy["rRMSE"],3),"\n")
,sep="; "));
}
else{
cat(paste(paste0("Bias: ",round(x$accuracy["cbias"],5)*100,"%"),
paste0("sMSE: ",round(x$accuracy["sMSE"],5)*100,"%"),
paste0("rRMSE: ",round(x$accuracy["rRMSE"],3)),
paste0("sPIS: ",round(x$accuracy["sPIS"],5)*100,"%"),
paste0("sCE: ",round(x$accuracy["sCE"],5)*100,"%\n"),sep="; "));
}
}
}
#### Coefficients ####
confint.adam <- function(object, parm, level=0.95, ...){
adamVcov <- vcov(object);
adamSD <- sqrt(abs(diag(adamVcov)));
# adamCoef <- coef(object);
adamCoefBounds <- matrix(0,length(adamSD),2);
adamCoefBounds[,1] <- qnorm((1-level)/2, 0, adamSD);
adamCoefBounds[,2] <- qnorm((1+level)/2, 0, adamSD);
adamReturn <- cbind(adamSD,adamCoefBounds);
colnames(adamReturn) <- c("S.E.",
paste0((1-level)/2*100,"%"), paste0((1+level)/2*100,"%"));
return(adamReturn);
}
#' @export
coef.adam <- function(object, ...){
return(object$B);
}
#' @importFrom stats sigma
#' @export
sigma.adam <- function(object, ...){
df <- (nobs(object, all=FALSE)-nparam(object));
# If the sample is too small, then use biased estimator
if(df<=0){
df[] <- nobs(object);
}
return(sqrt(switch(object$distribution,
"dnorm"=,
"dlaplace"=,
"ds"=,
"dgnorm"=,
"dt"=,
"dlogis"=,
"dalaplace"=sum(residuals(object)^2),
"dlnorm"=,
"dllaplace"=,
"dls"=,
"dlgnorm"=sum(log(residuals(object))^2),
"dinvgauss"=sum((residuals(object)-1)^2))
/df));
}
#' @export
summary.adam <- function(object, level=0.95, ...){
ourReturn <- list(model=object$model,responseName=all.vars(formula(object))[1]);
occurrence <- NULL;
if(is.occurrence(object$occurrence)){
occurrence <- switch(object$occurrence$occurrence,
"f"=,
"fixed"="Fixed probability",
"o"=,
"odds-ratio"="Odds ratio",
"i"=,
"inverse-odds-ratio"="Inverse odds ratio",
"d"=,
"direct"="Direct",
"g"=,
"general"="General");
}
ourReturn$occurrence <- occurrence;
ourReturn$distribution <- object$distribution;
# Collect parameters and their standard errors
parametersConfint <- confint(object, level=level);
parametersValues <- coef(object);
if(is.null(parametersValues)){
if(!is.null(object$xreg) && all(object$persistenceXreg!=0)){
parametersValues <- c(object$persistence,object$persistenceXreg,object$initial,object$initialXreg);
}
else{
parametersValues <- c(object$persistence,object$initial);
}
warning(paste0("Parameters are not available. You have probably provided them in the model, ",
"so there was nothing to estimate. We extracted smoothing parameters and initials."),
call.=FALSE);
}
parametersConfint[,2:3] <- parametersValues + parametersConfint[,2:3];
parametersTable <- cbind(parametersValues,parametersConfint);
rownames(parametersTable) <- rownames(parametersConfint);
colnames(parametersTable) <- c("Estimate","Std. Error",
paste0("Lower ",(1-level)/2*100,"%"),
paste0("Upper ",(1+level)/2*100,"%"));
ourReturn$coefficients <- parametersTable;
ourReturn$loss <- object$loss;
ourReturn$lossValue <- object$lossValue;
ourReturn$nobs <- nobs(object);
ourReturn$nparam <- nparam(object);
ourReturn$nParam <- object$nParam;
if(object$loss=="likelihood" ||
(any(object$loss==c("MSE","MSEh","MSCE")) & any(object$distribution==c("dnorm","dlnorm"))) ||
(any(object$loss==c("MAE","MAEh","MACE")) & any(object$distribution==c("dlaplace","dllaplace"))) ||
(any(object$loss==c("HAM","HAMh","CHAM")) & any(object$distribution==c("ds","dls")))){
ICs <- c(AIC(object),AICc(object),BIC(object),BICc(object));
names(ICs) <- c("AIC","AICc","BIC","BICc");
ourReturn$ICs <- ICs;
}
return(structure(ourReturn, class="summary.adam"));
}
#' @export
summary.adamCombined <- function(object, ...){
return(print.adamCombined(object, ...));
}
#' @export
print.summary.adam <- function(x, ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="digits")){
digits <- 4;
}
else{
digits <- ellipsis$digits;
}
cat(paste0("Model estimated: ",x$model));
cat(paste0("\nResponse variable: ", paste0(x$responseName,collapse="")));
if(!is.null(x$occurrence)){
cat(paste0("\nOccurrence model type: ",x$occurrence));
}
distrib <- switch(x$distribution,
"dnorm" = "Normal",
"dlaplace" = "Laplace",
"ds" = "S",
"dgnorm" = "Generalised Normal",
"dlogis" = "Logistic",
"dt" = paste0("Student t with df=",round(x$lambda, digits)),
"dalaplace" = paste0("Asymmetric Laplace with lambda=",round(x$lambda,digits)),
"dlnorm" = "Log Normal",
"dllaplace" = "Log Laplace",
"dls" = "Log S",
"dlgnorm" = "Log Generalised Normal",
# "dbcnorm" = paste0("Box-Cox Normal with lambda=",round(x$other$lambda,2)),
"dinvgauss" = "Inverse Gaussian"
);
if(!is.null(x$occurrence)){
distrib <- paste0("\nMixture of Bernoulli and ", distrib);
}
cat(paste0("\nDistribution used in the estimation: ", distrib));
cat(paste0("\nLoss function type: ",x$loss));
if(!is.null(x$lossValue)){
cat(paste0("; Loss function value: ",round(x$lossValue,digits)));
if(any(x$loss==c("LASSO","RIDGE"))){
cat(paste0("; lambda=",x$lambda));
}
}
if(!is.null(x$coefficients)){
cat("\nCoefficients:\n");
print(round(x$coefficients,digits));
}
cat("\nSample size: "); cat(x$nobs);
cat("\nNumber of estimated parameters: "); cat(x$nparam);
cat("\nNumber of degrees of freedom: "); cat(x$nobs-x$nparam);
if(x$nParam[2,4]>0){
cat("\nNumber of provided parameters: "); cat(x$nParam[2,4]);
}
if(x$loss=="likelihood" ||
(any(x$loss==c("MSE","MSEh","MSCE")) & any(x$distribution==c("dnorm","dlnorm"))) ||
(any(x$loss==c("MAE","MAEh","MACE")) & any(x$distribution==c("dlaplace","dllaplace"))) ||
(any(x$loss==c("HAM","HAMh","CHAM")) & any(x$distribution==c("ds","dls")))){
cat("\nInformation criteria:\n");
print(round(x$ICs,digits));
}
else{
cat("\nInformation criteria are unavailable for the chosen loss & distribution.\n");
}
}
#' @export
vcov.adam <- function(object, ...){
# If the forecast is in numbers, then use its length as a horizon
if(any(!is.na(object$forecast))){
h <- length(object$forecast)
}
else{
h <- 0;
}
y <- actuals(object);
modelReturn <- suppressWarnings(adam(y, h=0, model=object, FI=TRUE));
vcovMatrix <- try(chol2inv(chol(modelReturn$FI)), silent=TRUE);
if(inherits(vcovMatrix,"try-error")){
vcovMatrix <- try(solve(modelReturn$FI, diag(ncol(modelReturn$FI)), tol=1e-20), silent=TRUE);
if(inherits(vcovMatrix,"try-error")){
warning(paste0("Sorry, but the hessian is singular, so we could not invert it.\n",
"We failed to produce the covariance matrix of parameters."),
call.=FALSE);
vcovMatrix <- diag(1e+100,ncol(modelReturn$FI));
}
}
colnames(vcovMatrix) <- rownames(vcovMatrix) <- colnames(modelReturn$FI);
return(vcovMatrix);
}
#### Residuals and actuals functions ####
#' @importFrom greybox actuals
#' @export
actuals.adam <- function(object, all=TRUE, ...){
if(all){
return(object$y);
}
else{
return(object$y[object$y!=0]);
}
}
#' @export
nobs.adam <- function(object, ...){
return(length(actuals(object, ...)));
}
#' @export
residuals.adam <- function(object, ...){
return(switch(object$distribution,
"dlnorm"=,
"dllaplace"=,
"dls"=,
"dlgnorm"=,
"dinvgauss"=switch(errorType(object),
# abs() is needed in order to avoid the weird cases
"A"=abs(1+object$residuals/fitted(object)),
"M"=1+object$residuals),
"dnorm"=,
"dlaplace"=,
"ds"=,
"dgnorm"=,
"dlogis"=,
"dt"=,
"dalaplace"=,
object$residuals));
}
#' Multiple steps ahead forecast errors
#'
#' The function extracts 1 to h steps ahead forecast errors from the model.
#'
#' 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.
#'
#' @template ssAuthor
#' @template ssKeywords
#'
#' @param object Model estimated using one of the forecasting functions.
#' @param h The forecasting horizon to use.
#' @param ... Currently nothing is accepted via ellipsis.
#' @return 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.
#' @seealso \link[stats]{residuals}, \link[stats]{rstandard}, \link[stats]{rstudent}
#' @examples
#'
#' x <- rnorm(100,0,1)
#' ourModel <- adam(x)
#' rmultistep(ourModel, h=13)
#'
#' @export rmultistep
rmultistep <- function(object, h=10, ...) UseMethod("rmultistep")
#' @export
rmultistep.default <- function(object, h=10, ...){
return(NULL);
}
#' @export
rmultistep.adam <- function(object, h=10, ...){
yClasses <- class(actuals(object));
# Model type
model <- modelType(object);
Etype <- errorType(object);
Ttype <- substr(model,2,2);
Stype <- substr(model,nchar(model),nchar(model));
# Technical parameters
lagsModelAll <- modelLags(object);
lagsModelMax <- max(lagsModelAll);
lagsOriginal <- lags(object);
if(Ttype!="N"){
lagsOriginal <- c(1,lagsOriginal);
}
if(!is.null(object$initial$seasonal)){
if(is.list(object$initial$seasonal)){
componentsNumberETSSeasonal <- length(object$initial$seasonal);
}
else{
componentsNumberETSSeasonal <- 1;
}
}
else{
componentsNumberETSSeasonal <- 0;
}
componentsNumberETS <- length(object$initial$level) + length(object$initial$trend) + componentsNumberETSSeasonal;
componentsNumberARIMA <- sum(substr(colnames(object$states),1,10)=="ARIMAState");
if(!is.null(object$xreg)){
xregNumber <- ncol(object$xreg);
}
else{
xregNumber <- 0;
}
obsInSample <- nobs(object);
# Function returns the matrix with multi-step errors
if(is.occurrence(object$occurrence)){
ot <- matrix(actuals(object$occurrence),obsInSample,1);
}
else{
ot <- matrix(1,obsInSample,1);
}
# Return multi-step errors matrix
if(any(yClasses=="ts")){
return(ts(adamErrorerWrap(t(object$states), object$measurement, object$transition,
lagsModelAll, Etype, Ttype, Stype,
componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, h,
matrix(actuals(object),obsInSample,1), ot),
start=start(actuals(object)), frequency=frequency(actuals(object))));
}
else{
return(zoo(adamErrorerWrap(t(object$states), object$measurement, object$transition,
lagsModelAll, Etype, Ttype, Stype,
componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber, h,
matrix(actuals(object),obsInSample,1), ot),
order.by=time(actuals(object))));
}
}
#' @importFrom stats rstandard
#' @export
rstandard.adam <- function(model, ...){
obs <- nobs(model);
df <- obs - nparam(model);
errors <- residuals(model);
# If this is an occurrence model, then only modify the non-zero obs
# Also, if there are NAs in actuals, consider them as occurrence
if(is.occurrence(model$occurrence)){
residsToGo <- which(actuals(model$occurrence)!=0 & !is.na(actuals(model)));
}
else{
residsToGo <- c(1:obs);
}
if(any(model$distribution==c("dt","dnorm"))){
return((errors - mean(errors[residsToGo])) / sqrt(model$scale^2 * obs / df));
}
else if(model$distribution=="ds"){
return((errors - mean(errors[residsToGo])) / (model$scale * obs / df)^2);
}
else if(model$distribution=="dls"){
errors[] <- log(errors);
return(exp((errors - mean(errors[residsToGo])) / (model$scale * obs / df)^2));
}
else if(model$distribution=="dgnorm"){
return((errors - mean(errors[residsToGo])) / (model$scale^model$lambda * obs / df)^{1/model$lambda});
}
else if(model$distribution=="dlgnorm"){
errors[] <- log(errors);
return(exp((errors - mean(errors[residsToGo])) / (model$scale^model$lambda * obs / df)^{1/model$lambda}));
}
else if(model$distribution=="dinvgauss"){
return(errors / mean(errors[residsToGo]));
}
else if(model$distribution=="dlnorm"){
errors[] <- log(errors);
return(exp((errors - mean(errors[residsToGo])) / sqrt(model$scale^2 * obs / df)));
}
else if(model$distribution=="dllaplace"){
errors[] <- log(errors);
return(exp((errors - mean(errors[residsToGo])) / model$scale * obs / df));
}
else{
return(errors / model$scale * obs / df);
}
}
#' @importFrom stats rstudent
#' @export
rstudent.adam <- function(model, ...){
obs <- nobs(model);
df <- obs - nparam(model) - 1;
rstudentised <- errors <- residuals(model);
# If this is an occurrence model, then only modify the non-zero obs
# Also, if there are NAs in actuals, consider them as occurrence
if(is.occurrence(model$occurrence)){
residsToGo <- which(actuals(model$occurrence)!=0 & !is.na(actuals(model)));
}
else{
residsToGo <- c(1:obs);
}
if(any(model$distribution==c("dt","dnorm"))){
errors[] <- errors - mean(errors);
for(i in residsToGo){
rstudentised[i] <- errors[i] / sqrt(sum(errors[-i]^2,na.rm=TRUE) / df);
}
}
else if(model$distribution=="dlaplace"){
errors[] <- errors - mean(errors);
for(i in residsToGo){
rstudentised[i] <- errors[i] / (sum(abs(errors[-i]),na.rm=TRUE) / df);
}
}
else if(model$distribution=="dlnorm"){
errors[] <- log(errors) - mean(log(errors));
for(i in residsToGo){
rstudentised[i] <- exp(errors[i] / sqrt(sum(errors[-i]^2,na.rm=TRUE) / df));
}
}
else if(model$distribution=="dllaplace"){
errors[] <- log(errors) - mean(log(errors));
for(i in residsToGo){
rstudentised[i] <- exp(errors[i] / (sum(abs(errors[-i]),na.rm=TRUE) / df));
}
}
else if(model$distribution=="ds"){
errors[] <- errors - mean(errors);
for(i in residsToGo){
rstudentised[i] <- errors[i] / (sum(sqrt(abs(errors[-i])),na.rm=TRUE) / (2*df))^2;
}
}
else if(model$distribution=="dls"){
errors[] <- log(errors) - mean(log(errors));
for(i in residsToGo){
rstudentised[i] <- exp(errors[i] / (sum(sqrt(abs(errors[-i])),na.rm=TRUE) / (2*df))^2);
}
}
else if(model$distribution=="dgnorm"){
errors[] <- errors - mean(errors);
for(i in residsToGo){
rstudentised[i] <- errors[i] / (sum(abs(errors[-i])^model$lambda) * (model$lambda/df))^{1/model$lambda};
}
}
else if(model$distribution=="dlgnorm"){
errors[] <- log(errors) - mean(log(errors));
for(i in residsToGo){
rstudentised[i] <- errors[i] / (sum(abs(errors[-i])^model$lambda) * (model$lambda/df))^{1/model$lambda};
}
}
else if(model$distribution=="dalaplace"){
for(i in residsToGo){
rstudentised[i] <- errors[i] / (sum(errors[-i] * (model$lambda - (errors[-i]<=0)*1),na.rm=TRUE) / df);
}
}
else if(model$distribution=="dlogis"){
errors[] <- errors - mean(errors);
for(i in residsToGo){
rstudentised[i] <- errors[i] / (sqrt(sum(errors[-i]^2,na.rm=TRUE) / df) * sqrt(3) / pi);
}
}
else if(model$distribution=="dinvgauss"){
for(i in residsToGo){
rstudentised[i] <- errors[i] / mean(errors[residsToGo][-i],na.rm=TRUE);
}
}
else{
for(i in residsToGo){
rstudentised[i] <- errors[i] / sqrt(sum(errors[-i]^2,na.rm=TRUE) / df);
}
}
return(rstudentised);
}
#' @importFrom greybox outlierdummy
#' @export
outlierdummy.adam <- function(object, level=0.999, type=c("rstandard","rstudent"), ...){
# Function returns the matrix of dummies with outliers
type <- match.arg(type);
errors <- switch(type,"rstandard"=rstandard(object),"rstudent"=rstudent(object));
statistic <- switch(object$distribution,
"dlaplace"=,
"dllaplace"=qlaplace(c((1-level)/2, (1+level)/2), 0, 1),
"dalaplace"=qalaplace(c((1-level)/2, (1+level)/2), 0, 1, object$other$alpha),
"dlogis"=qlogis(c((1-level)/2, (1+level)/2), 0, 1),
"dt"=qt(c((1-level)/2, (1+level)/2), nobs(object)-nparam(object)),
"dgnorm"=,
"dlgnorm"=qgnorm(c((1-level)/2, (1+level)/2), 0, 1, object$other$beta),
"ds"=,
"dls"=qs(c((1-level)/2, (1+level)/2), 0, 1),
# In the next one, the scale is debiased, taking n-k into account
"dinvgauss"=qinvgauss(c((1-level)/2, (1+level)/2), mean=1,
dispersion=object$scale * nobs(object) /
(nobs(object)-nparam(object))),
qnorm(c((1-level)/2, (1+level)/2), 0, 1));
if(any(object$distribution==c("dlnorm","dllaplace","dls","dlgnorm"))){
errors[] <- log(errors);
}
outliersID <- which(errors>statistic[2] | errors<statistic[1]);
outliersNumber <- length(outliersID);
if(outliersNumber>0){
outliers <- matrix(0, nobs(object), outliersNumber,
dimnames=list(rownames(object$data),
paste0("outlier",c(1:outliersNumber))));
outliers[cbind(outliersID,c(1:outliersNumber))] <- 1;
}
else{
outliers <- NULL;
}
return(structure(list(outliers=outliers, statistic=statistic, id=outliersID,
level=level, type=type),
class="outlierdummy"));
}
#### Predict and forecast functions ####
#' @export
predict.adam <- function(object, newxreg=NULL, interval=c("none", "confidence", "prediction"),
level=0.95, side=c("both","upper","lower"), ...){
interval <- match.arg(interval);
obsInSample <- nobs(object);
# Check if newxreg is provided
if(!is.null(newxreg)){
# If this is not a matrix / data.frame, then convert to one
if(!is.data.frame(newxreg) && !is.matrix(newxreg)){
newxreg <- as.data.frame(newxreg);
colnames(newxreg) <- "xreg";
}
h <- nrow(newxreg);
# If the newxreg is provided, then just do forecasts for that part
if(any(interval==c("none","prediction","confidence"))){
if(interval==c("prediction")){
interval[] <- "simulated";
}
return(forecast(object, h=h, newxreg=newxreg,
interval=interval,
level=level, side=side, ...));
}
}
else{
# If there are no newxreg, then we need to produce fitted with / without interval
if(interval=="none"){
return(structure(list(mean=fitted(object), lower=NA, upper=NA, model=object,
level=level, interval=interval, side=side),
class=c("adam.predict","adam.forecast")));
}
# Otherwise we do one-step-ahead prediction / confidence interval
else{
yForecast <- fitted(object);
}
}
##### Prediction interval for in sample only! #####
side <- match.arg(side);
if(length(level)>1){
warning(paste0("Sorry, but we only support scalar for the level, ",
"when constructing in-sample prediction interval. ",
"Using the first provided value."),
call.=FALSE);
level <- level[1];
}
# Fix just in case a silly user used 95 etc instead of 0.95
if(level>1){
level[] <- level / 100;
}
# Basic parameters
model <- modelType(object);
Etype <- errorType(object);
# ts structure
yForecastStart <- time(actuals(object))[obsInSample]+deltat(actuals(object));
yFrequency <- frequency(actuals(object));
# Extract variance and amend it in case of confidence interval
s2 <- sigma(object)^2;
if(interval=="confidence"){
warning(paste0("Note that the ETS assumes that the initial level is known, ",
"so the confidence interval depends on smoothing parameters only."),
call.=FALSE);
s2 <- s2 * object$measurement[1:obsInSample,1:length(object$persistence),drop=FALSE] %*% object$persistence;
}
yUpper <- yLower <- yForecast;
yUpper[] <- yLower[] <- NA;
# If this is a mixture model, produce forecasts for the occurrence
if(!is.null(object$occurrence)){
occurrenceModel <- TRUE;
pForecast <- fitted(object$occurrence);
}
else{
occurrenceModel <- FALSE;
pForecast <- rep(1, obsInSample);
}
# If this is an occurrence model, then take probability into account in the level.
if(occurrenceModel && (interval=="prediction")){
levelNew <- (level-(1-pForecast))/pForecast;
levelNew[levelNew<0] <- 0;
}
else{
levelNew <- level;
}
levelLow <- levelUp <- vector("numeric",obsInSample);
if(side=="both"){
levelLow[] <- (1-levelNew)/2;
levelUp[] <- (1+levelNew)/2;
}
else if(side=="upper"){
levelLow[] <- rep(0,length(levelNew));
levelUp[] <- levelNew;
}
else{
levelLow[] <- 1-levelNew;
levelUp[] <- rep(1,length(levelNew));
}
levelLow[levelLow<0] <- 0;
levelUp[levelUp<0] <- 0;
#### Produce the intervals for the data ####
if(object$distribution=="dnorm"){
if(Etype=="A"){
yLower[] <- qnorm(levelLow, 0, sqrt(s2));
yUpper[] <- qnorm(levelUp, 0, sqrt(s2));
}
else{
yLower[] <- qnorm(levelLow, 1, sqrt(s2));
yUpper[] <- qnorm(levelUp, 1, sqrt(s2));
}
}
else if(object$distribution=="dlaplace"){
if(Etype=="A"){
yLower[] <- qlaplace(levelLow, 0, sqrt(s2/2));
yUpper[] <- qlaplace(levelUp, 0, sqrt(s2/2));
}
else{
yLower[] <- qlaplace(levelLow, 1, sqrt(s2/2));
yUpper[] <- qlaplace(levelUp, 1, sqrt(s2/2));
}
}
else if(object$distribution=="ds"){
if(Etype=="A"){
yLower[] <- qs(levelLow, 0, (s2/120)^0.25);
yUpper[] <- qs(levelUp, 0, (s2/120)^0.25);
}
else{
yLower[] <- qs(levelLow, 1, (s2/120)^0.25);
yUpper[] <- qs(levelUp, 1, (s2/120)^0.25);
}
}
else if(object$distribution=="dgnorm"){
alpha <- sqrt(s2*(gamma(1/object$lambda)/gamma(3/object$lambda)));
if(Etype=="A"){
yLower[] <- suppressWarnings(qgnorm(levelLow, 0, alpha, object$lambda));
yUpper[] <- suppressWarnings(qgnorm(levelUp, 0, alpha, object$lambda));
}
else{
yLower[] <- suppressWarnings(qgnorm(levelLow, 1, alpha, object$lambda));
yUpper[] <- suppressWarnings(qgnorm(levelUp, 1, alpha, object$lambda));
}
}
else if(object$distribution=="dlogis"){
if(Etype=="A"){
yLower[] <- qlogis(levelLow, 0, sqrt(s2*3)/pi);
yUpper[] <- qlogis(levelUp, 0, sqrt(s2*3)/pi);
}
else{
yLower[] <- qlogis(levelLow, 1, sqrt(s2*3)/pi);
yUpper[] <- qlogis(levelUp, 1, sqrt(s2*3)/pi);
}
}
else if(object$distribution=="dt"){
df <- nobs(object) - nparam(object);
if(Etype=="A"){
yLower[] <- sqrt(s2)*qt(levelLow, df);
yUpper[] <- sqrt(s2)*qt(levelUp, df);
}
else{
yLower[] <- (1 + sqrt(s2)*qt(levelLow, df));
yUpper[] <- (1 + sqrt(s2)*qt(levelUp, df));
}
}
else if(object$distribution=="dalaplace"){
lambda <- object$lambda;
if(Etype=="A"){
yLower[] <- qalaplace(levelLow, 0,
sqrt(s2*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
yUpper[] <- qalaplace(levelUp, 0,
sqrt(s2*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
}
else{
yLower[] <- qalaplace(levelLow, 1,
sqrt(s2*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
yUpper[] <- qalaplace(levelUp, 1,
sqrt(s2*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
}
}
else if(object$distribution=="dlnorm"){
yLower[] <- qlnorm(levelLow, 0, sqrt(s2));
yUpper[] <- qlnorm(levelUp, 0, sqrt(s2));
}
else if(object$distribution=="dllaplace"){
yLower[] <- exp(qlaplace(levelLow, 0, sqrt(s2/2)));
yUpper[] <- exp(qlaplace(levelUp, 0, sqrt(s2/2)));
}
else if(object$distribution=="dls"){
yLower[] <- exp(qs(levelLow, 0, (s2/120)^0.25));
yUpper[] <- exp(qs(levelUp, 0, (s2/120)^0.25));
}
else if(object$distribution=="dlgnorm"){
alpha <- sqrt(s2*(gamma(1/object$lambda)/gamma(3/object$lambda)));
yLower[] <- suppressWarnings(exp(qgnorm(levelLow, 0, alpha, object$lambda)));
yUpper[] <- suppressWarnings(exp(qgnorm(levelUp, 0, alpha, object$lambda)));
}
else if(object$distribution=="dinvgauss"){
yLower[] <- qinvgauss(levelLow, 1, dispersion=s2);
yUpper[] <- qinvgauss(levelUp, 1, dispersion=s2);
}
#### Clean up the produced values for the interval ####
# Make sensible values out of those weird quantiles
if(Etype=="A"){
yLower[levelLow==0] <- -Inf;
}
else{
yLower[levelLow==0] <- 0;
}
yUpper[levelUp==1] <- Inf;
# Substitute NAs and NaNs with zeroes
if(any(is.nan(yLower)) || any(is.na(yLower))){
yLower[is.nan(yLower)] <- 0;
yLower[is.na(yLower)] <- 0;
}
if(any(is.nan(yUpper)) || any(is.na(yUpper))){
yUpper[is.nan(yUpper)] <- 0;
yUpper[is.na(yUpper)] <- 0;
}
if(Etype=="A"){
yLower[] <- yForecast + yLower;
yUpper[] <- yForecast + yUpper;
}
else{
yLower[] <- yForecast * yLower;
yUpper[] <- yForecast * yUpper;
}
return(structure(list(mean=yForecast, lower=yLower, upper=yUpper, model=object,
level=level, interval=interval, side=side),
class=c("adam.predict","adam.forecast")));
}
#' @export
plot.adam.predict <- function(x, ...){
ellipsis <- list(...);
if(is.null(ellipsis$ylim)){
ellipsis$ylim <- range(c(actuals(x$model),x$mean,x$lower,x$upper),na.rm=TRUE);
}
ellipsis$x <- actuals(x$model);
do.call(plot, ellipsis);
lines(x$mean,col="purple",lwd=2,lty=2);
if(x$interval!="none"){
lines(x$lower,col="grey",lwd=3,lty=2);
lines(x$upper,col="grey",lwd=3,lty=2);
}
}
# Work in progress...
#' @param nsim Number of iterations to do in case of \code{interval="simulated"}.
#' @param interval What type of mechanism to use for interval construction. The
#' most statistically correct one is \code{interval="simulated"} (this is
#' recommended method), but it is the slowest method. \code{interval="approximate"}
#' (aka \code{interval="parametric"}) uses analytical formulae for conditional
#' h-steps ahead variance, but is approximate for the non-additive error models.
#' \code{interval="semiparametric"} relies on the multiple steps ahead forecast
#' error and an assumed distribution of the error term. \code{interval="nonparametric"}
#' uses Taylor & Bunn (1999) approach with quantile regressions. Finally,
#' \code{interval="confidence"} tries to generate the confidence intervals for the
#' point forecast. This relies on the assumption that the parameters of ETS are known
#' (unrealistic, but a standard one). The function also accepts
#' \code{interval="parametric"} and \code{interval="prediction"}, which are equivalent
#' to \code{interval="approximate"}.
#' @param occurrence The vector containing the future occurrence variable
#' (values in [0,1]), if it is known.
#' @rdname forecast.smooth
#' @importFrom stats rnorm rlogis rt rlnorm qnorm qlogis qt qlnorm
#' @importFrom statmod rinvgauss qinvgauss
#' @importFrom greybox rlaplace rs ralaplace qlaplace qs qalaplace
#' @export
forecast.adam <- function(object, h=10, newxreg=NULL, occurrence=NULL,
interval=c("none", "simulated", "approximate", "semiparametric", "nonparametric", "confidence"),
level=0.95, side=c("both","upper","lower"), cumulative=FALSE, nsim=10000, ...){
ellipsis <- list(...);
interval <- match.arg(interval[1],c("none", "simulated", "approximate", "semiparametric",
"nonparametric", "confidence", "parametric","prediction"));
# If the horizon is zero, just construct fitted and potentially confidence interval thingy
if(h<=0){
if(all(interval!=c("none","confidence"))){
interval[] <- "prediction";
}
return(predict(object, newxreg=newxreg,
interval=interval,
level=level, side=side, ...));
}
if(any(interval==c("parametric","prediction"))){
# warning("The parameter 'interval' does not accept 'parametric' anymore. We use 'approximate' value instead.",
# call.=FALSE);
interval <- "approximate";
}
else if(interval=="confidence"){
warning(paste0("Note that the ETS assumes that the initial level is known, ",
"so the confidence interval depends on smoothing parameters only."),
call.=FALSE);
}
side <- match.arg(side);
# Model type
model <- modelType(object);
Etype <- errorType(object);
Ttype <- substr(model,2,2);
Stype <- substr(model,nchar(model),nchar(model));
# Technical parameters
lagsModelAll <- modelLags(object);
lagsModelMax <- max(lagsModelAll);
if(!is.null(object$initial$seasonal)){
if(is.list(object$initial$seasonal)){
componentsNumberETSSeasonal <- length(object$initial$seasonal);
}
else{
componentsNumberETSSeasonal <- 1;
}
}
else{
componentsNumberETSSeasonal <- 0;
}
componentsNumberETS <- length(object$initial$level) + length(object$initial$trend) + componentsNumberETSSeasonal;
componentsNumberARIMA <- sum(substr(colnames(object$states),1,10)=="ARIMAState");
obsStates <- nrow(object$states);
obsInSample <- nobs(object);
yClasses <- class(actuals(object));
if(any(yClasses=="ts")){
# ts structure
yForecastStart <- time(actuals(object))[obsInSample]+deltat(actuals(object));
yFrequency <- frequency(actuals(object));
}
else{
# zoo thingy
yIndex <- time(actuals(object));
yForecastIndex <- yIndex[obsInSample]+diff(tail(yIndex,2))*c(1:h);
}
# All the important matrices
matVt <- t(object$states[obsStates-(lagsModelMax:1)+1,,drop=FALSE]);
matWt <- tail(object$measurement,h);
# If the forecast horizon is higher than the in-sample, duplicate the last value in matWt
if(nrow(matWt)<h){
matWt <- matrix(tail(matWt,1), nrow=h, ncol=ncol(matWt), dimnames=list(NULL,colnames(matWt)), byrow=TRUE);
}
vecG <- matrix(object$persistence, ncol=1);
# Deal with explanatory variables
if(!is.null(object$xreg)){
xregNumber <- ncol(object$xreg);
if(is.null(newxreg) && (nrow(object$xreg)<(obsInSample+h))){
warning("The newxreg is not provided. Predicting the explanatory variables based on what we have in-sample.",
call.=FALSE);
newxreg <- matrix(NA,h,xregNumber);
for(i in 1:xregNumber){
newxreg[,i] <- adam(object$xreg[,i],h=h,silent=TRUE)$forecast;
}
}
else if(is.null(newxreg) && (nrow(object$xreg)>=(obsInSample+h))){
newxreg <- object$xreg[obsInSample+c(1:h),,drop=FALSE];
}
else{
# If this is not a matrix / data.frame, then convert to one
if(!is.data.frame(newxreg) && !is.matrix(newxreg)){
newxreg <- as.data.frame(newxreg);
colnames(newxreg) <- "xreg";
}
if(nrow(newxreg)<h){
warning(paste0("The newxreg has ",nrow(newxreg)," observations, while ",h," are needed. ",
"Using the last available values as future ones."),
call.=FALSE);
newnRows <- h-nrow(newxreg);
xreg <- rbind(newxreg,matrix(rep(tail(newxreg,1),each=newnRows),newnRows,ncol(newxreg)));
}
else if(nrow(newxreg)>h){
warning(paste0("The newxreg has ",nrow(newxreg)," observations, while only ",h," are needed. ",
"Using the last ",h," of them."),
call.=FALSE);
xreg <- tail(newxreg,h);
}
else{
xreg <- newxreg;
}
xregNames <- colnames(object$xreg);
if(is.data.frame(xreg)){
testFormula <- formula(object);
testFormula[[2]] <- NULL;
# Expand the variables and use only those that are in the model
newxreg <- model.frame(testFormula, xreg);
newxreg <- model.matrix(newxreg,data=newxreg)[,xregNames];
}
else{
newxreg <- xreg[,xregNames];
}
rm(xreg);
}
matWt[,componentsNumberETS+componentsNumberARIMA+c(1:xregNumber)] <- newxreg;
# If this is not "adapt", then fill in the matrix with zeroes
if(object$xregDo!="adapt"){
vecG <- matrix(c(vecG,rep(0,xregNumber)),ncol=1);
}
}
else{
xregNumber <- 0;
}
matF <- object$transition;
# Produce point forecasts for non-multiplicative trend / seasonality
if(Ttype!="M" && Stype!="M"){
adamForecast <- adamForecasterWrap(matVt, matWt, matF,
lagsModelAll, Etype, Ttype, Stype,
componentsNumberETS, componentsNumberETSSeasonal,
componentsNumberARIMA, xregNumber,
h);
}
else{
# If we do simulations, leave it for later
if(interval=="simulated"){
adamForecast <- rep(0, h);
}
# If we don't do simulations to get mean
else{
adamForecast <- forecast(object, h=h, newxreg=newxreg, occurrence=occurrence,
interval="simulated",
level=level, side="both", cumulative=cumulative, nsim=nsim, ...)$mean;
}
}
#### Make safety checks
# If there are NaN values
if(any(is.nan(adamForecast))){
adamForecast[is.nan(adamForecast)] <- 0;
}
# If the occurrence values are provided for the holdout
if(!is.null(occurrence) && is.numeric(occurrence)){
pForecast <- occurrence;
}
else{
# If this is a mixture model, produce forecasts for the occurrence
if(is.occurrence(object$occurrence)){
occurrenceModel <- TRUE;
if(is.alm(object$occurrence)){
pForecast <- forecast(object$occurrence,h=h,newdata=newxreg)$mean;
}
else{
pForecast <- forecast(object$occurrence,h=h,newxreg=newxreg)$mean;
}
}
else{
occurrenceModel <- FALSE;
# If this was provided occurrence, then use provided values
if(!is.null(object$occurrence) && !is.null(object$occurrence$occurrence) &&
(object$occurrence$occurrence=="provided")){
pForecast <- object$occurrence$forecast;
}
else{
pForecast <- rep(1, h);
}
}
}
# Make sure that the values are of the correct length
if(h<length(pForecast)){
pForecast <- pForecast[1:h];
}
else if(h>length(pForecast)){
pForecast <- c(pForecast, rep(tail(pForecast,1), h-length(pForecast)));
}
# How many levels did user asked to produce
nLevels <- length(level);
# Cumulative forecasts have only one observation
if(cumulative){
# hFinal is the number of elements we will have in the final forecast
hFinal <- 1;
}
else{
hFinal <- h;
}
# Create necessary matrices for the forecasts
if(any(yClasses=="ts")){
yForecast <- ts(vector("numeric", hFinal), start=yForecastStart, frequency=yFrequency);
yUpper <- yLower <- ts(matrix(0,hFinal,nLevels), start=yForecastStart, frequency=yFrequency);
}
else{
yForecast <- zoo(vector("numeric", hFinal), order.by=yForecastIndex);
yUpper <- yLower <- zoo(matrix(0,hFinal,nLevels), order.by=yForecastIndex);
}
# Fill in the point forecasts
if(cumulative){
yForecast[] <- sum(adamForecast * pForecast);
}
else{
yForecast[] <- adamForecast * pForecast;
}
if(interval!="none"){
# Fix just in case a silly user used 95 etc instead of 0.95
if(any(level>1)){
level[] <- level / 100;
}
levelLow <- levelUp <- matrix(0,hFinal,nLevels);
levelNew <- matrix(level,nrow=hFinal,ncol=nLevels,byrow=TRUE);
# If this is an occurrence model, then take probability into account in the level.
# This correction is only needed for approximate, because the others contain zeroes
if(occurrenceModel && interval=="approximate"){
levelNew[] <- (levelNew-(1-as.vector(pForecast)))/as.vector(pForecast);
levelNew[levelNew<0] <- 0;
}
if(side=="both"){
levelLow[] <- (1-levelNew)/2;
levelUp[] <- (1+levelNew)/2;
}
else if(side=="upper"){
levelLow[] <- 0;
levelUp[] <- levelNew;
}
else{
levelLow[] <- 1-levelNew;
levelUp[] <- 1;
}
levelLow[levelLow<0] <- 0;
levelUp[levelUp<0] <- 0;
}
#### Simulated interval ####
if(interval=="simulated"){
arrVt <- array(NA, c(componentsNumberETS+componentsNumberARIMA+xregNumber, h+lagsModelMax, nsim));
arrVt[,1:lagsModelMax,] <- rep(matVt,nsim);
sigmaValue <- sigma(object);
matErrors <- matrix(switch(object$distribution,
"dnorm"=rnorm(h*nsim, 0, sigmaValue),
"dlaplace"=rlaplace(h*nsim, 0, sigmaValue/2),
"ds"=rs(h*nsim, 0, (sigmaValue^2/120)^0.25),
"dgnorm"=rgnorm(h*nsim, 0, sigmaValue*sqrt(gamma(1/object$lambda)/gamma(3/object$lambda))),
"dlogis"=rlogis(h*nsim, 0, sigmaValue*sqrt(3)/pi),
"dt"=rt(h*nsim, obsInSample-nparam(object)),
"dalaplace"=ralaplace(h*nsim, 0,
sqrt(sigmaValue^2*object$lambda^2*(1-object$lambda)^2/
(object$lambda^2+(1-object$lambda)^2)),
object$lambda),
"dlnorm"=rlnorm(h*nsim, 0, sigmaValue)-1,
"dinvgauss"=rinvgauss(h*nsim, 1, dispersion=sigmaValue^2)-1,
"dllaplace"=exp(rlaplace(h*nsim, 0, sigmaValue/2))-1,
"dls"=exp(rs(h*nsim, 0, (sigmaValue^2/120)^0.25))-1,
"dlgnorm"=exp(rgnorm(h*nsim, 0, sigmaValue*sqrt(gamma(1/object$lambda)/gamma(3/object$lambda))))-1
),
h,nsim);
# This stuff is needed in order to produce adequate values for weird models
EtypeModified <- Etype;
if(Etype=="A" && any(object$distribution==c("dlnorm","dinvgauss","dls","dllaplace"))){
EtypeModified[] <- "M";
}
# States, Errors, Ot, Transition, Measurement, Persistence
ySimulated <- adamSimulatorwrap(arrVt, matErrors,
# matrix(rep(1,h*nsim), h, nsim),
matrix(rbinom(h*nsim, 1, pForecast), h, nsim),
array(matF,c(dim(matF),nsim)), matWt,
matrix(vecG, componentsNumberETS+componentsNumberARIMA+xregNumber, nsim),
EtypeModified, Ttype, Stype, lagsModelAll,
componentsNumberETSSeasonal, componentsNumberETS,
componentsNumberARIMA, xregNumber)$matrixYt;
#### Note that the cumulative doesn't work with oes at the moment!
if(cumulative){
yForecast[] <- mean(colSums(ySimulated,na.rm=T));
yLower[] <- quantile(colSums(ySimulated,na.rm=T),levelLow,type=7);
yUpper[] <- quantile(colSums(ySimulated,na.rm=T),levelUp,type=7);
}
else{
for(i in 1:h){
if(Ttype=="M" || Stype=="M"){
yForecast[i] <- mean(ySimulated[i,],na.rm=T);
}
yLower[i,] <- quantile(ySimulated[i,],levelLow[i,],na.rm=T,type=7);
yUpper[i,] <- quantile(ySimulated[i,],levelUp[i,],na.rm=T,type=7);
}
}
# This step is needed in order to make intervals similar between the different methods
if(Etype=="A"){
yLower[] <- yLower - yForecast;
yUpper[] <- yUpper - yForecast;
}
else{
yLower[] <- yLower / yForecast;
yUpper[] <- yUpper / yForecast;
}
}
else{
#### Approximate and confidence interval ####
# Produce covatiance matrix and use it
if(any(interval==c("approximate","confidence"))){
s2 <- sigma(object)^2;
# IG and Lnorm can use approximations from the multiplications
if(any(object$distribution==c("dinvgauss","dlnorm","dllaplace","dls","dlgnorm")) && Etype=="M"){
vcovMulti <- adamVarAnal(lagsModelAll, h, matWt[1,,drop=FALSE], matF, vecG, s2);
if(any(object$distribution==c("dlnorm","dls","dllaplace","dlgnorm"))){
vcovMulti[] <- log(1+vcovMulti);
}
# The confidence interval relies on the assumption that initial level is known
if(interval=="confidence"){
vcovMulti[] <- vcovMulti - s2;
}
# We don't do correct cumulatives in this case...
if(cumulative){
vcovMulti <- sum(vcovMulti);
}
}
else{
vcovMulti <- covarAnal(lagsModelAll, h, matWt[1,,drop=FALSE], matF, vecG, s2);
# The confidence interval relies on the assumption that initial level is known
if(interval=="confidence"){
vcovMulti[] <- vcovMulti - s2;
}
# Do either the variance of sum, or a diagonal
if(cumulative){
vcovMulti <- sum(vcovMulti);
}
else{
vcovMulti <- diag(vcovMulti);
}
}
}
#### Semiparametric and nonparametric interval ####
# Extract multistep errors and calculate the covariance matrix
else if(any(interval==c("semiparametric","nonparametric"))){
if(h>1){
adamErrors <- as.matrix(rmultistep(object, h=h));
if(any(object$distribution==c("dinvgauss","dlnorm","dls","dllaplace","dlgnorm")) && (Etype=="A")){
yFittedMatrix <- adamErrors;
for(i in 1:h){
yFittedMatrix[,i] <- fitted(object)[1:(obsInSample-h)+i];
}
adamErrors[] <- adamErrors/yFittedMatrix;
}
if(interval=="semiparametric"){
# Do either the variance of sum, or a diagonal
if(cumulative){
vcovMulti <- sum(t(adamErrors) %*% adamErrors / (obsInSample-h));
}
else{
vcovMulti <- diag(t(adamErrors) %*% adamErrors / (obsInSample-h));
}
}
# For nonparametric and cumulative...
else{
if(cumulative){
adamErrors <- matrix(apply(adamErrors, 2, sum),obsInSample-h,1);
}
}
}
else{
vcovMulti <- sigma(object)^2;
adamErrors <- as.vector(resid(object));
}
}
# Calculate interval for approximate and semiparametric
if(any(interval==c("approximate","confidence","semiparametric"))){
if(object$distribution=="dnorm"){
if(Etype=="A"){
yLower[] <- qnorm(levelLow, 0, sqrt(vcovMulti));
yUpper[] <- qnorm(levelUp, 0, sqrt(vcovMulti));
}
else{
yLower[] <- qnorm(levelLow, 1, sqrt(vcovMulti));
yUpper[] <- qnorm(levelUp, 1, sqrt(vcovMulti));
}
}
else if(object$distribution=="dlaplace"){
if(Etype=="A"){
yLower[] <- qlaplace(levelLow, 0, sqrt(vcovMulti/2));
yUpper[] <- qlaplace(levelUp, 0, sqrt(vcovMulti/2));
}
else{
yLower[] <- qlaplace(levelLow, 1, sqrt(vcovMulti/2));
yUpper[] <- qlaplace(levelUp, 1, sqrt(vcovMulti/2));
}
}
else if(object$distribution=="ds"){
if(Etype=="A"){
yLower[] <- qs(levelLow, 0, (vcovMulti/120)^0.25);
yUpper[] <- qs(levelUp, 0, (vcovMulti/120)^0.25);
}
else{
yLower[] <- qs(levelLow, 1, (vcovMulti/120)^0.25);
yUpper[] <- qs(levelUp, 1, (vcovMulti/120)^0.25);
}
}
else if(object$distribution=="dgnorm"){
alpha <- sqrt(vcovMulti*(gamma(1/object$lambda)/gamma(3/object$lambda)));
if(Etype=="A"){
yLower[] <- suppressWarnings(qgnorm(levelLow, 0, alpha, object$lambda));
yUpper[] <- suppressWarnings(qgnorm(levelUp, 0, alpha, object$lambda));
}
else{
yLower[] <- suppressWarnings(qgnorm(levelLow, 1, alpha, object$lambda));
yUpper[] <- suppressWarnings(qgnorm(levelUp, 1, alpha, object$lambda));
}
}
else if(object$distribution=="dlogis"){
if(Etype=="A"){
yLower[] <- qlogis(levelLow, 0, sqrt(vcovMulti*3)/pi);
yUpper[] <- qlogis(levelUp, 0, sqrt(vcovMulti*3)/pi);
}
else{
yLower[] <- qlogis(levelLow, 1, sqrt(vcovMulti*3)/pi);
yUpper[] <- qlogis(levelUp, 1, sqrt(vcovMulti*3)/pi);
}
}
else if(object$distribution=="dt"){
df <- nobs(object) - nparam(object);
if(Etype=="A"){
yLower[] <- sqrt(vcovMulti)*qt(levelLow, df);
yUpper[] <- sqrt(vcovMulti)*qt(levelUp, df);
}
else{
yLower[] <- (1 + sqrt(vcovMulti)*qt(levelLow, df));
yUpper[] <- (1 + sqrt(vcovMulti)*qt(levelUp, df));
}
}
else if(object$distribution=="dalaplace"){
lambda <- object$lambda;
if(Etype=="A"){
yLower[] <- qalaplace(levelLow, 0,
sqrt(vcovMulti*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
yUpper[] <- qalaplace(levelUp, 0,
sqrt(vcovMulti*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
}
else{
yLower[] <- qalaplace(levelLow, 1,
sqrt(vcovMulti*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
yUpper[] <- qalaplace(levelUp, 1,
sqrt(vcovMulti*lambda^2*(1-lambda)^2/(lambda^2+(1-lambda)^2)), lambda);
}
}
else if(object$distribution=="dlnorm"){
yLower[] <- qlnorm(levelLow, 0, sqrt(vcovMulti));
yUpper[] <- qlnorm(levelUp, 0, sqrt(vcovMulti));
if(Etype=="A"){
yLower[] <- (yLower-1)*yForecast;
yUpper[] <-(yUpper-1)*yForecast;
}
}
else if(object$distribution=="dllaplace"){
yLower[] <- exp(qlaplace(levelLow, 0, sqrt(vcovMulti/2)));
yUpper[] <- exp(qlaplace(levelUp, 0, sqrt(vcovMulti/2)));
if(Etype=="A"){
yLower[] <- (yLower-1)*yForecast;
yUpper[] <-(yUpper-1)*yForecast;
}
}
else if(object$distribution=="dls"){
yLower[] <- exp(qs(levelLow, 0, (vcovMulti/120)^0.25));
yUpper[] <- exp(qs(levelUp, 0, (vcovMulti/120)^0.25));
if(Etype=="A"){
yLower[] <- (yLower-1)*yForecast;
yUpper[] <-(yUpper-1)*yForecast;
}
}
else if(object$distribution=="dlgnorm"){
alpha <- sqrt(vcovMulti*(gamma(1/object$lambda)/gamma(3/object$lambda)));
yLower[] <- suppressWarnings(exp(qgnorm(levelLow, 0, alpha, object$lambda)));
yUpper[] <- suppressWarnings(exp(qgnorm(levelUp, 0, alpha, object$lambda)));
}
else if(object$distribution=="dinvgauss"){
yLower[] <- qinvgauss(levelLow, 1, dispersion=vcovMulti);
yUpper[] <- qinvgauss(levelUp, 1, dispersion=vcovMulti);
if(Etype=="A"){
yLower[] <- (yLower-1)*yForecast;
yUpper[] <-(yUpper-1)*yForecast;
}
}
}
# Use Taylor & Bunn approach for the nonparametric ones
#### Nonparametric intervals, regression ####
else if(interval=="nonparametric"){
if(h>1){
# This is needed in order to see if quant regression can be used
if(all(levelLow==levelLow[1,])){
levelLow <- levelLow[1,,drop=FALSE];
}
if(all(levelUp==levelUp[1,])){
levelUp <- levelUp[1,,drop=FALSE];
}
# Do quantile regression for h>1 and scalars for the level (no change across h)
# transpose is needed in order to compare correctly
if(all(t(levelNew)==levelNew[1,])){
# Quantile regression function
intervalQuantile <- function(A, alpha){
ee[] <- adamErrors - (A[1]*xe^A[2]);
return((1-alpha)*sum(abs(ee[ee<0]))+alpha*sum(abs(ee[ee>=0])));
}
ee <- adamErrors;
xe <- matrix(c(1:h),nrow=nrow(ee),ncol=ncol(ee),byrow=TRUE);
for(i in 1:nLevels){
# lower quantiles
A <- nlminb(rep(1,2),intervalQuantile,alpha=levelLow[1,i])$par;
yLower[,i] <- A[1]*c(1:h)^A[2];
# upper quantiles
A[] <- nlminb(rep(1,2),intervalQuantile,alpha=levelUp[1,i])$par;
yUpper[,i] <- A[1]*c(1:h)^A[2];
}
}
# Otherwise just return quantiles of errors
else{
if(cumulative){
yLower[] <- quantile(adamErrors,levelLow,type=7);
yUpper[] <- quantile(adamErrors,levelUp,type=7);
}
else{
for(i in 1:h){
yLower[i] <- quantile(adamErrors[,i],levelLow[i],na.rm=T,type=7);
yUpper[i] <- quantile(adamErrors[,i],levelUp[i],na.rm=T,type=7);
}
}
}
}
else{
yLower[] <- quantile(adamErrors,levelLow,type=7);
yUpper[] <- quantile(adamErrors,levelUp,type=7);
}
if(Etype=="M"){
yLower[] <- 1+yLower;
yUpper[] <- 1+yUpper;
}
else if(Etype=="A" & any(object$distribution==c("dinvgauss","dlnorm","dllaplace","dls","dlgnorm"))){
yLower[] <- yLower*yForecast;
yUpper[] <- yUpper*yForecast;
}
}
else{
yUpper[] <- yLower[] <- NA;
}
}
# Fix of prediction intervals depending on what has happened
if(interval!="none"){
# Make sensible values out of those weird quantiles
if(!cumulative){
if(Etype=="A"){
yLower[levelLow==0] <- -Inf;
}
else{
yLower[levelLow==0] <- 0;
}
if(any(levelUp==1)){
yUpper[levelUp==1] <- Inf;
}
}
else{
if(Etype=="A" && any(levelLow==0)){
yLower[] <- -Inf;
}
else if(Etype=="M" && any(levelLow==0)){
yLower[] <- 0;
}
if(any(levelUp==1)){
yUpper[] <- Inf;
}
}
# Substitute NAs and NaNs with zeroes
if(any(is.nan(yLower)) || any(is.na(yLower))){
yLower[is.nan(yLower)] <- switch(Etype,"A"=0,"M"=1);
yLower[is.na(yLower)] <- switch(Etype,"A"=0,"M"=1);
}
if(any(is.nan(yUpper)) || any(is.na(yUpper))){
yUpper[is.nan(yUpper)] <- switch(Etype,"A"=0,"M"=1);
yUpper[is.na(yUpper)] <- switch(Etype,"A"=0,"M"=1);
}
# Do intervals around the forecasts...
if(Etype=="A"){
yLower[] <- yForecast + yLower;
yUpper[] <- yForecast + yUpper;
}
else{
yLower[] <- yForecast*yLower;
yUpper[] <- yForecast*yUpper;
}
# Check what we have from the occurrence model
if(occurrenceModel){
# If there are NAs, then there's no variability and no intervals.
if(any(is.na(yUpper))){
yUpper[is.na(yUpper)] <- (yForecast/pForecast)[is.na(yUpper)];
}
if(any(is.na(yLower))){
yLower[is.na(yLower)] <- 0;
}
}
colnames(yLower) <- switch(side,
"both"=paste0("Lower bound (",(1-level)/2*100,"%)"),
"lower"=paste0("Lower bound (",(1-level)*100,"%)"),
"upper"=rep("Lower 0%",nLevels));
colnames(yUpper) <- switch(side,
"both"=paste0("Upper bound (",(1+level)/2*100,"%)"),
"lower"=rep("Upper 100%",nLevels),
"upper"=paste0("Upper bound (",level*100,"%)"));
}
return(structure(list(mean=yForecast, lower=yLower, upper=yUpper, model=object,
level=level, interval=interval, side=side, cumulative=cumulative),
class=c("adam.forecast","smooth.forecast","forecast")));
}
#' @export
forecast.adamCombined <- function(object, h=10, newxreg=NULL,
interval=c("none", "simulated", "approximate", "semiparametric", "nonparametric"),
level=0.95, side=c("both","upper","lower"), cumulative=FALSE, nsim=5000, ...){
interval <- match.arg(interval[1],c("none", "simulated", "approximate", "semiparametric",
"nonparametric", "confidence", "parametric","prediction"));
side <- match.arg(side);
yClasses <- class(actuals(object));
obsInSample <- nobs(object);
if(any(yClasses=="ts")){
# ts structure
yForecastStart <- time(actuals(object))[obsInSample]+deltat(actuals(object));
yFrequency <- frequency(actuals(object));
}
else{
# zoo thingy
yIndex <- time(actuals(object));
yForecastIndex <- yIndex[obsInSample]+diff(tail(yIndex,2))*c(1:h);
}
# How many levels did user asked to produce
nLevels <- length(level);
# Cumulative forecasts have only one observation
if(cumulative){
# hFinal is the number of elements we will have in the final forecast
hFinal <- 1;
}
else{
hFinal <- h;
}
# Create necessary matrices for the forecasts
if(any(yClasses=="ts")){
yForecast <- ts(vector("numeric", hFinal), start=yForecastStart, frequency=yFrequency);
yUpper <- yLower <- ts(matrix(0,hFinal,nLevels), start=yForecastStart, frequency=yFrequency);
}
else{
yForecast <- zoo(vector("numeric", hFinal), order.by=yForecastIndex);
yUpper <- yLower <- zoo(matrix(0,hFinal,nLevels), order.by=yForecastIndex);
}
# The list contains 8 elements
adamForecasts <- vector("list",8);
names(adamForecasts)[c(1:3)] <- c("mean","lower","upper");
for(i in 1:length(object$models)){
adamForecasts[] <- forecast.adam(object$models[[i]], h=h, newxreg=newxreg,
interval=interval,
level=level, side=side, cumulative=cumulative, nsim=nsim, ...);
yForecast[] <- yForecast + adamForecasts$mean * object$ICw[i];
yUpper[] <- yUpper + adamForecasts$upper * object$ICw[i];
yLower[] <- yLower + adamForecasts$lower * object$ICw[i];
}
# Fix the names of the columns
if(interval!="none"){
colnames(yLower) <- colnames(adamForecasts$lower);
colnames(yUpper) <- colnames(adamForecasts$upper);
}
# Fix the content of upper / lower bounds
if(side=="upper"){
yLower[] <- -Inf;
}
else if(side=="lower"){
yUpper[] <- Inf;
}
# Get rid of specific models
object$models <- NULL;
return(structure(list(mean=yForecast, lower=yLower, upper=yUpper, model=object,
level=level, interval=interval, side=side, cumulative=cumulative),
class=c("adam.forecast","smooth.forecast","forecast")));
}
#' @export
print.adam.forecast <- function(x, ...){
if(x$interval!="none"){
returnedValue <- switch(x$side,
"both"=cbind(x$mean,x$lower,x$upper),
"lower"=cbind(x$mean,x$lower),
"upper"=cbind(x$mean,x$upper));
colnames(returnedValue) <- switch(x$side,
"both"=c("Point forecast",colnames(x$lower),colnames(x$upper)),
"lower"=c("Point forecast",colnames(x$lower)),
"upper"=c("Point forecast",colnames(x$upper)))
}
else{
returnedValue <- x$mean;
}
print(returnedValue);
}
#' @export
plot.adam.forecast <- function(x, ...){
yClasses <- class(actuals(x));
ellipsis <- list(...);
if(is.null(ellipsis$ylim)){
vectorOfValues <- switch(x$side,
"both"=c(as.vector(actuals(x$model)),as.vector(x$mean),
as.vector(x$lower),as.vector(x$upper)),
"lower"=c(as.vector(actuals(x$model)),as.vector(x$mean),
as.vector(x$lower)),
"upper"=c(as.vector(actuals(x$model)),as.vector(x$mean),
as.vector(x$upper)));
ellipsis$ylim <- range(vectorOfValues[is.finite(vectorOfValues)],na.rm=TRUE);
}
if(is.null(ellipsis$legend)){
ellipsis$legend <- FALSE;
ellipsis$parReset <- FALSE;
}
if(is.null(ellipsis$main)){
distrib <- switch(x$model$distribution,
"dnorm" = "Normal",
"dlogis" = "Logistic",
"dlaplace" = "Laplace",
"ds" = "S",
"dgnorm" = "Generalised Normal",
"dalaplace" = paste0("Asymmetric Laplace with lambda=",round(x$model$lambda,digits)),
"dt" = paste0("Student t with df=",round(x$model$lambda, digits)),
"dlnorm" = "Log Normal",
"dllaplace" = "Log Laplace",
"dls" = "Log S",
"dgnorm" = "Log Generalised Normal",
# "dbcnorm" = paste0("Box-Cox Normal with lambda=",round(x$other$lambda,2)),
"dinvgauss" = "Inverse Gaussian",
"default"
);
ellipsis$main <- paste0("Forecast from ",x$model$model," with ",distrib," distribution");
}
if(!is.null(x$model$holdout)){
if(any(yClasses=="ts")){
ellipsis$actuals <- ts(c(actuals(x$model),x$model$holdout),
start=start(actuals(x$model)),
frequency=frequency(actuals(x$model)));
}
else{
ellipsis$actuals <- zoo(c(as.vector(actuals(x$model)),as.vector(x$model$holdout)),
order.by=c(time(actuals(x$model)),time(x$model$holdout)));
}
}
else{
ellipsis$actuals <- actuals(x$model);
}
ellipsis$forecast <- x$mean;
ellipsis$fitted <- fitted(x);
ellipsis$lower <- x$lower;
ellipsis$upper <- x$upper;
ellipsis$level <- x$level;
do.call(graphmaker, ellipsis);
}
#### Other methods ####
#' @export
multicov.adam <- function(object, type=c("analytical","empirical","simulated"), ...){
type <- match.arg(type);
# Model type
Ttype <- substr(modelType(object),2,2);
h <- length(object$holdout);
lagsModelAll <- modelLags(object);
lagsModelMax <- max(lagsModelAll);
lagsOriginal <- lags(object);
if(Ttype!="N"){
lagsOriginal <- c(1,lagsOriginal);
}
componentsNumberETS <- length(lagsOriginal);
componentsNumberETSSeasonal <- sum(lagsOriginal>1);
componentsNumberARIMA <- sum(substr(colnames(object$states),1,10)=="ARIMAState");
s2 <- sigma(object)^2;
matWt <- tail(object$measurement,h);
vecG <- matrix(object$persistence, ncol=1);
if(!is.null(object$xreg)){
xregNumber <- ncol(object$xreg);
# matWt[,componentsNumberETS+componentsNumberARIMA+c(1:xregNumber)] <- newxreg[1:h,];
}
else{
xregNumber <- 0;
}
matF <- object$transition;
if(type=="analytical"){
covarMat <- covarAnal(lagsModelAll, h, matWt[1,,drop=FALSE], matF, vecG, s2);
}
else if(type=="empirical"){
adamErrors <- rmultistep(object, h=h);
covarMat <- t(adamErrors) %*% adamErrors / (nobs(object) - h);
}
return(covarMat);
}
#' @export
pointLik.adam <- function(object, ...){
distribution <- object$distribution;
yInSample <- actuals(object);
obsInSample <- nobs(object);
if(is.occurrence(object$occurrence)){
otLogical <- yInSample!=0;
yFitted <- fitted(object) / fitted(object$occurrence);
}
else{
otLogical <- rep(TRUE, obsInSample);
yFitted <- fitted(object);
}
scale <- object$scale;
lambda <- object$lambda;
Etype <- errorType(object);
likValues <- vector("numeric",obsInSample);
likValues[otLogical] <- switch(distribution,
"dnorm"=switch(Etype,
"A"=dnorm(x=yInSample[otLogical], mean=yFitted[otLogical],
sd=scale, log=TRUE),
"M"=dnorm(x=yInSample[otLogical], mean=yFitted[otLogical],
sd=scale*yFitted[otLogical], log=TRUE)),
"dlaplace"=switch(Etype,
"A"=dlaplace(q=yInSample[otLogical], mu=yFitted[otLogical],
scale=scale, log=TRUE),
"M"=dlaplace(q=yInSample[otLogical], mu=yFitted[otLogical],
scale=scale*yFitted[otLogical], log=TRUE)),
"ds"=switch(Etype,
"A"=ds(q=yInSample[otLogical],mu=yFitted[otLogical],
scale=scale, log=TRUE),
"M"=ds(q=yInSample[otLogical],mu=yFitted[otLogical],
scale=scale*sqrt(yFitted[otLogical]), log=TRUE)),
"dgnorm"=switch(Etype,
"A"=dgnorm(x=yInSample[otLogical], mu=yFitted[otLogical],
alpha=scale, beta=lambda, log=TRUE),
"M"=suppressWarnings(dgnorm(x=yInSample[otLogical], mu=yFitted[otLogical],
alpha=scale*yFitted[otLogical], beta=lambda,
log=TRUE))),
"dlogis"=switch(Etype,
"A"=dlogis(x=yInSample[otLogical], location=yFitted[otLogical],
scale=scale, log=TRUE),
"M"=dlogis(x=yInSample[otLogical], location=yFitted[otLogical],
scale=scale*yFitted[otLogical], log=TRUE)),
"dt"=switch(Etype,
"A"=dt(adamFitted$errors[otLogical], df=abs(lambda), log=TRUE),
"M"=dt(adamFitted$errors[otLogical]*yFitted[otLogical],
df=abs(lambda), log=TRUE)),
"dalaplace"=switch(Etype,
"A"=dalaplace(q=yInSample[otLogical], mu=yFitted[otLogical],
scale=scale, alpha=lambda, log=TRUE),
"M"=dalaplace(q=yInSample[otLogical], mu=yFitted[otLogical],
scale=scale*yFitted[otLogical], alpha=lambda, log=TRUE)),
"dlnorm"=dlnorm(x=yInSample[otLogical], meanlog=log(yFitted[otLogical]),
sdlog=scale, log=TRUE),
"dllaplace"=dlaplace(q=log(yInSample[otLogical]), mu=log(yFitted[otLogical]),
scale=scale, log=TRUE),
"dls"=ds(q=log(yInSample[otLogical]), mu=log(yFitted[otLogical]),
scale=scale, log=TRUE),
"dgnorm"=dgnorm(x=log(yInSample[otLogical]), mu=log(yFitted[otLogical]),
alpha=scale, beta=lambda, log=TRUE),
"dinvgauss"=dinvgauss(x=yInSample[otLogical], mean=yFitted[otLogical],
dispersion=scale/yFitted[otLogical], log=TRUE));
if(any(distribution==c("dllaplace","dls","dlgnorm"))){
likValues[otLogical] <- likValues[otLogical] - log(yInSample[otLogical]);
}
# If this is a mixture model, take the respective probabilities into account (differential entropy)
if(is.occurrence(object$occurrence)){
likValues[!otLogical] <- -switch(distribution,
"dnorm" =,
"dlnorm" = (log(sqrt(2*pi)*scale)+0.5),
"dlogis" = 2,
"dlaplace" =,
"dllaplace" =,
"dalaplace" = (1 + log(2*scale)),
"dt" = ((scale+1)/2 * (digamma((scale+1)/2)-digamma(scale/2)) +
log(sqrt(scale) * beta(scale/2,0.5))),
"ds" =,
"dls" = (2 + 2*log(2*scale)),
"dgnorm" =,
"dlgnorm" = 1/lambda-log(lambda/(2*scale*gamma(1/lambda))),
"dinvgauss" = (0.5*(log(pi/2)+1+log(scale))));
likValues[] <- likValues + pointLik(object$occurrence);
}
likValues <- ts(likValues, start=start(yFitted), frequency=frequency(yFitted));
return(likValues);
}
##### Other methods to implement #####
# accuracy.adam <- function(object, holdout, ...){}
# simulate.adam <- function(object, nsim=1, seed=NULL, obs=NULL, ...){}
#' @export
modelType.adam <- function(object, ...){
etsModel <- any(unlist(gregexpr("ETS",object$model))!=-1);
if(etsModel){
modelType <- substring(object$model,
unlist(gregexpr("\\(",object$model))+1,
unlist(gregexpr("\\)",object$model))-1)[1];
}
else{
modelType <- "NNN";
}
return(modelType)
}
#' @export
errorType.adam <- function(object, ...){
model <- modelType(object);
if(model=="NNN"){
return(switch(object$distribution,
"dnorm"=,"dlaplace"=,"ds"=,"dgnorm"=,"dlogis"=,"dt"=,"dalaplace"="A",
"dlnorm"=,"dllaplace"=,"dls"=,"dlgnorm"=,"dinvgauss"="M"));
}
else{
return(substr(model,1,1));
}
}
# This is an internal function, no need to export it
# modelLags <- function(object, ...) UseMethod("modelLags")
modelLags.adam <- function(object, ...){
return(object$lagsAll);
}
#' @export
orders.adam <- function(object, ...){
return(object$orders);
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.