V0.1.1 release

Author: Resul Akay

CRAN-SUBMISSION

@@ -0,0 +1,3 @@
+Version: 0.1.1
+Date: 2022-10-23 19:55:28 UTC
+SHA: 8f9d65c189202cc31ce9d95940002e2dbfc29d4c

DESCRIPTION

@@ -6,11 +6,10 @@ Authors@R:
            family = "Akay",
            role = c("aut", "cre"),
            email = "resulakay1@gmail.com")
-Description: The caretForecast package makes it easy to create conformal time 
-  series forecasts using the caret infrastructure. The primary function ARml 
-  provides access to state-of-the-art machine learning models for forecasting 
-  applications. The hyperparameter of each model is selected based on time series
-  cross-validation, and forecasting is done recursively.
+Description: Conformal time series forecasting using the caret infrastructure. 
+  It provides access to state-of-the-art machine learning models for forecasting
+  applications. The hyperparameter of each model is selected based on time 
+  series cross-validation, and forecasting is done recursively.
 License: GPL (>= 3)
 URL: https://github.com/Akai01/caretForecast
 BugReports: https://github.com/Akai01/caretForecast/issues

NAMESPACE

@@ -1,37 +1,41 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(forecast,ARml)
+S3method(predict,conformalRegressor)
+S3method(residuals,ARml)
 export("%<>%")
 export("%>%")
 export(ARml)
 export(accuracy)
 export(autolayer)
 export(autoplot)
+export(conformalRegressor)
 export(forecast)
 export(get_var_imp)
 export(split_ts)
 export(suggested_methods)
 importFrom(caret,train)
 importFrom(caret,trainControl)
 importFrom(caret,varImp)
+importFrom(dplyr,bind_cols)
 importFrom(forecast,BoxCox)
 importFrom(forecast,BoxCox.lambda)
 importFrom(forecast,InvBoxCox)
-importFrom(forecast,accuracy)
 importFrom(forecast,autolayer)
 importFrom(forecast,autoplot)
-importFrom(forecast,bld.mbb.bootstrap)
 importFrom(forecast,fourier)
 importFrom(forecast,is.constant)
 importFrom(forecast,na.interp)
+importFrom(generics,accuracy)
+importFrom(generics,forecast)
 importFrom(magrittr,"%<>%")
 importFrom(magrittr,"%>%")
 importFrom(methods,is)
 importFrom(stats,"tsp<-")
 importFrom(stats,frequency)
 importFrom(stats,is.ts)
 importFrom(stats,predict)
-importFrom(stats,quantile)
+importFrom(stats,residuals)
 importFrom(stats,sd)
 importFrom(stats,start)
 importFrom(stats,time)

R/arml.R

