R/scenario.R

Defines functions scenario

Documented in scenario

#' Title
#' 
#' Description
#' 
#' Details
#' 
#' @param data text
#' @param y text
#' @param X text
#' @param date_y text
#' @param n_ahead text
#' @param variant text
#' @param nnar_model text
#' @param auto_ic text
#' @param accuracy_crit text
#' @param diff_test text
#' @param output text
#' @param write_report text
#' @param dir text
#' @param verbose text
#' @param parallel text
#' 
#' @return text
#' 
#' @author Fagan Ahmedzade
#' 
#' @keywords text
#' 
#' @examples text
#' 
#' @export
scenario <- function(
    data = NULL,
    y = NULL,
    X = NULL,
    date_y = c(1, NROW(data)),
    n_ahead = 3,
    variant = 1,
    nnar_model = FALSE,
    auto_ic = FALSE,
    accuracy_crit = c("RMSE", "MASE", "MAPE", "MAE", "MPE"),
    diff_test = c("KPSS", "ADF", "PP"),
    output = c("best", "consensus", "all_output"),
    write_report = FALSE,
    dir = "Output",
    verbose = FALSE,
    parallel = FALSE
) {
    
    # !!! & and && list (...)
    
    output <- match.arg(output)
    diff_test <- match.arg(diff_test)
    accuracy_crit <- match.arg(accuracy_crit)
    
    start <- date_y[1]
    end <- date_y[2] + n_ahead
    index_y <- c(date_y[1]:date_y[2])
    index_x <- c(1:length(start:end))
    naive_forecast <- NULL
    
    if (n_ahead <= 0 | is.null(n_ahead)) {
        n_ahead <- 5
    }
    
    forecast_param <- names(data)[1]
    if (is.null(forecast_param)) {
        forecast_param <- "Y"
    }
    
    if (is.null(data)) {
        y <- y
        X <- X
    } else {
        data <- as.data.frame(data)
        if (NCOL(data) > 1) {
            y <- ts(data[index_y,  1, ], start = date_y[1], end = date_y)
            X <- data.frame(data[index_x, -1, ], row.names = index_x)
        }
        
        if (NCOL(data) == 1) {
            y <- ts(data[index_y,  1, ], start = date_y[1], end = date_y)
            X <- NULL
        }
    }
    
    if (!is.null(X)) {
        scenarios <- TRUE
        factors <- 1:NCOL(X)
    } else {
        scenarios <- FALSE
    }
    
    if (anyNA(y)) {
        tryCatch(
            y <- forecast::na.interp(y),
            error = function(e) {
                y[is.na(y)] <- mean(y, na.rm = TRUE)
                cat("Message: missing values in", forecast_param,
                    "were not interpolated => interpolated mean.\n"
                )
            }
        )
    }
    
    # naive
    # write_report return
    
    y_d <- set_d(y, test = diff_test)
    y_pq <- set_lag(y, diff = y_d)
    
    if (scenarios) {
        for (i in factors) {
            if (anyNA(X[1:NROW(y), i])) {
                tryCatch(
                    X[[i]][1:NROW(y)] <- forecast::na.interp(X[[i]][1:NROW(y)]),
                    error = function(e) {
                        X[[i]][is.na(X[[i]][1:NROW(y)])] <- mean(X[[i]][1:NROW(y)], na.rm = TRUE)
                        cat(
                            "Message: missing values in", forecast_param,
                            "were not interpolated => interpolated mean.\n"
                        )
                    }
                )
            }
        }
        
        x_d <- vector("integer", NCOL(X))
        x_pq <- vector("list",  NCOL(X))
        x_lag <- vector("list", NCOL(X))
        for (i in factors) {
            x_d[i] <- set_d(X[[i]], test = diff_test)
            x_pq[[i]] <- set_lag(X[[i]], diff = x_d[i])
            x_lag[[i]] <- set_lag(X[[i]], z = y)[3] # 3 - CCF
        }
        
        if (anyNA(X[setdiff(index_x, index_y), ])) {
            X <- na_extrap(X, index_x, index_y, data, x_d, x_pq)
        }
        
        lageed_X <- X
        for (i in factors) {
            lageed_X[, i] <- lagged(X[[i]], k = x_lag[[i]])
        }
        lageed_X <- na.omit(lageed_X)
        names(lageed_X) <- paste0(names(lageed_X), paste0(".lag.", x_lag))
        x_max_lag <- max(sapply(x_lag, function(x) `[`(x)))
        disp_y <- y[(1 + x_max_lag):NROW(y)]
        disp_y <- ts(disp_y, start = start + x_max_lag, end = end - n_ahead)
    }
    
    models <- list(
        func = c(
            if (scenarios) "lmr",
            if (scenarios) "rmlv", 
            if (scenarios) "adl",
            if (scenarios) "arimax", 
            if (auto_ic) "arima_auto", 
            if (scenarios & auto_ic) "arimax_auto", 
            if (scenarios & auto_ic) "nnar", 
            if (nnar_model & auto_ic) "nnar_auto", 
            if (scenarios & nnar_model) "nnar_xreg",
            if (scenarios & auto_ic & nnar_model) "nnar_xreg_auto",
            "arima_sc",
            "ets_sc",
            "lmtrend"
        ),
        args = list(
            y = y, 
            X = if (scenarios) X else NULL, 
            lageed_X = if (scenarios) lageed_X else NULL, 
            disp_y = if (scenarios) disp_y else NULL, 
            y_pq = y_pq, 
            y_d = y_d,
            index_x = index_x, 
            index_y = index_y,
            n_ahead = n_ahead,
            method = c("ML"), 
            ic = c("aic")
        )
    )
    
    models_list <- list()
    for (i in models$func) {
        models_list[[i]] <- NULL
        tryCatch(
            models_list[[i]] <- do.call(
                what = i, 
                args = models$args
            ),
            error = function(e) {
                cat(
                    "Message:", toupper(i), "was not fitted for",
                    paste0("'", forecast_param, "'"),
                    "\b.\n"
                )
            }
        )
    }
    
    summary_models <- list()
    exc_models <- c("nnar", "nnar_xreg", "nnar_auto", "nnar_xreg_auto", "ets_sc")
    suppressWarnings(
        for (i in names(models_list)) {
            if (any(i == exc_models)) {
                summary_models[i] <- models_list[[i]]$model
            } else {
                tryCatch(
                    expr = summary_models[[i]] <- model_stat(
                        fit = models_list[[i]],
                        data = data, y = y, X = X,
                        lag = if (any(i == c("adl", "rmlv", "lmr"))) {
                            TRUE
                        } else {
                            FALSE
                        }
                    ),
                    error = function(e) {
                        cat(
                            "Message: statistics were not calculated for the model",
                            toupper(i), "in the forecast",
                            forecast_param,
                            "\b.\n"
                        )
                    }
                )
            }
        }
    )
    
    if (length(models_list) == 0) {
        naive_forecast <- data.frame(
            row.names = 1:n_ahead,
            forecast = c(
                rep(x = y[max(which(y >= 0))], n_ahead)
            )
        )
        output_stats <- as.vector(naive_forecast)
        cat("naive\n")
        return(output_stats)
    }
    
    ranking <- accur(models_list, accuracy_crit)
    
    forecast_best <- data.frame(
        row.names = start:end,
        y = c(y, models_list[[row.names(ranking)[1]]]$mean)
    )
    
    consensus <- forecast_consensus(
        models_list = models_list, 
        y = y,
        ranking = ranking, 
        accuracy_crit = accuracy_crit,
        start = start, 
        end = end, 
        n_ahead = n_ahead
    )
    
    tables <- forecast_tables(
        models_list = models_list, 
        consensus = consensus,
        end = end, 
        n_ahead = n_ahead, 
        index_x = index_x, 
        index_y = index_y
    )
    
    output_stats <- list(
        models = models_list,
        forecast_best = forecast_best,
        consensus = consensus$forecast,
        summary = summary_models,
        ranking_score = ranking,
        models_weight = t(data.frame(weight = consensus$models_weight * 100)),
        tables = tables
    )
    
    forecast_verbose <- function(data, all_forecasts, dir = '') {
        cat(
            readLines(
                con = paste0(
                    "Output/Reports_variant_1/", 
                    paste0(forecast_param, ".txt")
                )
            )[1:18],
            fill = 18 - 1
        )
        cat("MODELS\n")
        print(tables$all_forecasts)
    }
    
    if (write_report) {
        forecast_report(
            forecast_param = forecast_param, 
            start = start, 
            end = end, 
            n_ahead = n_ahead, 
            X = X, 
            index_x = index_x, 
            index_y = index_y, 
            ranking = ranking,
            accuracy_crit = accuracy_crit, 
            models_list = models_list,
            summary_models = summary_models, 
            variant = variant, 
            output = output,
            models_weight = consensus$models_weight, 
            forecasts_lower = tables$forecasts_lower,
            forecasts_upper = tables$forecasts_upper,
            all_forecasts = tables$all_forecasts
        )
        
        if (write_report & verbose) {
            forecast_verbose(data, all_forecasts)
        }
    }
    
    if (output == "best") {
        return(forecast_best)
    }
    
    if (output == "consensus") {
        return(consensus$forecast)
    }
    
    if (output == "all_output") {
        return(structure(output_stats, class = "scenario"))
    }
}
faganok/scenario documentation built on Nov. 28, 2017, 4:06 p.m.