Nothing
#' Functions that extract values from the fitted model
#'
#' These functions allow extracting orders and lags for \code{ssarima()}, \code{gum()} and \code{sma()},
#' type of model from \code{es()} and \code{ces()} and name of model.
#'
#' \code{orders()} and \code{lags()} are useful only for SSARIMA, GUM and SMA. They return \code{NA} for other functions.
#' This can also be applied to \code{arima()}, \code{Arima()} and \code{auto.arima()} functions from stats and forecast packages.
#' \code{modelType()} is useful only for ETS and CES. They return \code{NA} for other functions.
#' This can also be applied to \code{ets()} function from forecast package. \code{errorType}
#' extracts the type of error from the model (either additive or multiplicative). Finally, \code{modelName}
#' returns the name of the fitted model. For example, "ARIMA(0,1,1)". This is purely descriptive and
#' can also be applied to non-smooth classes, so that, for example, you can easily extract the name
#' of the fitted AR model from \code{ar()} function from \code{stats} package.
#'
#' @template ssAuthor
#' @template ssKeywords
#'
#' @aliases orders
#' @param object Model estimated using one of the functions of smooth package.
#' @param ... Currently nothing is accepted via ellipsis.
#' @return Either vector, scalar or list with values is returned.
#' \code{orders()} in case of ssarima returns list of values:
#' \itemize{
#' \item \code{ar} - AR orders.
#' \item \code{i} - I orders.
#' \item \code{ma} - MA orders.
#' }
#' \code{lags()} returns the vector of lags of the model.
#' All the other functions return strings of character.
#' @seealso \link[smooth]{ssarima}, \link[smooth]{msarima}
#' @examples
#'
#' x <- rnorm(100,0,1)
#'
#' # Just as example. orders and lags do not return anything for ces() and es(). But modelType() does.
#' ourModel <- ces(x, h=10)
#' orders(ourModel)
#' lags(ourModel)
#' modelType(ourModel)
#' modelName(ourModel)
#'
#' # And as another example it does the opposite for gum() and ssarima()
#' ourModel <- gum(x, h=10, orders=c(1,1), lags=c(1,4))
#' orders(ourModel)
#' lags(ourModel)
#' modelType(ourModel)
#' modelName(ourModel)
#'
#' # Finally these values can be used for simulate functions or original functions.
#' ourModel <- auto.ssarima(x)
#' ssarima(x, orders=orders(ourModel), lags=lags(ourModel), constant=ourModel$constant)
#' sim.ssarima(orders=orders(ourModel), lags=lags(ourModel), constant=ourModel$constant)
#'
#' ourModel <- es(x)
#' es(x, model=modelType(ourModel))
#' sim.es(model=modelType(ourModel))
#'
#' @rdname orders
#' @export orders
orders <- function(object, ...) UseMethod("orders")
#' @aliases lags
#' @rdname orders
#' @export lags
lags <- function(object, ...) UseMethod("lags")
#' @aliases modelName
#' @rdname orders
#' @export modelName
modelName <- function(object, ...) UseMethod("modelName")
#' @aliases modelType
#' @rdname orders
#' @export modelType
modelType <- function(object, ...) UseMethod("modelType")
##### Likelihood function and stuff #####
#' @importFrom greybox AICc
#' @export
AICc.smooth <- function(object, ...){
llikelihood <- logLik(object);
nParamAll <- nparam(object);
llikelihood <- llikelihood[1:length(llikelihood)];
if(!is.null(object$occurrence)){
obs <- sum(object$fitted!=0);
nParamSizes <- nParamAll - object$nParam[1,3];
IC <- (2*nParamAll - 2*llikelihood +
2*nParamSizes*(nParamSizes + 1) / (obs - nParamSizes - 1));
}
else{
obs <- nobs(object);
IC <- 2*nParamAll - 2*llikelihood + 2 * nParamAll * (nParamAll + 1) / (obs - nParamAll - 1);
}
return(IC);
}
#' @importFrom greybox BICc
#' @export
BICc.smooth <- function(object, ...){
llikelihood <- logLik(object);
nParamAll <- nparam(object);
llikelihood <- llikelihood[1:length(llikelihood)];
if(!is.null(object$occurrence)){
obs <- sum(object$fitted!=0);
nParamSizes <- nParamAll - object$nParam[1,3];
IC <- - 2*llikelihood + (nParamSizes * log(obs) * obs) / (obs - nParamSizes - 1);
}
else{
obs <- nobs(object);
IC <- - 2*llikelihood + (nParamAll * log(obs) * obs) / (obs - nParamAll - 1);
}
return(IC);
}
#' Function returns the multiple steps ahead covariance matrix of forecast errors
#'
#' This function extracts covariance matrix of 1 to h steps ahead forecast errors for
#' \code{ssarima()}, \code{gum()}, \code{sma()}, \code{es()} and \code{ces()} models.
#'
#' The function returns either scalar (if it is a non-smooth model)
#' or the matrix of (h x h) size with variances and covariances of 1 to h steps ahead
#' forecast errors. This is currently done based on empirical values. The analytical ones
#' are more complicated.
#'
#' @template ssAuthor
#' @template ssKeywords
#'
#' @param object Model estimated using one of the functions of smooth package.
#' @param type What method to use in order to produce covariance matrix:
#' \enumerate{
#' \item \code{analytical} - based on the state space structure of the model and the
#' one-step-ahead forecast error. This works for pure additive and pure multiplicative
#' models. The values for the mixed models might be off.
#' \item \code{empirical} - based on the in-sample 1 to h steps ahead forecast errors
#' (works fine on larger samples);
#' \item \code{simulated} - the data is simulated from the estimated model, then the
#' same model is applied to it and then the empirical 1 to h steps ahead forecast
#' errors are produced;
#' }
#' @param h Forecast horizon to use in the calculations.
#' @param nsim Number of iterations to produce in the simulation. Only needed if
#' \code{type="simulated"}
#' @param ... Other parameters passed to simulate function (if \code{type="simulated"}
#' is used). These are \code{obs} and \code{seed}. By default \code{obs=1000}.
#' This approach increases the accuracy of covariance matrix on small samples
#' and intermittent data;
#' @return Scalar in cases of non-smooth functions. (h x h) matrix otherwise.
#'
#' @seealso \link[smooth]{orders}
#' @examples
#'
#' x <- rnorm(100,0,1)
#'
#' # A simple example with a 5x5 covariance matrix
#' ourModel <- ces(x, h=5)
#' multicov(ourModel)
#'
#' @rdname multicov
#' @export multicov
multicov <- function(object, type=c("analytical","empirical","simulated"), h=10, nsim=1000,
...) UseMethod("multicov")
#' @export
multicov.default <- function(object, type=c("analytical","empirical","simulated"), h=10, nsim=1000,
...){
# Function extracts the conditional variances from the model
return(sigma(object)^2);
}
#' @aliases multicov.smooth
#' @rdname multicov
#' @export
multicov.smooth <- function(object, type=c("analytical","empirical","simulated"), h=10, nsim=1000,
...){
# Function extracts the conditional variances from the model
if(is.smoothC(object)){
stop("Sorry, but covariance matrix is not available for the combinations.",
call.=FALSE)
}
type <- substr(type[1],1,1);
if(is.null(object$persistence) & any(type==c("a","s"))){
warning(paste0("The provided model does not contain the components necessary for the ",
"derivation of the covariance matrix.\n",
"Did you combine forecasts? Switching to 'empirical'"),
call.=FALSE);
type <- "e";
}
if(!is.null(object$occurrence) & type=="e"){
warning(paste0("Empirical covariance matrix can be very inaccurate in cases of ",
"intemittent models.\nWe recommend using type='s' or type='a' instead."),
call.=FALSE);
}
# Empirical covariance matrix
if(type=="e"){
if(errorType(object)=="A"){
errors <- rmultistep(object, h=h);
}
else{
errors <- log(1 + rmultistep(object, h=h));
}
if(!is.null(object$occurrence)){
obs <- t((errors!=0)*1) %*% (errors!=0)*1;
obs[obs==0] <- 1;
df <- obs - nparam(object);
df[df<=0] <- obs[df<=0];
}
else{
obs <- matrix(nobs(object),ncol(errors),ncol(errors));
df <- obs - nparam(object);
df[df<=0] <- obs[df<=0];
}
covarMat <- t(errors) %*% errors / df;
}
# Simulated covariance matrix
else if(type=="s"){
smoothType <- smoothType(object);
ellipsis <- list(...);
if(any(names(ellipsis)=="obs")){
obs <- ellipsis$obs;
}
else{
obs <- length(actuals(object));
}
if(any(names(ellipsis)=="seed")){
seed <- ellipsis$seed;
}
else{
seed <- NULL;
}
h <- length(object$forecast);
# if(smoothType=="ETS"){
# smoothFunction <- es;
# }
# # GUM models
# else if(smoothType=="GUM"){
# smoothFunction <- gum;
# }
# # SSARIMA models
# else if(smoothType=="ARIMA"){
# smoothFunction <- ssarima;
# }
# # CES models
# else if(smoothType=="CES"){
# smoothFunction <- ces;
# }
# # SMA models
# else if(smoothType=="SMA"){
# smoothFunction <- sma;
# }
covarArray <- array(NA,c(h,h,nsim));
newData <- simulate(object, nsim=nsim, obs=obs, seed=seed);
for(i in 1:nsim){
# Apply the model to the simulated data
# smoothModel <- smoothFunction(newData$data[,i], model=object, h=h);
# Remove first h-1 and last values.
# errors <- smoothModel$errors;
errors <- rmultistep(object, h=h)
# Transform errors if needed
if(errorType(object)=="M"){
# Remove zeroes if they are present
errors <- errors[!apply(errors==0,1,any),];
errors <- log(1 + errors);
}
# Calculate covariance matrix
if(!is.null(object$occurrence)){
obsInSample <- t((errors!=0)*1) %*% (errors!=0)*1;
obsInSample[obsInSample==0] <- 1;
}
else{
obsInSample <- matrix(nrow(errors),ncol(errors),ncol(errors));
}
covarArray[,,i] <- t(errors) %*% errors / obsInSample;
}
covarMat <- apply(covarArray,c(1,2),mean);
}
# Analytical covariance matrix
else if(type=="a"){
if(!is.null(object$occurrence)){
ot <- (residuals(object)!=0)*1;
}
else{
ot <- rep(1,nobs(object));
}
# h <- length(object$forecast);
lagsModel <- modelLags(object);
s2 <- sigma(object)^2;
persistence <- matrix(object$persistence,length(object$persistence),1);
transition <- object$transition;
measurement <- object$measurement;
covarMat <- covarAnal(lagsModel, h, measurement, transition, persistence, s2);
}
return(covarMat);
# correlation matrix: multicov(test) / sqrt(diag(multicov(test)) %*% t(diag(multicov(test))))
}
#' @importFrom stats logLik
#' @export
logLik.smooth <- function(object, ...){
if(is.null(object$logLik)){
warning("The likelihood of this model is unavailable. Hint: did you use combinations?");
return(NULL);
}
else{
return(structure(object$logLik,nobs=nobs(object),df=nparam(object),class="logLik"));
}
}
#' @export
logLik.smooth.sim <- function(object, ...){
obs <- nobs(object);
return(structure(object$logLik,nobs=obs,df=0,class="logLik"));
}
#' @importFrom stats nobs
#' @export
nobs.smooth <- function(object, ...){
return(length(object$fitted));
}
#' @method nobs smooth.sim
#' @export
nobs.smooth.sim <- function(object, ...){
if(is.null(dim(object$data))){
return(length(object$data));
}
else{
return(nrow(object$data));
}
}
#' @importFrom greybox nparam
#' @export
nparam.smooth <- function(object, ...){
if(is.null(object$nParam)){
warning("Number of parameters of the model is unavailable. Hint: did you use combinations?",
call.=FALSE);
return(NULL);
}
else{
return(object$nParam[1,ncol(object$nParam)]);
}
}
#' Prediction Likelihood Score
#'
#' Function estimates Prediction Likelihood Score for the provided model
#'
#' Prediction likelihood score (PLS) is based on either normal or log-normal
#' distribution of errors. This is extracted from the provided model. The likelihood
#' based on the distribution of 1 to h steps ahead forecast errors is used in the process.
#'
#' @template ssAuthor
#' @template ssKeywords
#'
#' @param object The model estimated using smooth functions. This thing also accepts
#' other models (e.g. estimated using functions from forecast package), but may not always
#' work properly with them.
#' @param holdout The values for the holdout part of the sample. If the model was fitted
#' on the data with the \code{holdout=TRUE}, then the parameter is not needed.
#' @param ... Parameters passed to multicov function. The function is called in order to get
#' the covariance matrix of 1 to h steps ahead forecast errors.
#'
#' @return A value of the log-likelihood.
#' @references \itemize{
#' %\item Eltoft, T., Taesu, K., Te-Won, L. (2006). On the multivariate Laplace
#' distribution. IEEE Signal Processing Letters. 13 (5): 300-303.
#' \doi{10.1109/LSP.2006.870353} - this is not yet used in the function.
#' \item Snyder, R. D., Ord, J. K., Beaumont, A., 2012. Forecasting the intermittent
#' demand for slow-moving inventories: A modelling approach. International
#' Journal of Forecasting 28 (2), 485-496.
#' \item Kolassa, S., 2016. Evaluating predictive count data distributions in retail
#' sales forecasting. International Journal of Forecasting 32 (3), 788-803..
#' }
#' @examples
#'
#' # Generate data, apply es() with the holdout parameter and calculate PLS
#' x <- rnorm(100,0,1)
#' ourModel <- es(x, h=10, holdout=TRUE)
#' pls(ourModel, type="a")
#' pls(ourModel, type="e")
#' pls(ourModel, type="s", obs=100, nsim=100)
#'
#' @rdname pls
#' @export pls
pls <- function(object, holdout=NULL, ...) UseMethod("pls")
# Function calculates PLS based on the provided model
#' @importFrom stats dnorm
#' @export
pls.default <- function(object, holdout=NULL, ...){
if(is.null(holdout)){
stop("We need the values from the holdout in order to proceed.",
call.=FALSE);
}
h <- length(holdout);
yForecast <- forecast(object, h=h)$mean;
return(sum(dnorm(holdout,yForecast,sigma(object),log=TRUE)));
}
#' @rdname pls
#' @aliases pls.smooth
#' @export
pls.smooth <- function(object, holdout=NULL, ...){
if(is.smoothC(object)){
stop("Sorry, but PLS is not available for the combinations.",
call.=FALSE)
}
# If holdout is provided, check it and use it. Otherwise try extracting from the model
yForecast <- object$forecast;
covarMat <- multicov(object, ...);
if(!is.null(holdout)){
if(length(yForecast)!=length(holdout)){
if(is.null(object$holdout)){
stop("The forecast of the model does not correspond to the provided holdout.",
call.=FALSE);
}
else{
holdout <- object$holdout;
}
}
}
else{
if(all(is.na(object$holdout))){
stop("No values for the holdout are available. Cannot proceed.",
call.=FALSE);
}
holdout <- object$holdout;
}
h <- length(holdout);
Etype <- errorType(object);
loss <- object$loss;
if(any(loss==c("MAE","MAEh","TMAE","GTMAE","MACE"))){
loss <- "MAE";
}
else if(any(loss==c("HAM","HAMh","THAM","GTHAM","CHAM"))){
loss <- "HAM";
}
else{
loss <- "MSE";
}
densityFunction <- function(loss, ...){
if(loss=="MAE"){
# This is a simplification. The real multivariate Laplace is bizarre!
scale <- sqrt(diag(covarMat)/2);
plsValue <- sum(dlaplace(errors, 0, scale, log=TRUE));
}
else if(loss=="HAM"){
# This is a simplification. We don't have multivariate HAM yet.
scale <- (diag(covarMat)/120)^0.25;
plsValue <- sum(ds(errors, 0, scale, log=TRUE));
}
else{
if(is.infinite(det(covarMat))){
plsValue <- -as.vector((log(2*pi)+(abs(determinant(covarMat)$modulus)))/2 +
(t(errors) %*% solve(covarMat) %*% errors) / 2);
}
# This is the case with overfitting the data
else if(det(covarMat)==0){
# If there is any non-zero error, then it means that the model is completely wrong (because it predicts that sigma=0)
if(any(errors!=0)){
plsValue <- -Inf;
}
else{
plsValue <- 0;
}
}
else{
# Here and later in the code the abs() is needed for weird cases of wrong covarMat
plsValue <- -as.vector(log(2*pi*abs(det(covarMat)))/2 +
(t(errors) %*% solve(covarMat) %*% errors) / 2);
}
}
return(plsValue);
}
# Additive models
if(Etype=="A"){
# Non-intermittent data
if(is.null(object$occurrence)){
errors <- holdout - yForecast;
plsValue <- densityFunction(loss, errors, covarMat);
}
# Intermittent data
else{
ot <- holdout!=0;
pForecast <- object$occurrence$forecast;
errors <- holdout - yForecast / pForecast;
if(all(ot)){
plsValue <- densityFunction(loss, errors, covarMat) + sum(log(pForecast));
}
else if(all(!ot)){
plsValue <- sum(log(1-pForecast));
}
else{
errors[!ot] <- 0;
plsValue <- densityFunction(loss, errors, covarMat);
plsValue <- plsValue + sum(log(pForecast[ot])) + sum(log(1-pForecast[!ot]));
}
}
}
# Multiplicative models
else{
# Non-intermittent data
if(is.null(object$occurrence)){
errors <- log(holdout) - log(yForecast);
plsValue <- densityFunction(loss, errors, covarMat) - sum(log(holdout));
}
# Intermittent data
else{
ot <- holdout!=0;
pForecast <- object$occurrence$forecast;
errors <- log(holdout) - log(yForecast / pForecast);
if(all(ot)){
plsValue <- (densityFunction(loss, errors, covarMat) - sum(log(holdout)) +
sum(log(pForecast)));
}
else if(all(!ot)){
plsValue <- sum(log(1-pForecast));
}
else{
errors[!ot] <- 0;
plsValue <- densityFunction(loss, errors, covarMat) - sum(log(holdout[ot]));
plsValue <- plsValue + sum(log(pForecast[ot])) + sum(log(1-pForecast[!ot]));
}
}
}
return(plsValue);
}
#' @importFrom stats sigma
#' @export
sigma.smooth <- function(object, ...){
if(!is.null(object$s2)){
return(sqrt(object$s2));
}
else{
return(NULL);
}
}
#' @export
sigma.smooth.sim <- function(object, ...){
return(sqrt(mean(residuals(object)^2)));
}
#### pointLik for smooth ####
#' @importFrom greybox pointLik
#' @export
pointLik.smooth <- function(object, log=TRUE, ...){
obs <- nobs(object);
errors <- residuals(object);
likValues <- vector("numeric",obs);
loss <- object$loss;
if(errorType(object)=="M"){
likValues <- likValues - log(actuals(object));
}
if(any(loss==c("MAE","MAEh","TMAE","GTMAE","MACE"))){
likValues <- likValues + dlaplace(errors, 0, mean(abs(errors)), log=log);
}
else if(any(loss==c("HAM","HAMh","THAM","GTHAM","CHAM"))){
likValues <- likValues + ds(errors, 0, mean(sqrt(abs(errors))/2), log=log);
}
else{
likValues <- likValues + dnorm(errors, 0, sqrt(mean(abs(errors)^2)), log=log);
}
likValues <- ts(as.vector(likValues), start=start(errors), frequency=frequency(errors));
return(likValues);
}
#' @export
pointLik.oes <- function(object, log=TRUE, ...){
ot <- actuals(object);
pFitted <- fitted(object);
likValues <- vector("numeric",nobs(object));
likValues[ot==1] <- log(pFitted[ot==1]);
likValues[ot==0] <- log(1-pFitted[ot==0]);
likValues <- ts(likValues, start=start(ot), frequency=frequency(ot));
if(!log){
likValues[] <- exp(likValues);
}
return(likValues);
}
#### Extraction of parameters of models ####
#' @export
coef.smooth <- function(object, ...)
{
smoothType <- smoothType(object);
if(smoothType=="CES"){
parameters <- c(object$a,object$b);
}
else if(smoothType=="ETS"){
if(any(unlist(gregexpr("C",modelType(object)))==-1)){
# If this was normal ETS, return values
parameters <- c(object$persistence,object$initial,object$initialSeason,object$initialX);
}
else{
# If we did combinations, we cannot return anything
message("Combination of models was done, so there are no coefficients to return");
parameters <- NULL;
}
}
else if(smoothType=="GUM"){
parameters <- c(object$measurement,object$transition,object$persistence,object$initial);
names(parameters) <- c(paste0("Measurement ",c(1:length(object$measurement))),
paste0("Transition ",c(1:length(object$transition))),
paste0("Persistence ",c(1:length(object$persistence))),
paste0("Initial ",c(1:length(object$initial))));
}
else if(smoothType=="ARIMA"){
if(any(unlist(gregexpr("combine",object$model))==-1)){
# If this was normal ARIMA, return values
namesConstant <- NamesMA <- NamesAR <- parameters <- NULL;
if(any(object$AR!=0)){
parameters <- c(parameters,object$AR);
NamesAR <- paste(rownames(object$AR),rep(colnames(object$AR),each=nrow(object$AR)),sep=", ");
}
if(any(object$MA!=0)){
parameters <- c(parameters,object$MA);
NamesMA <- paste(rownames(object$MA),rep(colnames(object$MA),each=nrow(object$MA)),sep=", ")
}
if(object$constant!=0){
parameters <- c(parameters,object$constant);
namesConstant <- "Constant";
}
names(parameters) <- c(NamesAR,NamesMA,namesConstant);
parameters <- parameters[parameters!=0];
}
else{
# If we did combinations, we cannot return anything
message("Combination of models was done, so there are no coefficients to return");
parameters <- NULL;
}
}
else if(smoothType=="SMA"){
parameters <- object$persistence;
}
return(parameters);
}
#### Fitted, forecast and actual values ####
#' @export
fitted.smooth <- function(object, ...){
return(object$fitted);
}
#' @export
fitted.smooth.forecast <- function(object, ...){
return(fitted(object$model));
}
#' Forecasting time series using smooth functions
#'
#' Function produces conditional expectation (point forecasts) and prediction
#' intervals for the estimated model.
#'
#' By default the function will generate conditional expectations from the
#' estimated model and will also produce a variety of prediction intervals
#' based on user preferences.
#'
#' @aliases forecast forecast.smooth
#' @param object Time series model for which forecasts are required.
#' @param h Forecast horizon.
#' @param level Confidence level. Defines width of prediction interval.
#' @param side Defines, whether to provide \code{"both"} sides of prediction
#' interval or only \code{"upper"}, or \code{"lower"}.
#' @param ... Other arguments accepted by either \link[smooth]{es},
#' \link[smooth]{ces}, \link[smooth]{gum} or \link[smooth]{ssarima}.
#' @return Returns object of class "smooth.forecast", which contains:
#'
#' \itemize{
#' \item \code{model} - the estimated model (ES / CES / GUM / SSARIMA).
#' \item \code{method} - the name of the estimated model (ES / CES / GUM / SSARIMA).
#' \item \code{forecast} aka \code{mean} - point forecasts of the model
#' (conditional mean).
#' \item \code{lower} - lower bound of prediction interval.
#' \item \code{upper} - upper bound of prediction interval.
#' \item \code{level} - confidence level.
#' \item \code{interval} - binary variable (whether interval were produced or not).
#' \item \code{scenarios} - in case of \code{forecast.adam()} and
#' \code{interval="simulated"} returns matrix with scenarios (future paths) that were
#' used in simulations.
#' }
#' @template ssAuthor
#' @seealso \code{\link[generics]{forecast}}
#' @references Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008)
#' Forecasting with exponential smoothing: the state space approach,
#' Springer-Verlag.
#' @keywords ts univar
#' @examples
#'
#' ourModel <- ces(rnorm(100,0,1),h=10)
#'
#' forecast(ourModel,h=10)
#' forecast(ourModel,h=10,interval=TRUE)
#' plot(forecast(ourModel,h=10,interval=TRUE))
#'
#' @rdname forecast.smooth
#' @importFrom generics forecast
#' @export
forecast.smooth <- function(object, h=10,
interval=c("parametric","semiparametric","nonparametric","none"),
level=0.95, side=c("both","upper","lower"), ...){
smoothType <- smoothType(object);
interval <- interval[1];
side <- match.arg(side);
# This correction is needed in order to reduce the level and then just use one bound
if(any(side==c("upper","lower"))){
levelNew <- level*2-1;
}
else{
levelNew <- level;
}
# Do calculations
if(smoothType=="ETS"){
newModel <- es_old(actuals(object),model=object,h=h,interval=interval,level=levelNew,silent="all",...);
}
else if(smoothType=="CES"){
newModel <- ces(actuals(object),model=object,h=h,interval=interval,level=levelNew,silent="all",...);
}
else if(smoothType=="GUM"){
newModel <- gum(actuals(object),model=object,type=errorType(object),h=h,interval=interval,level=levelNew,silent="all",...);
}
else if(smoothType=="ARIMA"){
if(any(unlist(gregexpr("combine",object$model))==-1)){
if(is.msarima(object)){
newModel <- msarima(actuals(object),model=object,h=h,interval=interval,level=levelNew,silent="all",...);
}
else{
newModel <- ssarima(actuals(object),model=object,h=h,interval=interval,level=levelNew,silent="all",...);
}
}
else{
stop(paste0("Sorry, but in order to produce forecasts for this ARIMA we need to recombine it.\n",
"You will have to use auto.ssarima() function instead."),call.=FALSE);
}
}
else if(smoothType=="SMA"){
newModel <- sma(actuals(object),model=object,h=h,interval=interval,level=levelNew,silent="all",...);
}
else{
stop("Wrong object provided. This needs to be either 'ETS', or 'CES', or 'GUM', or 'SSARIMA', or 'SMA' model.",call.=FALSE);
}
# Remove the redundant values, if they were produced
if(side=="upper"){
newModel$lower[] <- NA;
newModel$level <- level;
}
else if(side=="lower"){
newModel$upper[] <- NA;
newModel$level <- level;
}
output <- list(model=object,method=object$model,
forecast=newModel$forecast,lower=newModel$lower,upper=newModel$upper,level=newModel$level,
interval=interval,mean=newModel$forecast,side=side);
return(structure(output,class=c("smooth.forecast","forecast")));
}
#' @rdname forecast.smooth
#' @export
forecast.oes <- function(object, h=10,
# interval=c("parametric","semiparametric","nonparametric","none"),
# level=0.95, side=c("both","upper","lower"),
...){
# side <- match.arg(side);
# This correction is needed in order to reduce the level and then just use one bound
# if(any(side==c("upper","lower"))){
# levelNew <- level*2-1;
# }
# else{
# levelNew <- level;
# }
if(is.oesg(object)){
newModel <- oesg(actuals(object),modelA=object$modelA,modelB=object$modelB,
h=h,silent="all",...);
}
else{
newModel <- oes(actuals(object),model=object,h=h,silent="all",...);
}
# Remove the redundant values, if they were produced
# if(side=="upper"){
# newModel$lower[] <- NA;
# newModel$level <- level;
# }
# else if(side=="lower"){
# newModel$upper[] <- NA;
# newModel$level <- level;
# }
output <- list(model=object, method=object$model, mean=newModel$forecast,
forecast=newModel$forecast, interval="none"
# lower=newModel$lower,upper=newModel$upper,level=levelNew,
# interval=interval,side=side
);
return(structure(output,class=c("smooth.forecast","forecast")));
}
#' @importFrom stats window
#' @importFrom greybox actuals
#' @export
actuals.smooth <- function(object, ...){
return(window(object$y,start=start(object$y),end=end(object$fitted)));
}
#' @export
actuals.smooth.forecast <- function(object, ...){
return(window(actuals(object$model),start=start(actuals(object$model)),end=end(fitted(object$model))));
}
#### Function extracts lags of provided model ####
#' @export
lags.default <- function(object, ...){
return(NA);
}
#' @export
lags.ets <- function(object, ...){
modelName <- modelType(object);
lags <- c(1);
if(substr(modelName,nchar(modelName),nchar(modelName))!="N"){
lags <- c(lags,frequency(actuals(object)));
}
return(lags);
}
#' @export
lags.ar <- function(object, ...){
return(1);
}
#' @export
lags.Arima <- function(object, ...){
return(c(1,object$arma[5]));
}
#' @export
lags.smooth <- function(object, ...){
model <- object$model;
smoothType <- smoothType(object);
if(!is.null(model)){
if(smoothType=="GUM"){
lags <- as.numeric(substring(model,unlist(gregexpr("\\[",model))+1,unlist(gregexpr("\\]",model))-1));
}
else if(smoothType=="ARIMA"){
if(any(unlist(gregexpr("combine",object$model))==-1)){
if(any(unlist(gregexpr("\\[",model))!=-1)){
lags <- as.numeric(substring(model,unlist(gregexpr("\\[",model))+1,unlist(gregexpr("\\]",model))-1));
}
else{
lags <- 1;
}
}
else{
warning("ARIMA was combined and we cannot extract lags anymore. Sorry!",call.=FALSE);
}
}
else if(smoothType=="SMA"){
lags <- 1;
}
else if(smoothType=="ETS"){
modelName <- modelType(object);
lags <- c(1);
if(substr(modelName,nchar(modelName),nchar(modelName))!="N"){
lags <- c(lags,frequency(actuals(object)));
}
}
else if(smoothType=="CES"){
modelName <- modelType(object);
dataFreq <- frequency(actuals(object));
if(modelName=="none"){
lags <- c(1);
}
else if(modelName=="simple"){
lags <- c(dataFreq);
}
else if(modelName=="partial"){
lags <- c(1,dataFreq);
}
else if(modelName=="full"){
lags <- c(1,dataFreq);
}
else{
stop("Sorry, but we cannot identify the type of the provided model.",
call.=FALSE);
}
}
else{
lags <- NA;
}
}
else{
lags <- NA;
}
return(lags);
}
#' @export
lags.smooth.sim <- lags.smooth;
#### Function extracts type of error in the model: "A" or "M" ####
#' @importFrom greybox errorType
#' @export
errorType.smooth <- function(object, ...){
smoothType <- smoothType(object);
# ETS models
if(smoothType=="ETS"){
if(any(substr(modelType(object),1,1)==c("A","X"))){
Etype <- "A";
}
else if(any(substr(modelType(object),1,1)==c("M","Y"))){
Etype <- "M";
}
else{
stop("Sorry, but we cannot define error type for this type of model",
call.=FALSE);
}
}
# GUM models
else if(smoothType=="GUM"){
if(substr(modelName(object),1,1)=="M"){
Etype <- "M";
}
else{
Etype <- "A";
}
}
# SSARIMA models
else if(smoothType=="ARIMA"){
Etype <- "A";
}
# CES models
else if(smoothType=="CES"){
Etype <- "A";
}
# SMA models
else if(smoothType=="SMA"){
Etype <- "A";
}
# SMA models
else if(smoothType=="CMA"){
Etype <- "A";
}
else{
stop(paste0("Sorry but we cannot identify error type for the model '",object$model),
call.=FALSE);
}
return(Etype);
}
#' @export
errorType.smooth.sim <- errorType.smooth;
##### Function returns the modelLags from the model - internal function #####
modelLags <- function(object, ...){
modelLags <- NA;
if(is.msarima(object)){
modelLags <- object$modelLags;
}
else if(is.adam(object)){
modelLags <- object$lagsAll;
}
else{
smoothType <- smoothType(object);
if(smoothType=="ETS"){
modelLags <- matrix(rep(lags(object),times=orders(object)),ncol=1);
}
else if(smoothType=="GUM"){
modelLags <- matrix(rep(lags(object),times=orders(object)),ncol=1);
}
else if(smoothType=="ARIMA"){
ordersARIMA <- orders(object);
nComponents <- max(ordersARIMA$ar %*% lags(object) + ordersARIMA$i %*% lags(object),
ordersARIMA$ma %*% lags(object));
modelLags <- matrix(rep(1,times=nComponents),ncol=1);
if(is.numeric(object$constant)){
modelLags <- rbind(modelLags,1);
}
}
else if(smoothType=="CES"){
modelLags <- matrix(rep(lags(object),times=orders(object)),ncol=1);
}
else if(smoothType=="SMA"){
modelLags <- matrix(rep(1,times=orders(object)),ncol=1);
}
}
return(modelLags);
}
#### Function extracts the full name of a model. For example "ETS(AAN)" ####
#' @export
modelName.default <- function(object, ...){
return(NA);
}
#' @export
modelName.ar <- function(object, ...){
return(paste0("AR(",object$order,")"));
}
#' @export
modelName.lm <- function(object, ...){
return("Regression");
}
#' @export
modelName.Arima <- function(object, ...){
ordersArma <- object$arma;
if(all(ordersArma[c(3,4,7)]==0)){
return(paste0("ARIMA(",paste(ordersArma[c(1,6,2)],collapse=","),")"));
}
else{
return(paste0("SARIMA(",paste(ordersArma[c(1,6,2)],collapse=","),")(",
paste(ordersArma[c(3,7,4)],collapse=","),")[",ordersArma[5],"]"));
}
}
#' @export
modelName.ets <- function(object, ...){
return(object$method);
}
#' @export
modelName.forecast <- function(object, ...){
return(object$method);
}
#' @export
modelName.smooth <- function(object, ...){
return(object$model);
}
#### Function extracts type of model. For example "AAN" from ets ####
#' @export
modelType.default <- function(object, ...){
return(NA);
}
#' @export
modelType.smooth <- function(object, ...){
model <- object$model;
smoothType <- smoothType(object);
if(smoothType=="ETS"){
modelType <- substring(model,unlist(gregexpr("\\(",model))+1,unlist(gregexpr("\\)",model))-1);
}
else if(smoothType=="CES"){
modelType <- substring(model,unlist(gregexpr("\\(",model))+1,unlist(gregexpr("\\)",model))-1);
if(modelType=="n"){
modelType <- "none";
}
else if(modelType=="s"){
modelType <- "simple";
}
else if(modelType=="p"){
modelType <- "partial";
}
else{
modelType <- "full";
}
}
else{
modelType <- NA;
}
return(modelType);
}
#' @export
modelType.smooth.sim <- modelType.smooth;
#' @export
modelType.oesg <- function(object, ...){
return(modelType(object$modelA));
}
#' @export
modelType.ets <- function(object, ...){
return(gsub(",","",substring(object$method,5,nchar(object$method)-1)));
}
#### Function extracts orders of provided model ####
#' @export
orders.default <- function(object, ...){
return(NA);
}
#' @export
orders.smooth <- function(object, ...){
smoothType <- smoothType(object);
model <- object$model;
if(!is.null(model)){
if(smoothType=="GUM"){
orders <- as.numeric(substring(model,unlist(gregexpr("\\[",model))-1,unlist(gregexpr("\\[",model))-1));
}
else if(smoothType=="ARIMA"){
if(any(unlist(gregexpr("combine",object$model))==-1)){
arima.orders <- paste0(c("",substring(model,unlist(gregexpr("\\(",model))+1,unlist(gregexpr("\\)",model))-1),"")
,collapse=";");
comas <- unlist(gregexpr("\\,",arima.orders));
semicolons <- unlist(gregexpr("\\;",arima.orders));
ar.orders <- as.numeric(substring(arima.orders,semicolons[-length(semicolons)]+1,comas[2*(1:(length(comas)/2))-1]-1));
i.orders <- as.numeric(substring(arima.orders,comas[2*(1:(length(comas)/2))-1]+1,comas[2*(1:(length(comas)/2))-1]+1));
ma.orders <- as.numeric(substring(arima.orders,comas[2*(1:(length(comas)/2))]+1,semicolons[-1]-1));
orders <- list(ar=ar.orders,i=i.orders,ma=ma.orders);
}
else{
warning("ARIMA was combined and we cannot extract orders anymore. Sorry!",call.=FALSE);
}
}
else if(smoothType=="SMA"){
orders <- as.numeric(substring(model,unlist(gregexpr("\\(",model))+1,unlist(gregexpr("\\)",model))-1));
}
else if(smoothType=="ETS"){
modelName <- modelType(object);
orders <- 1;
if(substr(modelName,2,2)!="N"){
orders <- 2;
}
if(substr(modelName,nchar(modelName),nchar(modelName))!="N"){
orders <- c(orders, 1);
}
}
else if(smoothType=="CES"){
modelName <- modelType(object);
if(modelName=="none"){
orders <- 2;
}
else if(modelName=="simple"){
orders <- 2;
}
else if(modelName=="partial"){
orders <- c(2,1);
}
else if(modelName=="full"){
orders <- c(2,2);
}
else{
stop("Sorry, but we cannot identify the type of the provided model.",
call.=FALSE);
}
}
else{
orders <- NA;
}
}
else{
orders <- NA;
}
return(orders);
}
#' @export
orders.smooth.sim <- orders.smooth;
#' @export
orders.ar <- function(object, ...){
return(list(ar=object$order));
}
#' @export
orders.Arima <- function(object, ...){
model <- object$arma;
ar.orders <- c(model[1],model[3]);
i.orders <- c(model[6],model[7]);
ma.orders <- c(model[2],model[4]);
orders <- list(ar=ar.orders,i=i.orders,ma=ma.orders);
return(orders);
}
#### Plots of smooth objects ####
#' Plots for the fit and states
#'
#' The function produces diagnostics plots for a \code{smooth} model
#'
#' The list of produced plots includes:
#' \enumerate{
#' \item Actuals vs Fitted values. Allows analysing, whether there are any issues in the fit.
#' Does the variability of actuals increase with the increase of fitted values? Is the relation
#' well captured? They grey line on the plot corresponds to the perfect fit of the model.
#' \item Standardised residuals vs Fitted. Plots the points and the confidence bounds
#' (red lines) for the specified confidence \code{level}. Useful for the analysis of outliers;
#' \item Studentised residuals vs Fitted. This is similar to the previous plot, but with the
#' residuals divided by the scales with the leave-one-out approach. Should be more sensitive
#' to outliers;
#' \item Absolute residuals vs Fitted. Useful for the analysis of heteroscedasticity;
#' \item Squared residuals vs Fitted - similar to (3), but with squared values;
#' \item Q-Q plot with the specified distribution. Can be used in order to see if the
#' residuals follow the assumed distribution. The type of distribution depends on the one used
#' in the estimation (see \code{distribution} parameter in \link[greybox]{alm});
#' \item ACF of the residuals. Are the residuals autocorrelated? See \link[stats]{acf} for
#' details;
#' \item Fitted over time. Plots actuals (black line), fitted values (purple line), point forecast
#' (blue line) and prediction interval (grey lines). Can be used in order to make sure that the model
#' did not miss any important events over time;
#' \item Standardised residuals vs Time. Useful if you want to see, if there is autocorrelation or
#' if there is heteroscedasticity in time. This also shows, when the outliers happen;
#' \item Studentised residuals vs Time. Similar to previous, but with studentised residuals;
#' \item PACF of the residuals. No, really, are they autocorrelated? See pacf function from stats
#' package for details;
#' \item Plot of the states of the model. It is not recommended to produce this plot together with
#' the others, because there might be several states, which would cause the creation of a different
#' canvas. In case of "msdecompose", this will produce the decomposition of the series into states
#' on a different canvas;
#' \item Absolute standardised residuals vs Fitted. Similar to the previous, but with absolute
#' values. This is more relevant to the models where scale is calculated as an absolute value of
#' something (e.g. Laplace);
#' \item Squared standardised residuals vs Fitted. This is an additional plot needed to diagnose
#' heteroscedasticity in a model with varying scale. The variance on this plot will be constant if
#' the adequate model for \code{scale} was constructed. This is more appropriate for normal and
#' the related distributions.
#' }
#' Which of the plots to produce, is specified via the \code{which} parameter.
#'
#' @param x Estimated smooth model.
#' @param which Which of the plots to produce. The possible options (see details for explanations):
#' \enumerate{
#' \item Actuals vs Fitted values;
#' \item Standardised residuals vs Fitted;
#' \item Studentised residuals vs Fitted;
#' \item Absolute residuals vs Fitted;
#' \item Squared residuals vs Fitted;
#' \item Q-Q plot with the specified distribution;
#' \item Fitted over time;
#' \item Standardised residuals vs Time;
#' \item Studentised residuals vs Time;
#' \item ACF of the residuals;
#' \item PACF of the residuals;
#' \item Plot of states of the model;
#' \item Absolute standardised residuals vs Fitted;
#' \item Squared standardised residuals vs Fitted;
#' \item ACF of the squared residuals;
#' \item PACF of the squared residuals.
#' }
#' @param level Confidence level. Defines width of confidence interval. Used in plots (2), (3), (7), (8),
#' (9), (10) and (11).
#' @param legend If \code{TRUE}, then the legend is produced on plots (2), (3) and (7).
#' @param ask Logical; if \code{TRUE}, the user is asked to press Enter before each plot.
#' @param lowess Logical; if \code{TRUE}, LOWESS lines are drawn on scatterplots, see \link[stats]{lowess}.
#' @param ... The parameters passed to the plot functions. Recommended to use with separate plots.
#' @return The function produces the number of plots, specified in the parameter \code{which}.
#'
#' @template ssAuthor
#' @seealso \link[greybox]{plot.greybox}
#' @keywords ts univar
#' @examples
#'
#' ourModel <- es(c(rnorm(50,100,10),rnorm(50,120,10)), "ANN", h=10)
#' par(mfcol=c(3,4))
#' plot(ourModel, c(1:11))
#' plot(ourModel, 12)
#'
#' @importFrom stats ppoints qqnorm qqplot qqline acf pacf lowess sd na.pass
#' @importFrom grDevices dev.interactive devAskNewPage
#' @importFrom graphics plot text
#' @importFrom greybox is.occurrence
#' @rdname plot.smooth
#' @export
plot.smooth <- 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"))){
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(!any(names(ellipsis)=="main")){
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";
}
}
# Get the IDs of outliers and statistic
outliers <- outlierdummy(x, level=level, type=type);
outliersID <- outliers$id;
statistic <- outliers$statistic;
# Substitute zeroes with NAs if there was an occurrence
if(is.occurrence(x$occurrence)){
ellipsis$x[actuals(x$occurrence)==0] <- NA;
}
if(!any(names(ellipsis)=="ylim")){
ellipsis$ylim <- range(c(ellipsis$y,statistic), 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(statistic[1],statistic[1],statistic[2],statistic[2]),
col="lightgrey", border=NA, density=10);
abline(h=statistic, col="red", lty=2);
if(length(outliersID)>0){
points(ellipsis$x[outliersID], ellipsis$y[outliersID], pch=16);
text(ellipsis$x[outliersID], ellipsis$y[outliersID], labels=outliersID, pos=(ellipsis$y[outliersID]>0)*2+1);
}
if(lowess){
# Remove NAs
if(any(is.na(ellipsis$x))){
ellipsis$y <- ellipsis$y[!is.na(ellipsis$x)];
ellipsis$x <- ellipsis$x[!is.na(ellipsis$x)];
}
lines(lowess(ellipsis$x, 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));
if(type=="abs"){
ellipsis$y <- abs(as.vector(residuals(x)));
}
else{
ellipsis$y <- as.vector(residuals(x))^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"){
ellipsis$main <- "|Residuals| 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, 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$loss==c("MAEh","TMAE","GTMAE","MACE"))){
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(any(x$loss==c("HAMh","THAM","GTHAM","CHAM"))){
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(!any(names(ellipsis)=="main")){
ellipsis$main <- "QQ plot of normal distribution";
}
do.call(qqnorm, ellipsis);
qqline(ellipsis$y);
}
}
# 7. Basic plot over time
plot5 <- function(x, ...){
ellipsis <- list(...);
ellipsis$actuals <- actuals(x);
if(!is.null(x$holdout)){
yHoldout <- x$holdout;
if(is.zoo(ellipsis$actuals)){
ellipsis$actuals <- zoo(c(as.vector(ellipsis$actuals),as.vector(yHoldout)),
order.by=c(time(ellipsis$actuals),time(yHoldout)));
}
else{
ellipsis$actuals <- ts(c(ellipsis$actuals,yHoldout),
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;
if(!any(x$interval==c("none","n"))){
ellipsis$lower <- x$lower;
ellipsis$upper <- x$upper;
ellipsis$level <- x$level;
}
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(!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";
}
# Get the IDs of outliers and statistic
outliers <- outlierdummy(x, level=level, type=type);
outliersID <- outliers$id;
statistic <- outliers$statistic;
if(!any(names(ellipsis)=="ylim")){
ellipsis$ylim <- range(c(ellipsis$x,statistic),na.rm=TRUE)*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(is.occurrence(x$occurrence)){
points(ellipsis$x);
}
if(length(outliersID)>0){
points(time(ellipsis$x)[outliersID], ellipsis$x[outliersID], pch=16);
text(time(ellipsis$x)[outliersID], ellipsis$x[outliersID], labels=outliersID, pos=(ellipsis$x[outliersID]>0)*2+1);
}
if(lowess){
# Substitute NAs with the mean
if(any(is.na(ellipsis$x))){
ellipsis$x[is.na(ellipsis$x)] <- mean(ellipsis$x, na.rm=TRUE);
}
lines(lowess(c(1:length(ellipsis$x)),ellipsis$x), col="red");
}
abline(h=0, col="grey", lty=2);
abline(h=statistic[1], col="red", lty=2);
abline(h=statistic[2], col="red", lty=2);
polygon(c(1:nobs(x), c(nobs(x):1)),
c(rep(statistic[1],nobs(x)), rep(statistic[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", squared=FALSE, ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="main")){
if(type=="acf"){
if(squared){
ellipsis$main <- "Autocorrelation Function of Squared Residuals";
}
else{
ellipsis$main <- "Autocorrelation Function of Residuals";
}
}
else{
if(squared){
ellipsis$main <- "Partial Autocorrelation Function of Squared 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(squared){
if(type=="acf"){
theValues <- acf(as.vector(residuals(x)^2), plot=FALSE, na.action=na.pass)
}
else{
theValues <- pacf(as.vector(residuals(x)^2), plot=FALSE, na.action=na.pass);
}
}
else{
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 <- switch(type,
"acf"=theValues$acf[-1],
"pacf"=theValues$acf);
statistic <- 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=statistic, col="red", lty=2);
if(any(ellipsis$x>statistic[2] | ellipsis$x<statistic[1])){
outliersID <- which(ellipsis$x >statistic[2] | ellipsis$x <statistic[1]);
points(outliersID, ellipsis$x[outliersID], pch=16);
text(outliersID, ellipsis$x[outliersID], labels=outliersID, pos=(ellipsis$x[outliersID]>0)*2+1);
}
}
# 12. Plot of states
plot8 <- function(x, ...){
parDefault <- par(no.readonly=TRUE);
on.exit(par(parDefault));
smoothType <- smoothType(x);
if(smoothType=="ETS"){
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.");
}
}
else if(smoothType=="CMA"){
ellipsis$actuals <- actuals(x);
ellipsis$forecast <- x$forecast;
ellipsis$fitted <- x$fitted;
ellipsis$legend <- FALSE;
ellipsis$vline <- FALSE;
if(is.null(ellipsis$main)){
ellipsis$main <- x$model;
}
do.call(graphmaker, ellipsis);
}
else{
if(any(unlist(gregexpr("combine",x$model))!=-1)){
# If we did combinations, we cannot do anything
message("Combination of models was done. Sorry, but there is nothing to plot.");
}
else{
statesNames <- c(colnames(x$states),"residuals");
x$states <- cbind(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))];
do.call(plot, ellipsis);
}
}
else{
if(is.null(ellipsis$main)){
ellipsis$main <- paste0("States of ",x$model);
}
ellipsis$x <- x$states;
do.call(plot, ellipsis);
}
}
}
}
# 13 and 14. Fitted vs (std. Residuals)^2 or Fitted vs |std. Residuals|
plot9 <- function(x, type="abs", ...){
ellipsis <- list(...);
ellipsis$x <- as.vector(fitted(x));
ellipsis$y <- as.vector(rstandard(x));
if(any(x$distribution==c("dinvgauss","dgamma"))){
ellipsis$y[] <- log(ellipsis$y);
}
if(type=="abs"){
ellipsis$y[] <- abs(ellipsis$y);
}
else{
ellipsis$y[] <- ellipsis$y^2;
}
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)];
}
if(!any(names(ellipsis)=="main")){
if(type=="abs"){
if(any(x$distribution==c("dinvgauss","dgamma","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$main <- paste0("|log(Standardised Residuals)| vs Fitted");
}
else{
ellipsis$main <- "|Standardised Residuals| vs Fitted";
}
}
else{
if(any(x$distribution==c("dinvgauss","dgamma","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$main <- paste0("log(Standardised Residuals)^2 vs Fitted");
}
else{
ellipsis$main <- "Standardised Residuals^2 vs Fitted";
}
}
}
if(!any(names(ellipsis)=="xlab")){
ellipsis$xlab <- "Fitted";
}
if(!any(names(ellipsis)=="ylab")){
if(type=="abs"){
if(any(x$distribution==c("dinvgauss","dgamma","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$ylab <- "|log(Standardised Residuals)|";
}
else{
ellipsis$ylab <- "|Standardised Residuals|";
}
}
else{
if(any(x$distribution==c("dinvgauss","dgamma","dlnorm","dllaplace","dls","dlgnorm"))){
ellipsis$ylab <- "log(Standardised Residuals)^2";
}
else{
ellipsis$ylab <- "Standardised 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");
}
}
# Do plots
for(i in which){
if(any(i==1)){
plot1(x, ...);
}
else if(any(i==2)){
plot2(x, ...);
}
else if(any(i==3)){
plot2(x, "rstudent", ...);
}
else if(any(i==4)){
plot3(x, ...);
}
else if(any(i==5)){
plot3(x, type="squared", ...);
}
else if(any(i==6)){
plot4(x, ...);
}
else if(any(i==7)){
plot5(x, ...);
}
else if(any(i==8)){
plot6(x, ...);
}
else if(any(i==9)){
plot6(x, "rstudent", ...);
}
else if(any(i==10)){
plot7(x, type="acf", ...);
}
else if(any(i==11)){
plot7(x, type="pacf", ...);
}
else if(any(i==12)){
plot8(x, ...);
}
else if(any(i==13)){
plot9(x, ...);
}
else if(any(i==14)){
plot9(x, type="squared", ...);
}
else if(any(i==15)){
plot7(x, type="acf", squared=TRUE, ...);
}
else if(any(i==16)){
plot7(x, type="pacf", squared=TRUE, ...);
}
}
}
#' @export
plot.smoothC <- function(x, ...){
graphmaker(actuals(x), x$forecast, x$fitted, x$lower, x$upper, x$level,
main="Combined smooth forecasts");
}
#' @export
plot.smooth.sim <- function(x, ...){
ellipsis <- list(...);
if(is.null(ellipsis$main)){
ellipsis$main <- x$model;
}
if(is.null(dim(x$data))){
nsim <- 1;
}
else{
nsim <- dim(x$data)[2];
}
if(nsim==1){
if(is.null(ellipsis$ylab)){
ellipsis$ylab <- "Data";
}
ellipsis$x <- x$data;
do.call(plot, ellipsis);
}
else{
randomNumber <- ceiling(runif(1,1,nsim));
message(paste0("You have generated ",nsim," time series. Not sure which of them to plot.\n",
"Please use plot(ourSimulation$data[,k]) instead. Plotting randomly selected series N",randomNumber,"."));
if(is.null(ellipsis$ylab)){
ellipsis$ylab <- paste0("Series N",randomNumber);
}
ellipsis$x <- x$data[,randomNumber];
do.call(plot, ellipsis);
}
}
#' @method plot smooth.forecast
#' @export
plot.smooth.forecast <- function(x, ...){
yActuals <- actuals(x$model);
if(!is.null(x$model$holdout)){
yActuals <- ts(c(yActuals,x$model$holdout), start=start(yActuals), frequency=frequency(yActuals));
yActuals <- window(yActuals, start=start(yActuals), end=min(tail(time(x$mean),1),tail(time(yActuals),1)));
}
if(!any(x$interval==c("none","n"))){
graphmaker(yActuals,x$mean,fitted(x$model),x$lower,x$upper,x$level,...);
}
else{
graphmaker(yActuals,x$mean,fitted(x$model),...);
}
}
#' @export
plot.oes <- function(x, which=7, ...){
# This is needed, because diagnostics doesn't make sense in case of oes
plot.smooth(x, which=which, ...);
}
#' @export
plot.oes.sim <- function(x, ...){
ellipsis <- list(...);
if(is.null(ellipsis$main)){
ellipsis$main <- x$model;
}
if(is.null(dim(x$ot))){
nsim <- 1;
}
else{
nsim <- dim(x$ot)[2];
}
if(nsim==1){
if(is.null(ellipsis$ylab)){
ellipsis$ylab <- "Data";
}
ellipsis$x <- x$probability;
do.call(plot, ellipsis);
}
else{
randomNumber <- ceiling(runif(1,1,nsim));
message(paste0("You have generated ",nsim," time series. Not sure which of them to plot.\n",
"Please use plot(ourSimulation$probability[,k]) instead. Plotting randomly selected series N",randomNumber,"."));
if(is.null(ellipsis$ylab)){
ellipsis$ylab <- paste0("Series N",randomNumber);
}
ellipsis$x <- x$probability[,randomNumber];
do.call(plot, ellipsis);
}
}
#### Prints of smooth ####
#' @export
print.smooth <- function(x, ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="digits")){
digits <- 4;
}
else{
digits <- ellipsis$digits;
}
smoothType <- smoothType(x);
if(!is.list(x$model)){
if(smoothType=="CMA"){
holdout <- FALSE;
interval <- FALSE;
cumulative <- FALSE;
}
else{
holdout <- any(!is.na(x$holdout));
interval <- any(!is.na(x$lower));
cumulative <- x$cumulative;
}
}
else{
holdout <- any(!is.na(x$holdout));
interval <- any(!is.na(x$lower));
cumulative <- x$cumulative;
}
if(all(holdout,interval)){
if(!cumulative){
insideinterval <- sum((x$holdout <= x$upper) & (x$holdout >= x$lower)) / length(x$forecast) * 100;
}
else{
insideinterval <- NULL;
}
}
else{
insideinterval <- NULL;
}
intervalType <- x$interval;
if(!is.null(x$model)){
if(!is.list(x$model)){
if(any(smoothType==c("SMA","CMA"))){
x$probability <- 1;
x$initialType <- "b";
occurrence <- "n";
}
else if(smoothType=="ETS"){
# If cumulative forecast and Etype=="M", report that this was "parameteric" interval
if(cumulative & substr(modelType(x),1,1)=="M"){
intervalType <- "p";
}
}
}
}
if(is.occurrence(x$occurrence)){
occurrence <- x$occurrence$occurrence;
}
else{
occurrence <- "n";
}
ssOutput(x$timeElapsed, x$model, persistence=x$persistence, transition=x$transition, measurement=x$measurement,
phi=x$phi, ARterms=x$AR, MAterms=x$MA, constant=x$constant, a=x$a, b=x$b,initialType=x$initialType,
nParam=x$nParam, s2=x$s2, hadxreg=!is.null(x$xreg), wentwild=FALSE,
loss=x$loss, cfObjective=x$lossValue, interval=interval, cumulative=cumulative,
intervalType=intervalType, level=x$level, ICs=x$ICs,
holdout=holdout, insideinterval=insideinterval, errormeasures=x$accuracy,
occurrence=occurrence, obs=nobs(x), digits=digits);
}
#' @export
print.smooth.sim <- function(x, ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="digits")){
digits <- 4;
}
else{
digits <- ellipsis$digits;
}
smoothType <- smoothType(x);
if(is.null(dim(x$data))){
nsim <- 1
}
else{
nsim <- dim(x$data)[2]
}
cat(paste0("Data generated from: ",x$model,"\n"));
cat(paste0("Number of generated series: ",nsim,"\n"));
if(nsim==1){
if(smoothType=="ETS"){
cat(paste0("Persistence vector: \n"));
xPersistence <- as.vector(x$persistence);
names(xPersistence) <- rownames(x$persistence);
print(round(xPersistence,digits));
if(x$phi!=1){
cat(paste0("Phi: ",x$phi,"\n"));
}
if(any(x$occurrence!=1)){
cat(paste0("The data is produced based on an occurrence model.\n"));
}
cat(paste0("True likelihood: ",round(x$logLik,digits),"\n"));
}
else if(smoothType=="ARIMA"){
ar.orders <- orders(x)$ar;
i.orders <- orders(x)$i;
ma.orders <- orders(x)$ma;
lags <- lags(x);
# AR terms
if(any(ar.orders!=0)){
ARterms <- matrix(0,max(ar.orders),sum(ar.orders!=0),
dimnames=list(paste0("AR(",c(1:max(ar.orders)),")"),
paste0("Lag ",lags[ar.orders!=0])));
}
else{
ARterms <- matrix(0,1,1);
}
# Differences
if(any(i.orders!=0)){
Iterms <- matrix(0,1,length(i.orders),
dimnames=list("I(...)",paste0("Lag ",lags)));
Iterms[,] <- i.orders;
}
else{
Iterms <- 0;
}
# MA terms
if(any(ma.orders!=0)){
MAterms <- matrix(0,max(ma.orders),sum(ma.orders!=0),
dimnames=list(paste0("MA(",c(1:max(ma.orders)),")"),
paste0("Lag ",lags[ma.orders!=0])));
}
else{
MAterms <- matrix(0,1,1);
}
n.coef <- ar.coef <- ma.coef <- 0;
ar.i <- ma.i <- 1;
for(i in 1:length(ar.orders)){
if(ar.orders[i]!=0){
ARterms[1:ar.orders[i],ar.i] <- x$AR[ar.coef+(1:ar.orders[i])];
ar.coef <- ar.coef + ar.orders[i];
ar.i <- ar.i + 1;
}
if(ma.orders[i]!=0){
MAterms[1:ma.orders[i],ma.i] <- x$MA[ma.coef+(1:ma.orders[i])];
ma.coef <- ma.coef + ma.orders[i];
ma.i <- ma.i + 1;
}
}
if(!is.null(x$AR)){
cat(paste0("AR parameters: \n"));
print(round(ARterms,digits));
}
if(!is.null(x$MA)){
cat(paste0("MA parameters: \n"));
print(round(MAterms,digits));
}
if(!is.na(x$constant)){
cat(paste0("Constant value: ",round(x$constant,digits),"\n"));
}
cat(paste0("True likelihood: ",round(x$logLik,digits),"\n"));
}
else if(smoothType=="CES"){
cat(paste0("Smoothing parameter a: ",round(x$a,digits),"\n"));
if(!is.null(x$b)){
if(is.complex(x$b)){
cat(paste0("Smoothing parameter b: ",round(x$b,digits),"\n"));
}
else{
cat(paste0("Smoothing parameter b: ",round(x$b,digits),"\n"));
}
}
cat(paste0("True likelihood: ",round(x$logLik,digits),"\n"));
}
else if(smoothType=="SMA"){
cat(paste0("True likelihood: ",round(x$logLik,digits),"\n"));
}
}
}
#' @export
print.smooth.forecast <- function(x, ...){
if(!any(x$interval==c("none","n"))){
level <- x$level;
if(level>1){
level <- level/100;
}
if(x$side=="both"){
output <- cbind(x$mean,x$lower,x$upper);
colnames(output) <- c("Point forecast",
paste0("Lower bound (",(1-level)/2*100,"%)"),
paste0("Upper bound (",(1+level)/2*100,"%)"));
}
else if(x$side=="upper"){
output <- cbind(x$mean,x$upper);
colnames(output) <- c("Point forecast",
paste0("Upper bound (",level*100,"%)"));
}
else if(x$side=="lower"){
output <- cbind(x$mean,x$lower);
colnames(output) <- c("Point forecast",
paste0("Lower bound (",level*100,"%)"));
}
}
else{
output <- x$mean;
}
print(output);
}
#' @export
print.oes <- function(x, ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="digits")){
digits <- 4;
}
else{
digits <- ellipsis$digits;
}
occurrence <- x$occurrence
if(occurrence=="general"){
occurrence <- "General";
}
else if(occurrence=="direct"){
occurrence <- "Direct probability";
}
else if(occurrence=="fixed"){
occurrence <- "Fixed probability";
}
else if(occurrence=="inverse-odds-ratio"){
occurrence <- "Inverse odds ratio";
}
else if(occurrence=="odds-ratio"){
occurrence <- "Odds ratio";
}
else{
occurrence <- "None";
}
ICs <- round(c(AIC(x),AICc(x),BIC(x),BICc(x)),digits);
names(ICs) <- c("AIC","AICc","BIC","BICc");
cat(paste0("Occurrence state space model estimated: ",occurrence,"\n"));
if(!is.null(x$model)){
cat(paste0("Underlying ETS model: ",x$model,"\n"));
}
if(!is.null(x$persistence) && (x$occurrence!="fixed")){
cat("Smoothing parameters:\n");
print(round(x$persistence[,1],digits));
}
if(!is.null(x$initial)){
cat("Vector of initials:\n");
print(round(x$initial,digits));
}
if(!is.null(sigma(x))){
cat("\nError standard deviation: "); cat(round(sigma(x),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));
cat("\nInformation criteria: \n");
print(ICs);
}
#' @export
print.oes.sim <- function(x, ...){
ellipsis <- list(...);
if(!any(names(ellipsis)=="digits")){
digits <- 4;
}
else{
digits <- ellipsis$digits;
}
if(is.null(dim(x$ot))){
nsim <- 1;
obs <- length(x$ot)
}
else{
nsim <- dim(x$ot)[2];
obs <- dim(x$ot)[1];
}
cat(paste0("Data generated from: ",x$model,"\n"));
cat(paste0("Number of generated series: ",nsim,"\n"));
cat(paste0("Number of observations in each series: ",obs,"\n"));
if(nsim==1){
cat(paste0("True likelihood: ",round(x$logLik,digits),"\n"));
}
}
#### Residuals for provided object ####
#' @export
residuals.smooth <- function(object, ...){
if(errorType(object)=="A"){
return(object$residuals);
}
else{
return(log(1+object$residuals));
}
}
#' @importFrom stats rstandard
#' @export
rstandard.smooth <- function(model, ...){
obs <- nobs(model);
df <- obs - nparam(model);
# If this is an occurrence model, then only modify the non-zero obs
if(is.occurrence(model$occurrence)){
residsToGo <- (actuals(model$occurrence)!=0);
}
else{
residsToGo <- rep(TRUE,obs);
}
errors <- residuals(model, ...);
errors[] <- (errors - mean(errors[residsToGo], na.rm=TRUE)) / sqrt(sigma(model)^2 * obs / df);
# Fill in values with NAs if there is occurrence model
if(is.occurrence(model$occurrence)){
errors[!residsToGo] <- NA;
}
return(errors);
}
#' @importFrom stats rstudent
#' @export
rstudent.smooth <- function(model, ...){
obs <- nobs(model);
df <- obs - nparam(model) - 1;
# If this is an occurrence model, then only modify the non-zero obs
if(is.occurrence(model$occurrence)){
residsToGo <- (actuals(model$occurrence)!=0);
}
else{
residsToGo <- rep(TRUE,obs);
}
rstudentised <- errors <- residuals(model, ...);
errors[] <- errors - mean(errors, na.rm=TRUE);
# Prepare the residuals
if(errorType(model)=="M"){
for(i in which(residsToGo)){
rstudentised[i] <- errors[i] / sqrt(sum(errors[-i]^2, na.rm=TRUE) / df);
}
}
else{
for(i in which(residsToGo)){
rstudentised[i] <- errors[i] / sqrt(sum(errors[-i]^2, na.rm=TRUE) / df);
}
}
# Fill in values with NAs if there is occurrence model
if(is.occurrence(model$occurrence)){
rstudentised[!residsToGo] <- NA;
}
return(rstudentised);
}
#' @importFrom greybox outlierdummy
#' @export
outlierdummy.smooth <- 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$loss,
"MAE"=,"MAEh"=,"MACE"=,"TMAE"=,"GTMAE"=qlaplace(c((1-level)/2, (1+level)/2), 0, 1),
"HAM"=,"HAMh"=,"CHAM"=,"THAM"=,"GTHAM"=qs(c((1-level)/2, (1+level)/2), 0, 1),
qnorm(c((1-level)/2, (1+level)/2), 0, 1));
outliersID <- which(errors>statistic[2] | errors<statistic[1]);
outliersNumber <- length(outliersID);
if(outliersNumber>0){
outliers <- ts(matrix(0, nobs(object), outliersNumber,
dimnames=list(NULL,
paste0("outlier",c(1:outliersNumber)))),
start=start(actuals(object)), frequency=frequency(actuals(object)));
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"));
}
#' @export
rmultistep.smooth <- function(object, h=10, ...){
yInSample <- actuals(object);
model <- modelType(object);
Etype <- errorType(object);
Ttype <- substr(model,2,2);
Stype <- substr(model,nchar(model),nchar(model));
lagsModel <- modelLags(object);
obsInSample <- nobs(object);
matxt <- object$xreg;
if(is.null(object$xreg)){
matxt <- matrix(1,obsInSample,1);
}
nXreg <- ncol(matxt);
if(is.null(object$xreg)){
matat <- matrix(1,obsInSample,1);
}
else{
matat <- matrix(object$initialX,obsInSample,nXreg,byrow=TRUE);
}
matFX <- diag(nXreg);
dataStart <- start(yInSample);
dataFreq <- frequency(yInSample);
errors.mat <- ts(errorerwrap(object$states, object$transition, object$measurement, matrix(yInSample,obsInSample,1),
h, Etype, Ttype, Stype, lagsModel,
matxt, matat, matFX, matrix(1,obsInSample,1)),
start=dataStart,frequency=dataFreq);
colnames(errors.mat) <- paste0("Error",c(1:h));
return(errors.mat);
}
#### Simulate data using provided object ####
#' @importFrom utils tail
#' @export
simulate.smooth <- function(object, nsim=1, seed=NULL, obs=NULL, ...){
ellipsis <- list(...);
smoothType <- smoothType(object);
if(is.null(obs)){
obs <- nobs(object);
}
if(!is.null(seed)){
set.seed(seed);
}
# Start a list of arguments
args <- vector("list",0);
loss <- object$loss;
if(any(loss==c("MAE","MAEh","TMAE","GTMAE","MACE"))){
randomizer <- "rlaplace";
if(!is.null(ellipsis$mu)){
args$mu <- ellipsis$mu;
}
else{
args$mu <- 0;
}
if(!is.null(ellipsis$scale)){
args$scale <- ellipsis$scale;
}
else{
args$scale <- mean(abs(residuals(object)));
}
}
else if(any(loss==c("HAM","HAMh","THAM","GTHAM","CHAM"))){
randomizer <- "rs";
if(!is.null(ellipsis$mu)){
args$mu <- ellipsis$mu;
}
else{
args$mu <- 0;
}
if(!is.null(ellipsis$scale)){
args$scale <- ellipsis$scale;
}
else{
args$scale <- mean(sqrt(abs(residuals(object))));
}
}
else{
randomizer <- "rnorm";
if(!is.null(ellipsis$mean)){
args$mean <- ellipsis$mean;
}
else{
args$mean <- 0;
}
if(!is.null(ellipsis$sd)){
args$sd <- ellipsis$sd;
}
else{
args$sd <- sigma(object);
}
}
args$randomizer <- randomizer;
args$frequency <- frequency(actuals(object));
args$obs <- obs;
args$nsim <- nsim;
args$initial <- object$initial;
# If this is an occurrence model, use the fitted values for the probabilities
if(is.occurrence(object$occurrence)){
args$probability <- fitted(object$occurrence);
}
else{
args$probability <- 1;
}
if(smoothType=="ETS"){
model <- modelType(object);
if(any(unlist(gregexpr("C",model))==-1)){
args$model <- model;
args$phi <- object$phi;
args$persistence <- object$persistence;
if(is.list(args$initial)){
args$initialSeason <- object$initial$seasonal;
args$initial <- unlist(args$initial[c("level","trend")]);
}
if(errorType(object)=="M" && args$mean==0){
args$mean <- 1;
}
simulatedData <- do.call("sim.es",args);
}
else{
message("Sorry, but we cannot simulate data from combined model.");
simulatedData <- NA;
}
}
else if(smoothType=="ARIMA"){
if(any(unlist(gregexpr("combine",object$model))==-1)){
args$orders <- orders(object);
args$lags <- lags(object);
args$AR <- object$AR;
args$MA <- object$MA;
args$constant <- object$constant;
simulatedData <- do.call("sim.ssarima",args);
}
else{
message("Sorry, but we cannot simulate data from combined model.");
simulatedData <- NA;
}
}
else if(smoothType=="CES"){
args$seasonality <- modelType(object);
args$a <- object$a;
args$b <- object$b;
simulatedData <- do.call("sim.ces",args);
}
else if(smoothType=="GUM"){
args$orders <- orders(object);
args$lags <- lags(object);
args$measurement <- object$measurement;
args$transition <- object$transition;
args$persistence <- object$persistence;
simulatedData <- do.call("sim.gum",args);
}
else if(smoothType=="SMA"){
args$order <- orders(object)[1];
simulatedData <- do.call("sim.sma",args);
}
else{
model <- substring(object$model,1,unlist(gregexpr("\\(",object$model))[1]-1);
message(paste0("Sorry, but simulate is not yet available for the model ",model,"."));
simulatedData <- NA;
}
return(simulatedData);
}
#### Type of smooth model. Internal function ####
smoothType <- function(object, ...){
if(!is.list(object$model)){
if(gregexpr("ETS",object$model)!=-1){
smoothType <- "ETS";
}
else if(gregexpr("CES",object$model)!=-1){
smoothType <- "CES";
}
else if(gregexpr("ARIMA",object$model)!=-1){
smoothType <- "ARIMA";
}
else if(gregexpr("GUM",object$model)!=-1){
smoothType <- "GUM";
}
else if(gregexpr("SMA",object$model)!=-1){
smoothType <- "SMA";
}
else if(gregexpr("CMA",object$model)!=-1){
smoothType <- "CMA";
}
else if(gregexpr("VES",object$model)!=-1){
smoothType <- "VES";
}
else{
smoothType <- NA;
}
}
else{
smoothType <- "smoothCombine";
}
return(smoothType);
}
#### Summary of objects ####
#' @export
summary.smooth <- function(object, ...){
print(object);
}
#' @export
summary.smooth.forecast <- function(object, ...){
print(object);
}
#### Accuracy measures ####
#' @importFrom generics accuracy
#' @export
generics::accuracy
#' Error measures for an estimated model
#'
#' Function produces error measures for the provided object and the holdout values of the
#' response variable. Note that instead of parameters \code{x}, \code{test}, the function
#' accepts the vector of values in \code{holdout}. Also, the parameters \code{d} and \code{D}
#' are not supported - MASE is always calculated via division by first differences.
#'
#' The function is a wrapper for the \link[greybox]{measures} function and is implemented
#' for convenience.
#'
#' @template ssAuthor
#'
#' @param object The estimated model or a forecast from the estimated model generated via
#' either \code{predict()} or \code{forecast()} functions.
#' @param holdout The vector of values of the response variable in the holdout (test) set.
#' If not provided, then the function will return the in-sample error measures. If the
#' \code{holdout=TRUE} parameter was used in the estimation of a model, the holdout values
#' will be extracted automatically.
#' @param ... Other variables passed to the \code{forecast()} function (e.g. \code{newdata}).
#' @examples
#'
#' y <- rnorm(100, 100, 10)
#' ourModel <- adam(y, holdout=TRUE, h=10)
#' accuracy(ourModel)
#'
#' @rdname accuracy
#' @export
accuracy.smooth <- function(object, holdout=NULL, ...){
if(is.null(object$holdout)){
if(is.null(holdout)){
return(measures(actuals(object), fitted(object), actuals(object)));
}
else{
h <- length(holdout);
return(measures(holdout, forecast(object, h=h, ...)$mean, actuals(object)));
}
}
else{
return(object$accuracy);
}
}
#' @rdname accuracy
#' @export
accuracy.smooth.forecast <- function(object, holdout=NULL, ...){
if(is.null(holdout)){
if(is.null(object$model$holdout)){
return(measures(actuals(object), fitted(object), actuals(object)));
}
else{
return(object$model$accuracy);
}
}
else{
h <- length(holdout);
return(measures(holdout, object$mean, actuals(object)));
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.