@@ -100,13 +100,13 @@ ARml <- function(y,
                  pre_process = NULL,
                  cv = TRUE,
                  cv_horizon = 4,
-                 initial_window = length(y) - max_lag -  cv_horizon * 2,
+                 initial_window = NULL,
                  fixed_window = FALSE,
                  verbose = TRUE,
                  seasonal = TRUE,
                  K =  frequency(y) / 2,
                  tune_grid = NULL,
-                 lambda = "auto",
+                 lambda = NULL,
                  BoxCox_method = c("guerrero", "loglik"),
                  BoxCox_lower = -1,
                  BoxCox_upper = 2,
@@ -118,27 +118,52 @@ ARml <- function(y,
   if ("ts" %notin% class(y)) {
     stop("y must be a univariate time series")
   }
-
+  freq <- stats::frequency(y)
   length_y <- length(y)
 
-  freq <- stats::frequency(y)
+  if (c(length_y - freq - round(freq / 4)) < max_lag) {
+    if(length_y > 3){
+      max_lag <- length_y + 3 - length_y
+    } else {
+      max_lag <- 1
+    }
+    if(max_lag >= length_y - max_lag - 2){
+      max_lag <- 1
+    }
+
+    warning(paste("Input data is too short. setting max_lag =  ", max_lag))
+  }
 
-  if (length_y < freq) {
+  if (length_y < 3) {
     stop("Not enough data to fit a model")
   }
+  constant_data <- is.constant(na.interp(y))
+  if (constant_data) {
+    warning("Constant data, setting max_lag = 1, seasonal = FALSE, lambda = NULL,
+            pre_process = NULL")
+    pre_process <- NULL
+    lambda <- NULL
+    max_lag <- 1
+  }
+
+  if (!is.null(xreg)) {
+    constant_xreg <- any(apply(as.matrix(xreg), 2,
+                               function(x) {
+                                 is.constant(na.interp(x))
+                                 }
+                               )
+                         )
+    if (constant_xreg) {
+      warning("Constant xreg column, setting pre_process=NULL")
+      pre_process <- NULL
+    }
+  }
 
   if (max_lag <= 0) {
     warning("max_lag increased to 1. max_lag must be max_lag >= 1")
     max_lag <- 1
   }
 
-  if (c(length_y - freq - round(freq / 4)) < max_lag) {
-    warning(paste(
-      "Input data is too short. Reducing max_lags to ",
-      round(length_y - freq - round(freq / 4))
-    ))
-    max_lag <- round(length_y - freq - round(freq / 4))
-  }
 
   if (max_lag != round(max_lag)) {
     max_lag <- round(max_lag)
@@ -151,14 +176,6 @@ ARml <- function(y,
     }
   }
 
-  constant_y <- forecast::is.constant(forecast::na.interp(y))
-
-  if (constant_y) {
-    warning("Constant data, setting max_lag = 1, lambda = NULL")
-    lambda <- NULL
-    max_lag = 1
-  }
-
   if (!is.null(xreg))
   {
     ncolxreg <- ncol(xreg)
@@ -185,7 +202,11 @@ ARml <- function(y,
                      start = time(modified_y)[max_lag + 1],
                      frequency = freq)
 
-  if (seasonal == TRUE | freq > 1)
+  if(length_y - max_lag < freq + 1) {
+    seasonal <- FALSE
+  }
+
+  if (seasonal == TRUE & freq > 1)
   {
     if (K == freq / 2) {
       ncolx <- max_lag + K * 2 - 1
@@ -235,6 +256,22 @@ ARml <- function(y,
     }
   }
 
+  initial_window_setted <- FALSE
+  if(is.null(initial_window)){
+    initial_window <- length_y - max_lag -  cv_horizon * 2
+    message("initial_window = NULL. Setting initial_window = ", initial_window)
+    initial_window_setted <- TRUE
+  }
+  if(initial_window<1 | initial_window >= nrow(x)){
+    initial_window <- length_y - max_lag - 1
+    cv_horizon <- 1
+    if(initial_window_setted){
+      warning("Resetting initial_window = ", initial_window, " cv_horizon = 1")
+    } else {
+      warning("Setting initial_window = ", initial_window, " cv_horizon = 1")
+    }
+  }
+
   model <- caret::train(
     x = x,
     y = as.numeric(modified_y_2),
@@ -303,134 +340,3 @@ ARml <- function(y,
   class(output) <- "ARml"
   return(output)
 }
-
-
-#' @title Forecasting an ARml object
-#'
-#' @param object A list class of ARml
-#' @param h forecast horizon
-#' @param xreg Optionally, a numerical vector or matrix of future external
-#' regressors
-#' @param level Confidence level for prediction intervals.
-#' @param PI If TRUE, prediction intervals are produced, otherwise only point
-#' forecasts are calculated. If PI is FALSE, then level, fan, bootstrap and
-#' npaths are all ignored.
-#' @param num_bs Number of bootstrapped versions to generate.
-#' @param ... Other arguments pased to forecast::forecast()
-#' @return A list class of forecast containing the following elemets
-#' * x : The input time series
-#' * method : The name of the forecasting method as a character string
-#' * mean : Point forecasts as a time series
-#' * lower : Lower limits for prediction intervals
-#' * upper : Upper limits for prediction intervals
-#' * level : The confidence values associated with the prediction intervals
-#' * model : A list containing information about the fitted model
-#' * newxreg : A matrix containing regressors
-#' @author Resul Akay
-#'
-#' @examples
-#'
-#'library(caretForecast)
-#'
-#'train_data <- window(AirPassengers, end = c(1959, 12))
-#'
-#'test <- window(AirPassengers, start = c(1960, 1))
-#'
-#'ARml(train_data, caret_method = "lm", max_lag = 12) -> fit
-#'
-#'forecast(fit, h = length(test), level = c(80,95), PI = TRUE) -> fc
-#'
-#'autoplot(fc)+ autolayer(test)
-#'
-#'accuracy(fc, test)
-#'
-#' @export
-forecast.ARml <- function(object,
-                          h = frequency(object$y),
-                          xreg = NULL,
-                          level = c(80, 95),
-                          PI = FALSE,
-                          num_bs = 1000,
-                          ...) {
-  if (!is.null(object$xreg_fit)) {
-    ncolxreg <- ncol(object$xreg_fit)
-  }
-
-  if (is.null(xreg)) {
-    if (!is.null(object$xreg_fit)) {
-      stop("No regressors provided")
-    }
-  }
-
-  if (!is.null(xreg)) {
-    if (is.null(object$xreg_fit)) {
-      stop("No regressors provided to fitted model")
-    }
-
-    if (ncol(xreg) != ncolxreg) {
-      stop("Number of regressors does not match to fitted model")
-    }
-
-    h <- nrow(xreg)
-    newxreg1 <- xreg
-  }
-
-  if (is.null(h)) {
-    h <- ifelse(frequency(object$y) > 1, 2 * frequency(object$y), 10)
-  }
-
-  if (is.null(xreg)) {
-    newxreg1 <- NULL
-  }
-
-  if (PI) {
-    if (is.null(level)) {
-      warning("level was not provided. Prediction intervals will not be returned")
-      PI <- FALSE
-    }
-  }
-
-  lambda <- object$lambda
-  BoxCox_biasadj <- object$BoxCox_biasadj
-  BoxCox_fvar <- object$BoxCox_fvar
-
-  fc_x <- forecast_loop(object = object, xreg = newxreg1, h = h)
-  x <- fc_x$x
-  y <- fc_x$y
-
-  if (!is.null(lambda)) {
-    y <- forecast::InvBoxCox(y,
-                             lambda = lambda,
-                             biasadj = BoxCox_biasadj,
-                             fvar = BoxCox_fvar)
-  }
-
-  if (PI) {
-    bs_pi <- pi(
-      y = object$y,
-      fc = y,
-      num = num_bs,
-      block_size = NULL,
-      level = level
-    )
-    lower <- bs_pi[["lower"]]
-    upper <- bs_pi[["upper"]]
-  } else {
-    lower <- NULL
-    upper <- NULL
-  }
-
-  output <- list(
-    x = object$y,
-    mean = y,
-    lower = lower,
-    upper = upper,
-    fitted = object$fitted,
-    level = level,
-    newxreg = x,
-    method = object$method,
-    model = object$model
-  )
-  class(output) <- c("forecast", "forecastARml")
-  return(output)
-}

R/fit_conformal_reg.R

@@ -0,0 +1,30 @@
+#' Fit a conformal regressor.
+#'
+#' @param residuals Model residuals.
+#' @param sigmas A vector of difficulty estimates
+#' @author Resul Akay
+#' @return A conformalRegressor object
+#'
+#' @references
+#' Boström, H., 2022. crepes: a Python Package for Generating Conformal
+#' Regressors and Predictive Systems. In Conformal and Probabilistic Prediction
+#' and Applications. PMLR, 179.
+#' \url{https://copa-conference.com/papers/COPA2022_paper_11.pdf}
+#'
+#' @export
+conformalRegressor <- function(residuals, sigmas=NULL) {
+  abs_residuals <- abs(residuals)
+  if(is.null(sigmas)){
+    normalized <- FALSE
+    alphas <- rev(sort(abs_residuals))
+  } else {
+    normalized <- TRUE
+    alphas = rev(sort(abs_residuals/sigmas))
+  }
+  out <- list(
+    alphas = alphas,
+    normalized = normalized
+  )
+  out <- structure(out, class = "conformalRegressor")
+  return(out)
+}