R/simoes.R

In config-i1/smooth: Forecasting Using State Space Models

#' Simulate Occurrence Part of ETS model
#'
#' Function generates data using ETS with Single Source of Error as a data
#' generating process for the demand occurrence. As the main output it produces
#' probabilities of occurrence.
#'
#' For the information about the function, see the vignette:
#' \code{vignette("simulate","smooth")}
#'
#' @template ssAuthor
#' @template ssKeywords
#'
#' @template ssGeneralRef
#'
#' @param model Type of ETS model according to [Hyndman et. al., 2008]
#' taxonomy. Can consist of 3 or 4 chars: \code{ANN}, \code{AAN}, \code{AAdN},
#' \code{AAA}, \code{AAdA}, \code{MAdM} etc. The conventional oETS model assumes
#' that the error term is positive, so "MZZ" models are recommended for this.
#' If you use additive error models, then the function will exponentiate the
#' obtained values before transforming them and getting the probability. This
#' is the type of model A.
#' @param obs Number of observations in each generated time series.
#' @param nsim Number of series to generate (number of simulations to do).
#' @param frequency Frequency of generated data. In cases of seasonal models
#' must be greater than 1.
#' @param occurrence Type of occurrence model. See \code{vignette("oes","smooth")}
#' for details.
#' @param randomizer Type of random number generator function used for error
#' term. It is advised to use \code{rlnorm()} or \code{rinvgauss()} in case of
#' multiplicative error models. If a randomiser is used, it is advised to
#' specify the parameters in the ellipsis.
#' @param bounds Type of bounds to use for persistence vector if values are
#' generated. \code{"usual"} - bounds from p.156 by Hyndman et. al., 2008.
#' \code{"restricted"} - similar to \code{"usual"} but with upper bound equal
#' to 0.3. \code{"admissible"} - bounds from tables 10.1 and 10.2 of Hyndman
#' et. al., 2008. Using first letter of the type of bounds also works. These
#' bounds are also used for multiplicative models, but the models are much
#' more restrictive, so weird results might be obtained. Be careful!
#' @param persistence Persistence vector, which includes all the smoothing
#' parameters. Must correspond to the chosen model. The maximum length is 3:
#' level, trend and seasonal smoothing parameters. If \code{NULL}, values are
#' generated.
#' @param phi Value of damping parameter. If trend is not chosen in the model,
#' the parameter is ignored.
#' @param initial Vector of initial states of level and trend. The maximum
#' length is 2. If \code{NULL}, values are generated.
#' @param initialSeason Vector of initial states for seasonal coefficients.
#' Should have length equal to \code{frequency} parameter. If \code{NULL},
#' values are generated.
#' @param modelB Type of model B. This is used in \code{occurrence="general"}
#' and \code{occurrence="inverse-odds-ratio"}.
#' @param persistenceB The persistence vector for the model B.
#' @param phiB Value of damping parameter for the model B.
#' @param initialB Vector of initial states of level and trend for the model B.
#' @param initialSeasonB Vector of initial states for seasonal coefficients for
#' the model B.
#' @param ...  Additional parameters passed to the chosen randomizer. All the
#' parameters should be passed in the order they are used in chosen randomizer.
#' Both model A and model B share the same parameters for the randomizer.
#'
#' @return List of the following values is returned:
#' \itemize{
#' \item \code{model} - Name of ETS model.
#' \item \code{modelA} - Model A, generated using \code{sim.es()} function;
#' \item \code{modelB} - Model B, generated using \code{sim.es()} function;
#' \item \code{probability} - The value of probability generated by the model;
#' \item \code{occurrence} - Type of occurrence used in the model;
#' \item \code{logLik} - Log-likelihood of the constructed model.
#' }
#'
#' @seealso \code{\link[smooth]{oes}, \link[smooth]{sim.es}, \link[stats]{Distributions}}
#'
#' @examples
#'
#' # This example uses rinvgauss function from statmod package.
#' \donttest{oETSMNNIG <- sim.oes(model="MNN",frequency=12,obs=60,
#'                                randomizer="rinvgauss",mean=1,dispersion=0.5)}
#'
#' # A simpler example with log normal distribution
#' oETSMNNlogN <- sim.oes(model="MNN",frequency=12,obs=60,initial=1,
#'                        randomizer="rlnorm",meanlog=0,sdlog=0.1)
#'
#' @export sim.oes
sim.oes <- function(model="MNN", obs=10, nsim=1, frequency=1,
                    occurrence=c("odds-ratio","inverse-odds-ratio","direct","general"),
                    bounds=c("usual","admissible","restricted"),
                    randomizer=c("rlnorm","rinvgauss","rgamma","rnorm"),
                    persistence=NULL, phi=1,
                    initial=NULL, initialSeason=NULL,
                    modelB=model, persistenceB=persistence, phiB=phi,
                    initialB=initial, initialSeasonB=initialSeason, ...){
    # Function generates data using ETS with Single Source of Error as a data generating process.
    #    Copyright (C) 2020 - Inf Ivan Svetunkov

    occurrence <- match.arg(occurrence);

    modelAName <- model;
    modelBName <- modelB;
    modelA <- switch(occurrence,
                     "odds-ratio"=,
                     "direct"=,
                     "general"=sim.es(model=model, obs=obs, nsim=nsim, frequency=frequency,
                                      persistence=persistence, phi=phi,
                                      initial=initial, initialSeason=initialSeason,
                                      bounds=bounds,
                                      randomizer=randomizer, ...),
                     list(data=matrix(0,obs,nsim)));

    modelB <- switch(occurrence,
                     "inverse-odds-ratio"=,
                     "general"=sim.es(model=modelB, obs=obs, nsim=nsim, frequency=frequency,
                                      persistence=persistenceB, phi=phiB,
                                      initial=initialB, initialSeason=initialSeasonB,
                                      bounds=bounds,
                                      randomizer=randomizer, ...),
                     list(data=matrix(0,obs,nsim)));

    # Direct relies on only one model with max / min functions
    if(occurrence=="direct"){
        model <- list(modelA=modelA);
        model$probability <- ts(vector("numeric",obs), frequency=frequency);
        model$probability[] <- modelA$data;
        model$probability[model$probability>1] <- 1;
        model$probability[model$probability<0] <- 0;
    }
    # The others do division of A / (A + B)
    else{
        model <- list(modelA=modelA, modelB=modelB);
        model$probability <- ts(vector("numeric",obs), frequency=frequency);
        # This way we preserve the potential ts structure
        model$probability[] <- modelA$data;
        model$probability[] <- 1 / (1+switch(errorType(modelB), "M"=c(modelB$data), "A"=c(exp(modelB$data))) /
                                        switch(errorType(modelA), "M"=c(modelA$data), "A"=c(exp(modelA$data))));
    }


    # Likelihood value
    if(nsim==1){
        model$logLik <- sum(log(model$probability));
    }
    else{
        model$logLik <- colSums(log(model$probability));
    }
    model$occurrence <- occurrence;
    model$model <- switch(occurrence,
                          "odds-ratio"=paste0("oETS[O](",modelAName,")"),
                          "inverse-odds-ratio"=paste0("oETS[I](",modelBName,")"),
                          "direct"=paste0("oETS[D](",modelAName,")"),
                          "general"=paste0("oETS[G](",modelAName,")(",modelBName,")"));

    return(structure(model,class="oes.sim"));
}