Nothing
#' @aliases auto.msarima
#'
#' @examples
#'
#' # The best ARIMA for the data
#' ourModel <- auto.msarima(x,orders=list(ar=c(2,1),i=c(1,1),ma=c(2,1)),lags=c(1,12),
#' h=18,holdout=TRUE)
#'
#' # The other one using optimised states
#' \donttest{auto.msarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
#' h=18,holdout=TRUE)}
#'
#' # And now combined ARIMA
#' \donttest{auto.msarima(x,orders=list(ar=c(3,2),i=c(2,1),ma=c(3,2)),lags=c(1,12),
#' combine=TRUE,h=18,holdout=TRUE)}
#'
#' plot(forecast(ourModel, h=18, interval="simulated"))
#'
#' @rdname msarima
#' @export
auto.msarima <- function(y, orders=list(ar=c(3,3),i=c(2,1),ma=c(3,3)), lags=c(1,frequency(y)),
initial=c("optimal","backcasting","complete"), ic=c("AICc","AIC","BIC","BICc"),
loss=c("likelihood","MSE","MAE","HAM","MSEh","TMSE","GTMSE","MSCE"),
h=10, holdout=FALSE,
bounds=c("usual","admissible","none"),
silent=TRUE,
xreg=NULL, regressors=c("use","select","adapt"), initialX=NULL, ...){
# This is a wrapper function for adam with order selection
if(!is.null(xreg) && is.numeric(y)){
data <- cbind(y=as.data.frame(y),as.data.frame(xreg));
data <- as.matrix(data)
data <- ts(data, start=start(y), frequency=frequency(y));
colnames(data)[1] <- "y";
# Give name to the explanatory variables if they do not have them
if(is.null(names(xreg))){
if(!is.null(ncol(xreg))){
colnames(data)[-1] <- paste0("x",c(1:ncol(xreg)));
}
else{
colnames(data)[-1] <- "x";
}
}
}
else{
data <- y;
}
# Fix orders
orders$select <- TRUE;
return(adam(data, model="NNN", orders=orders, lags=lags, distribution="dnorm",
initial=initial, loss=loss, h=h, holdout=holdout, bounds=bounds,
silent=silent, regressors=regressors, ...))
# Function estimates several msarima models and selects the best one using the selected information criterion.
#
# Copyright (C) 2015 - 2016 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;
# }
# }
#
# interval <- "none";
# cumulative <- FALSE;
#
# ##### Set environment for ssInput and make all the checks #####
# environment(ssAutoInput) <- environment();
# ssAutoInput("auto.msarima",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 the 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 <- msarima(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,"%"));
# }
# # Originally, we only have a constant
# nParamOriginal <- 1;
# if(silent[1]=="d"){
# cat("I: ");cat(iOrders[d,]);cat(", ");
# }
# testModel <- msarima(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;
# 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;
#
# if(silent[1]=="d"){
# cat("MA: ");cat(maTest);cat(", ");
# }
# testModel <- msarima(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);
# # Exclude the variance from the number of parameters
# nParamMA <- nparam(testModel)-1;
# nParamNew <- nParamOriginal + nParamMA;
# ICValue <- icCorrector(testModel$ICs, 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;
#
# if(silent[1]=="d"){
# cat("AR: ");cat(arTest);cat(", ");
# }
# testModel <- msarima(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);
# # Exclude the variance from the number of parameters
# nParamAR <- nparam(testModel)-1;
# nParamNew <- nParamOriginal + nParamMA + nParamAR;
# ICValue <- icCorrector(testModel$ICs, 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 <- msarima(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, 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 <- msarima(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;
# 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("No constant: "); cat(ICValue); cat("\n");
# }
# if(ICValue < bestIC){
# bestModel <- testModel;
# constantValue <- FALSE;
# bestIC <- ICValue;
# }
# else{
# constantValue <- TRUE;
# }
# }
# }
# if(silent[1]=="d"){
# cat("Best IC: "); cat(bestIC); cat("\n");
# }
#
# 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=c("smooth","msarima"));
# }
# else{
# #### Reestimate the best model in order to get rid of bias ####
# bestModel <- msarima(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=FALSE);
# # }
# }
#
# 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.