Nothing
utils::globalVariables(c("silentText","silentGraph","silentLegend","initialType","ar.orders","i.orders","ma.orders"));
#' State Space ARIMA
#'
#' Function selects the best State Space ARIMA based on information criteria,
#' using fancy branch and bound mechanism. The resulting model can be not
#' optimal in IC meaning, but it is usually reasonable.
#'
#' The function constructs bunch of ARIMAs in Single Source of Error
#' state space form (see \link[smooth]{ssarima} documentation) and selects the
#' best one based on information criterion. The mechanism is described in
#' Svetunkov & Boylan (2019).
#'
#' Due to the flexibility of the model, multiple seasonalities can be used. For
#' example, something crazy like this can be constructed:
#' SARIMA(1,1,1)(0,1,1)[24](2,0,1)[24*7](0,0,1)[24*30], but the estimation may
#' take a lot of time... It is recommended to use \link[smooth]{auto.msarima} in
#' cases with more than one seasonality and high frequencies.
#'
#' For some more information about the model and its implementation, see the
#' vignette: \code{vignette("ssarima","smooth")}
#'
#' @template ssBasicParam
#' @template ssAdvancedParam
#' @template ssXregParam
#' @template ssIntervals
#' @template ssInitialParam
#' @template ssAuthor
#' @template ssKeywords
#'
#' @template ssGeneralRef
#' @template ssIntermittentRef
#' @template ssARIMARef
#'
#' @param orders List of maximum orders to check, containing vector variables
#' \code{ar}, \code{i} and \code{ma}. If a variable is not provided in the
#' list, then it is assumed to be equal to zero. At least one variable should
#' have the same length as \code{lags}.
#' @param lags Defines lags for the corresponding orders (see examples). The
#' length of \code{lags} must correspond to the length of \code{orders}. There
#' is no restrictions on the length of \code{lags} vector.
#' @param combine If \code{TRUE}, then resulting ARIMA is combined using AIC
#' weights.
#' @param fast If \code{TRUE}, then some of the orders of ARIMA are
#' skipped. This is not advised for models with \code{lags} greater than 12.
#' @param constant If \code{NULL}, then the function will check if constant is
#' needed. if \code{TRUE}, then constant is forced in the model. Otherwise
#' constant is not used.
#' @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 parameters of the model. Maximum orders to
#' check can also be specified separately, however \code{orders} variable must
#' be set to \code{NULL}: \code{ar.orders} - Maximum order of AR term. Can be
#' vector, defining max orders of AR, SAR etc. \code{i.orders} - Maximum order
#' of I. Can be vector, defining max orders of I, SI etc. \code{ma.orders} -
#' Maximum order of MA term. Can be vector, defining max orders of MA, SMA etc.
#' @return Object of class "smooth" is returned. See \link[smooth]{ssarima} for
#' details.
#' @seealso \code{\link[smooth]{es}, \link[smooth]{ces},
#' \link[smooth]{sim.es}, \link[smooth]{gum}, \link[smooth]{ssarima}}
#'
#' @examples
#'
#' \donttest{x <- rnorm(118,100,3)}
#'
#' # The best ARIMA for the data
#' \donttest{ourModel <- auto.ssarima(x,orders=list(ar=c(2,1),i=c(1,1),ma=c(2,1)),lags=c(1,12),
#' h=18,holdout=TRUE,interval="np")}
#'
#' # The other one using optimised states
#' \donttest{auto.ssarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
#' initial="o",h=18,holdout=TRUE)}
#'
#' # And now combined ARIMA
#' \donttest{auto.ssarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
#' combine=TRUE,h=18,holdout=TRUE)}
#'
#' \donttest{summary(ourModel)
#' forecast(ourModel)
#' plot(forecast(ourModel))}
#'
#'
#' @export auto.ssarima
auto.ssarima <- function(y, orders=list(ar=c(3,3),i=c(2,1),ma=c(3,3)), lags=c(1,frequency(y)),
combine=FALSE, fast=TRUE, constant=NULL,
initial=c("backcasting","optimal"), ic=c("AICc","AIC","BIC","BICc"),
loss=c("likelihood","MSE","MAE","HAM","MSEh","TMSE","GTMSE","MSCE"),
h=10, holdout=FALSE, cumulative=FALSE,
interval=c("none","parametric","likelihood","semiparametric","nonparametric"), level=0.95,
bounds=c("admissible","none"),
silent=c("all","graph","legend","output","none"),
xreg=NULL, regressors=c("use","select"), initialX=NULL, ...){
# Function estimates several ssarima models and selects the best one using the selected information criterion.
#
# Copyright (C) 2015 - Inf Ivan Svetunkov
# Start measuring the time of calculations
startTime <- Sys.time();
### Depricate the old parameters
ellipsis <- list(...)
ellipsis <- depricator(ellipsis, "xregDo", "regressors");
updateX <- FALSE;
persistenceX <- transitionX <- NULL;
occurrence <- "none";
oesmodel <- "MNN";
# Add all the variables in ellipsis to current environment
list2env(ellipsis,environment());
if(!is.null(orders)){
arMax <- orders$ar;
iMax <- orders$i;
maMax <- orders$ma;
}
# If orders are provided in ellipsis via arMax, write them down.
if(exists("ar.orders",inherits=FALSE)){
if(is.null(ar.orders)){
arMax <- 0;
}
else{
arMax <- ar.orders;
}
}
else{
if(is.null(orders)){
arMax <- 0;
}
}
if(exists("i.orders",inherits=FALSE)){
if(is.null(i.orders)){
iMax <- 0;
}
else{
iMax <- i.orders;
}
}
else{
if(is.null(orders)){
iMax <- 0;
}
}
if(exists("ma.orders",inherits=FALSE)){
if(is.null(ma.orders)){
maMax <- 0;
}
else{
maMax <- ma.orders
}
}
else{
if(is.null(orders)){
maMax <- 0;
}
}
##### Set environment for ssInput and make all the checks #####
environment(ssAutoInput) <- environment();
ssAutoInput("auto.ssarima",ParentEnvironment=environment());
if(is.null(constant)){
constantCheck <- TRUE;
constantValue <- TRUE;
}
else{
if(is.logical(constant)){
constantCheck <- FALSE;
constantValue <- constant;
}
else{
constant <- NULL;
constantCheck <- TRUE;
constantValue <- TRUE;
warning("Strange value of constant parameter. We changed it to default value.");
}
}
if(any(is.complex(c(arMax,iMax,maMax,lags)))){
stop("Come on! Be serious! This is ARIMA, not CES!",call.=FALSE);
}
if(any(c(arMax,iMax,maMax)<0)){
stop("Funny guy! How am I gonna construct a model with negative order?",call.=FALSE);
}
if(any(c(lags)<0)){
stop("Right! Why don't you try complex lags then, mister smart guy?",call.=FALSE);
}
# If there are zero lags, drop them
if(any(lags==0)){
arMax <- arMax[lags!=0];
iMax <- iMax[lags!=0];
maMax <- maMax[lags!=0];
lags <- lags[lags!=0];
}
# Define maxorder and make all the values look similar (for the polynomials)
maxorder <- max(length(arMax),length(iMax),length(maMax));
if(length(arMax)!=maxorder){
arMax <- c(arMax,rep(0,maxorder-length(arMax)));
}
if(length(iMax)!=maxorder){
iMax <- c(iMax,rep(0,maxorder-length(iMax)));
}
if(length(maMax)!=maxorder){
maMax <- c(maMax,rep(0,maxorder-length(maMax)));
}
# If zeroes are defined as orders for some lags, drop them.
if(any((arMax + iMax + maMax)==0)){
orders2leave <- (arMax + iMax + maMax)!=0;
if(all(!orders2leave)){
orders2leave <- lags==min(lags);
}
arMax <- arMax[orders2leave];
iMax <- iMax[orders2leave];
maMax <- maMax[orders2leave];
lags <- lags[orders2leave];
}
# Get rid of duplicates in lags
if(length(unique(lags))!=length(lags)){
if(dataFreq!=1){
warning(paste0("'lags' variable contains duplicates: (",paste0(lags,collapse=","),"). Getting rid of some of them."),call.=FALSE);
}
lagsNew <- unique(lags);
arMaxNew <- iMaxNew <- maMaxNew <- lagsNew;
for(i in 1:length(lagsNew)){
arMaxNew[i] <- max(arMax[which(lags==lagsNew[i])],na.rm=TRUE);
iMaxNew[i] <- max(iMax[which(lags==lagsNew[i])],na.rm=TRUE);
maMaxNew[i] <- max(maMax[which(lags==lagsNew[i])],na.rm=TRUE);
}
arMax <- arMaxNew;
iMax <- iMaxNew;
maMax <- maMaxNew;
lags <- lagsNew;
}
# Order things, so we would deal with the lowest level of seasonality first
arMax <- arMax[order(lags,decreasing=FALSE)];
iMax <- iMax[order(lags,decreasing=FALSE)];
maMax <- maMax[order(lags,decreasing=FALSE)];
lags <- sort(lags,decreasing=FALSE);
# 1 stands for constant, the other one stands for variance
nParamMax <- (1 + max(arMax %*% lags + iMax %*% lags,maMax %*% lags)
+ sum(arMax) + sum(maMax) + constantCheck);
# Try to figure out if the number of parameters can be tuned in order to fit something smaller on small samples
# Don't try to fix anything if the number of seasonalities is greater than 2
if(length(lags)<=2){
if(obsNonzero <= nParamMax){
armaLength <- length(arMax);
while(obsNonzero <= nParamMax){
if(any(c(arMax[armaLength],maMax[armaLength])>0)){
arMax[armaLength] <- max(0,arMax[armaLength] - 1);
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
if(obsNonzero <= nParamMax){
maMax[armaLength] <- max(0,maMax[armaLength] - 1);
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
}
}
else{
if(armaLength==2){
arMax[1] <- arMax[1] - 1;
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
if(obsNonzero <= nParamMax){
maMax[1] <- maMax[1] - 1;
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
}
}
else{
break;
}
}
if(all(c(arMax,maMax)==0)){
if(iMax[armaLength]>0){
iMax[armaLength] <- max(0,iMax[armaLength] - 1);
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
}
else if(iMax[1]>0){
if(obsNonzero <= nParamMax){
iMax[1] <- max(0,iMax[1] - 1);
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
}
}
else{
break;
}
}
}
nParamMax <- max(arMax %*% lags + iMax %*% lags,maMax %*% lags) + sum(arMax) + sum(maMax) + 1 + 1;
}
}
if(obsNonzero <= nParamMax){
message(paste0("Not enough observations for the reasonable fit. Number of possible parameters is ",
nParamMax," while the number of observations is ",obsNonzero,"!"));
stop("Redefine maximum orders and try again.",call.=FALSE)
}
# 1 stands for constant/no constant, another one stands for ARIMA(0,0,0)
if(all(maMax==0)){
nModels <- prod(iMax + 1) * (1 + sum(arMax)) + constantCheck;
}
else{
nModels <- prod(iMax + 1) * (1 + sum(maMax*(1 + sum(arMax)))) + constantCheck;
}
testModel <- list(NA);
# Array with elements x maxorders x horizon x point/lower/upper
if(combine){
testForecasts <- list(NA);
testFitted <- list(NA);
testICs <- list(NA);
testLevels <- list(NA);
testStates <- list(NA);
testTransition <- list(NA);
testPersistence <- list(NA);
}
ICValue <- 1E+100;
m <- 0;
# constant <- TRUE;
lagsTest <- maTest <- arTest <- rep(0,length(lags));
arBest <- maBest <- iBest <- rep(0,length(lags));
arBestLocal <- maBestLocal <- arBest;
#### Function corrects IC taking number of parameters on previous step ####
icCorrector <- function(icValue, nParam, obsNonzero, nParamNew){
if(ic=="AIC"){
llikelihood <- (2*nParam - icValue)/2;
correction <- 2*nParamNew - 2*llikelihood;
}
else if(ic=="AICc"){
llikelihood <- (2*nParam*obsNonzero/(obsNonzero-nParam-1) - icValue)/2;
correction <- 2*nParamNew*obsNonzero/(obsNonzero-nParamNew-1) - 2*llikelihood;
}
else if(ic=="BIC"){
llikelihood <- (nParam*log(obsNonzero) - icValue)/2;
correction <- nParamNew*log(obsNonzero) - 2*llikelihood;
}
else if(ic=="BICc"){
llikelihood <- ((nParam*log(obsNonzero)*obsNonzero)/(obsNonzero-nParam-1) - icValue)/2;
correction <- (nParamNew*log(obsNonzero)*obsNonzero)/(obsNonzero-nParamNew-1) - 2*llikelihood;
}
return(correction);
}
if(!silentText){
cat("Estimation progress: ");
}
### If for some reason we have model with zeroes for orders, return it.
if(all(c(arMax,iMax,maMax)==0)){
cat("\b\b\b\bDone!\n");
bestModel <- ssarima(y, orders=list(ar=arBest,i=(iBest),ma=(maBest)), lags=(lags),
constant=constantValue, initial=initialType, loss=loss,
h=h, holdout=holdout, cumulative=cumulative,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=xreg, regressors=regressors, initialX=initialX, FI=FI);
return(bestModel);
}
iOrders <- matrix(0,prod(iMax+1),ncol=length(iMax));
##### Loop for differences #####
if(any(iMax!=0)){
# Prepare table with differences
iOrders[,1] <- rep(c(0:iMax[1]),times=prod(iMax[-1]+1));
if(length(iMax)>1){
for(seasLag in 2:length(iMax)){
iOrders[,seasLag] <- rep(c(0:iMax[seasLag]),each=prod(iMax[1:(seasLag-1)]+1))
}
}
}
# Start the loop with differences
for(d in 1:nrow(iOrders)){
m <- m + 1;
if(!silentText){
cat(paste0(rep("\b",nchar(round(m/nModels,2)*100)+1),collapse=""));
cat(paste0(round((m)/nModels,2)*100,"%"));
}
nParamOriginal <- 1;
if(silent[1]=="d"){
cat("I: ");cat(iOrders[d,]);cat(", ");
}
testModel <- ssarima(y, orders=list(ar=0,i=iOrders[d,],ma=0), lags=lags,
constant=constantValue, initial=initialType, loss=loss,
h=h, holdout=holdout, cumulative=cumulative,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=xreg, regressors=regressors, initialX=initialX, FI=FI);
ICValue <- testModel$ICs[ic];
if(combine){
testForecasts[[m]] <- matrix(NA,h,3);
testForecasts[[m]][,1] <- testModel$forecast;
testForecasts[[m]][,2] <- testModel$lower;
testForecasts[[m]][,3] <- testModel$upper;
testFitted[[m]] <- testModel$fitted;
testICs[[m]] <- ICValue;
testLevels[[m]] <- 1;
testStates[[m]] <- testModel$states;
testTransition[[m]] <- testModel$transition;
testPersistence[[m]] <- testModel$persistence;
}
if(silent[1]=="d"){
cat(ICValue); cat("\n");
}
if(m==1){
bestIC <- ICValue;
dataMA <- dataI <- testModel$residuals;
iBest <- iOrders[d,];
bestICAR <- bestICI <- bestICMA <- bestIC;
}
else{
if(ICValue < bestICI){
bestICI <- ICValue;
dataMA <- dataI <- testModel$residuals;
if(ICValue < bestIC){
iBest <- iOrders[d,];
bestIC <- ICValue;
maBest <- arBest <- rep(0,length(arTest));
}
}
else{
if(fast){
m <- m + sum(maMax*(1 + sum(arMax)));
next;
}
}
}
##### Loop for MA #####
if(any(maMax!=0)){
bestICMA <- bestICI;
maBestLocal <- maTest <- rep(0,length(maTest));
for(seasSelectMA in 1:length(lags)){
if(maMax[seasSelectMA]!=0){
for(maSelect in 1:maMax[seasSelectMA]){
m <- m + 1;
if(!silentText){
cat(paste0(rep("\b",nchar(round(m/nModels,2)*100)+1),collapse=""));
cat(paste0(round((m)/nModels,2)*100,"%"));
}
maTest[seasSelectMA] <- maMax[seasSelectMA] - maSelect + 1;
nParamMA <- sum(maTest);
nParamNew <- nParamOriginal + nParamMA;
if(silent[1]=="d"){
cat("MA: ");cat(maTest);cat(", ");
}
testModel <- ssarima(dataI, orders=list(ar=0,i=0,ma=maTest), lags=lags,
constant=FALSE, initial=initialType, loss=loss,
h=h, holdout=FALSE,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=NULL, regressors="use", initialX=initialX, FI=FI);
ICValue <- icCorrector(testModel$ICs[ic], nParamMA, obsNonzero, nParamNew);
if(combine){
testForecasts[[m]] <- matrix(NA,h,3);
testForecasts[[m]][,1] <- testModel$forecast;
testForecasts[[m]][,2] <- testModel$lower;
testForecasts[[m]][,3] <- testModel$upper;
testFitted[[m]] <- testModel$fitted;
testICs[[m]] <- ICValue;
testLevels[[m]] <- 2;
testStates[[m]] <- testModel$states;
testTransition[[m]] <- testModel$transition;
testPersistence[[m]] <- testModel$persistence;
}
if(silent[1]=="d"){
cat(ICValue); cat("\n");
}
if(ICValue < bestICMA){
bestICMA <- ICValue;
maBestLocal <- maTest;
if(ICValue < bestIC){
bestIC <- bestICMA;
iBest <- iOrders[d,];
maBest <- maTest;
arBest <- rep(0,length(arTest));
}
dataMA <- testModel$residuals;
}
else{
if(fast){
m <- m + maTest[seasSelectMA] * (1 + sum(arMax)) - 1;
maTest <- maBestLocal;
break;
}
else{
maTest <- maBestLocal;
}
}
##### Loop for AR #####
if(any(arMax!=0)){
bestICAR <- bestICMA;
arBestLocal <- arTest <- rep(0,length(arTest));
for(seasSelectAR in 1:length(lags)){
lagsTest[seasSelectAR] <- lags[seasSelectAR];
if(arMax[seasSelectAR]!=0){
for(arSelect in 1:arMax[seasSelectAR]){
m <- m + 1;
if(!silentText){
cat(paste0(rep("\b",nchar(round(m/nModels,2)*100)+1),collapse=""));
cat(paste0(round((m)/nModels,2)*100,"%"));
}
arTest[seasSelectAR] <- arMax[seasSelectAR] - arSelect + 1;
nParamAR <- sum(arTest);
nParamNew <- nParamOriginal + nParamMA + nParamAR;
if(silent[1]=="d"){
cat("AR: ");cat(arTest);cat(", ");
}
testModel <- ssarima(dataMA, orders=list(ar=arTest,i=0,ma=0), lags=lags,
constant=FALSE, initial=initialType, loss=loss,
h=h, holdout=FALSE,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=NULL, regressors="use", initialX=initialX, FI=FI);
ICValue <- icCorrector(testModel$ICs[ic], nParamAR, obsNonzero, nParamNew);
if(combine){
testForecasts[[m]] <- matrix(NA,h,3);
testForecasts[[m]][,1] <- testModel$forecast;
testForecasts[[m]][,2] <- testModel$lower;
testForecasts[[m]][,3] <- testModel$upper;
testFitted[[m]] <- testModel$fitted;
testICs[[m]] <- ICValue;
testLevels[[m]] <- 3;
testStates[[m]] <- testModel$states;
testTransition[[m]] <- testModel$transition;
testPersistence[[m]] <- testModel$persistence;
}
if(silent[1]=="d"){
cat(ICValue); cat("\n");
}
if(ICValue < bestICAR){
bestICAR <- ICValue;
arBestLocal <- arTest;
if(ICValue < bestIC){
bestIC <- ICValue;
iBest <- iOrders[d,];
arBest <- arTest;
maBest <- maTest;
}
}
else{
if(fast){
m <- m + arTest[seasSelectAR] - 1;
arTest <- arBestLocal;
break;
}
else{
arTest <- arBestLocal;
}
}
}
}
}
}
}
}
}
}
else{
##### Loop for AR #####
if(any(arMax!=0)){
bestICAR <- bestICMA;
arBestLocal <- arTest <- rep(0,length(arTest));
for(seasSelectAR in 1:length(lags)){
lagsTest[seasSelectAR] <- lags[seasSelectAR];
if(arMax[seasSelectAR]!=0){
for(arSelect in 1:arMax[seasSelectAR]){
m <- m + 1;
if(!silentText){
cat(paste0(rep("\b",nchar(round(m/nModels,2)*100)+1),collapse=""));
cat(paste0(round((m)/nModels,2)*100,"%"));
}
arTest[seasSelectAR] <- arMax[seasSelectAR] - arSelect + 1;
nParamAR <- sum(arTest);
nParamNew <- nParamOriginal + nParamAR;
if(silent[1]=="d"){
cat("AR: ");cat(arTest);cat(", ");
}
testModel <- ssarima(dataMA, orders=list(ar=arTest,i=0,ma=0), lags=lags,
constant=FALSE, initial=initialType, loss=loss,
h=h, holdout=FALSE,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=NULL, regressors="use", initialX=initialX, FI=FI);
ICValue <- icCorrector(testModel$ICs[ic], nParamAR, obsNonzero, nParamNew);
if(combine){
testForecasts[[m]] <- matrix(NA,h,3);
testForecasts[[m]][,1] <- testModel$forecast;
testForecasts[[m]][,2] <- testModel$lower;
testForecasts[[m]][,3] <- testModel$upper;
testFitted[[m]] <- testModel$fitted;
testICs[[m]] <- ICValue;
testLevels[[m]] <- 3;
testStates[[m]] <- testModel$states;
testTransition[[m]] <- testModel$transition;
testPersistence[[m]] <- testModel$persistence;
}
if(silent[1]=="d"){
cat(ICValue); cat("\n");
}
if(ICValue < bestICAR){
bestICAR <- ICValue;
arBestLocal <- arTest;
if(ICValue < bestIC){
bestIC <- ICValue;
iBest <- iOrders[d,];
arBest <- arTest;
maBest <- maTest;
}
}
else{
if(fast){
m <- m + arTest[seasSelectAR] - 1;
arTest <- arBestLocal;
break;
}
else{
arTest <- arBestLocal;
}
}
}
}
}
}
}
}
#### Test the constant ####
if(constantCheck){
m <- m + 1;
if(!silentText){
cat(paste0(rep("\b",nchar(round(m/nModels,2)*100)+1),collapse=""));
cat(paste0(round((m)/nModels,2)*100,"%"));
}
if(any(c(arBest,iBest,maBest)!=0)){
testModel <- ssarima(y, orders=list(ar=(arBest),i=(iBest),ma=(maBest)), lags=(lags),
constant=FALSE, initial=initialType, loss=loss,
h=h, holdout=holdout, cumulative=cumulative,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=xreg, regressors=regressors, initialX=initialX, FI=FI);
ICValue <- testModel$ICs[ic];
if(combine){
testForecasts[[m]] <- matrix(NA,h,3);
testForecasts[[m]][,1] <- testModel$forecast;
testForecasts[[m]][,2] <- testModel$lower;
testForecasts[[m]][,3] <- testModel$upper;
testFitted[[m]] <- testModel$fitted;
testICs[[m]] <- ICValue;
testLevels[[m]] <- 1;
testStates[[m]] <- testModel$states;
testTransition[[m]] <- testModel$transition;
testPersistence[[m]] <- testModel$persistence;
}
# cat("Constant: ");print(ICValue);
if(ICValue < bestIC){
bestModel <- testModel;
constantValue <- FALSE;
}
else{
constantValue <- TRUE;
}
}
}
if(combine){
testICs <- unlist(testICs);
testLevels <- unlist(testLevels);
testForecasts <- array(unlist(testForecasts),c(h,3,length(testICs)));
testFitted <- matrix(unlist(testFitted),ncol=length(testICs));
icWeights <- exp(-0.5*(testICs-min(testICs)))/sum(exp(-0.5*(testICs-min(testICs))));
testForecastsNew <- testForecasts;
testFittedNew <- testFitted;
for(i in 1:length(testLevels)){
if(testLevels[i]==1){
j <- i;
}
else if(testLevels[i]==2){
k <- i;
testForecastsNew[,,i] <- testForecasts[,,j] + testForecasts[,,i];
testFittedNew[,i] <- testFitted[,j] + testFitted[,i];
}
else if(testLevels[i]==3){
testForecastsNew[,,i] <- testForecasts[,,j] + testForecasts[,,k] + testForecasts[,,i];
testFittedNew[,i] <- testFitted[,j] + testFitted[,k] + testFitted[,i];
}
}
yForecast <- ts(testForecastsNew[,1,] %*% icWeights,start=yForecastStart,frequency=dataFreq);
yLower <- ts(testForecastsNew[,2,] %*% icWeights,start=yForecastStart,frequency=dataFreq);
yUpper <- ts(testForecastsNew[,3,] %*% icWeights,start=yForecastStart,frequency=dataFreq);
yFitted <- ts(testFittedNew %*% icWeights,start=dataStart,frequency=dataFreq);
modelname <- "ARIMA combined";
errors <- ts(yInSample-c(yFitted),start=dataStart,frequency=dataFreq);
yHoldout <- ts(y[(obsNonzero+1):obsAll],start=yForecastStart,frequency=dataFreq);
s2 <- mean(errors^2);
errormeasures <- measures(yHoldout,yForecast,yInSample);
ICs <- c(t(testICs) %*% icWeights);
names(ICs) <- ic;
bestModel <- list(model=modelname,timeElapsed=Sys.time()-startTime,
initialType=initialType,
fitted=yFitted,forecast=yForecast,cumulative=cumulative,
lower=yLower,upper=yUpper,residuals=errors,s2=s2,interval=intervalType,level=level,
y=y,holdout=yHoldout,
xreg=xreg, regressors=regressors, initialX=initialX,
ICs=ICs,ICw=icWeights,lossValue=NULL,loss=loss,accuracy=errormeasures);
bestModel <- structure(bestModel,class="smooth");
}
else{
#### Reestimate the best model in order to get rid of bias ####
bestModel <- ssarima(y, orders=list(ar=(arBest),i=(iBest),ma=(maBest)), lags=(lags),
constant=constantValue, initial=initialType, loss=loss,
h=h, holdout=holdout, cumulative=cumulative,
interval=intervalType, level=level,
bounds=bounds, silent=TRUE,
xreg=xreg, regressors=regressors, initialX=initialX, FI=FI);
yFitted <- bestModel$fitted;
yForecast <- bestModel$forecast;
yUpper <- bestModel$upper;
yLower <- bestModel$lower;
modelname <- bestModel$model;
bestModel$timeElapsed <- Sys.time()-startTime;
}
if(!silentText){
cat("... Done! \n");
}
##### Make a plot #####
if(!silentGraph){
yForecastNew <- yForecast;
yUpperNew <- yUpper;
yLowerNew <- yLower;
if(cumulative){
yForecastNew <- ts(rep(yForecast/h,h),start=yForecastStart,frequency=dataFreq)
if(interval){
yUpperNew <- ts(rep(yUpper/h,h),start=yForecastStart,frequency=dataFreq)
yLowerNew <- ts(rep(yLower/h,h),start=yForecastStart,frequency=dataFreq)
}
}
if(interval){
graphmaker(actuals=y,forecast=yForecastNew,fitted=yFitted, lower=yLowerNew,upper=yUpperNew,
level=level,legend=!silentLegend,main=modelname,cumulative=cumulative);
}
else{
graphmaker(actuals=y,forecast=yForecastNew,fitted=yFitted,
legend=!silentLegend,main=modelname,cumulative=cumulative);
}
}
return(bestModel);
}
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.