update to allow for specific values of first-order and seasonal differencing outside of forecast::auto.arima

Author: elray1

Diff for: NAMESPACE

@@ -1,11 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(simulate,sarimaTD)
+export(do_difference)
 export(do_initial_transform)
-export(do_seasonal_difference)
 export(fit_sarima)
 export(interpolate_and_clean_missing)
 export(invert_bc_transform)
+export(invert_difference)
 export(invert_initial_transform)
-export(invert_seasonal_difference)
 export(sarimaTD_FF)

Diff for: R/estimation.R

@@ -13,8 +13,10 @@
 #'   in case of observations of 0, and also ensures that the de-transformed
 #'   values will always be at least -0.5, so that they round up to non-negative
 #'   values.
-#' @param seasonal_difference boolean; take a seasonal difference before passing
-#'   to auto.arima?
+#' @param sarimaTD_d integer order of first differencing done before passing to
+#'   auto.arima
+#' @param sarimaTD_D integer order of seasonal differencing done before passing
+#'   to auto.arima
 #' @param d order of first differencing argument to auto.arima.
 #' @param D order of seasonal differencing argument to auto.arima.
 #' @param ... arguments passed on to forecast::auto.arima
@@ -37,7 +39,8 @@ fit_sarima <- function(
   ts_frequency,
   transformation = "box-cox",
   bc_gamma = 0.5,
-  seasonal_difference = TRUE,
+  sarimaTD_d = 0,
+  sarimaTD_D = 1,
   d = NA,
   D = NA,
   ...) {
@@ -66,15 +69,10 @@ fit_sarima <- function(
     transformation = transformation,
     bc_params = est_bc_params)
 
-  ## Initial seasonal differencing, if necessary
-  if(seasonal_difference) {
-    differenced_y <- do_seasonal_difference(
-      y = transformed_y,
-      ts_frequency = ts_frequency)
-  } else {
-    differenced_y <- ts(transformed_y, frequency = ts_frequency)
-  }
-
+  ## Initial differencing, if necessary
+  differenced_y <- do_difference(transformed_y, d = sarimaTD_d, D = sarimaTD_D,
+    frequency = ts_frequency)
+  
   ## Get SARIMA fit
   if(identical(transformation, "forecast-box-cox")) {
     ## box-cox transformation done by auto.arima
@@ -98,7 +96,8 @@ fit_sarima <- function(
   }
 
   sarima_fit$sarimaTD_call <- match.call()
-  for(param_name in c("y", "ts_frequency", "transformation", "seasonal_difference", "d", "D")) {
+  for(param_name in c("y", "ts_frequency", "transformation", "sarimaTD_d",
+                      "sarimaTD_D", "d", "D")) {
     sarima_fit[[paste0("sarimaTD_used_", param_name)]] <- get(param_name)
   }
   if(identical(transformation, "box-cox")) {

Diff for: R/prediction.R

@@ -3,8 +3,8 @@
 #' Simulate predictive trajectories from an ARIMA model
 #'
 #' This is directly taken from the forecast.Arima function from the forecast
-#' package, but it returns the simulated trajectories which were not returned in
-#' the original function definition.  This function calls
+#' package, but it returns the simulated trajectories which were not returned
+#' in the original function definition.  This function calls
 #' forecast:::simulate.Arima, which is NOT exported from the forecast package!!
 #'
 #' @param object an Arima fit object (with class "Arima")
@@ -63,12 +63,12 @@ simulate.sarimaTD <- function(
 
   if(is.null(object$sarimaTD_call)) {
     transformation <- "none"
-    seasonal_difference <- FALSE
+    sarimaTD_d <- sarimaTD_D <- 0
     ts_frequency <- object$arma[5]
   } else {
     transformation <- object$sarimaTD_used_transformation
-    seasonal_difference <-
-      object$sarimaTD_used_seasonal_difference
+    sarimaTD_d <- object$sarimaTD_used_sarimaTD_d
+    sarimaTD_D <- object$sarimaTD_used_sarimaTD_D
     ts_frequency <- object$arma[5]
   }
 
@@ -85,14 +85,9 @@ simulate.sarimaTD <- function(
     bc_params = est_bc_params)
 
   ## Initial seasonal differencing, if necessary
-  if(seasonal_difference) {
-    differenced_y <- do_seasonal_difference(
-      y = transformed_y,
-      ts_frequency = ts_frequency)
-  } else {
-    differenced_y <- ts(transformed_y, frequency = ts_frequency)
-  }
-
+  differenced_y <- do_difference(transformed_y, d = sarimaTD_d, D = sarimaTD_D,
+    frequency = ts_frequency)
+  
   ## Drop leading missing values, fill in internal missing values via linear
   ## interpolation.  This is necessary to ensure non-missing predictions if the
   ## sarima model has a moving average component.
@@ -120,27 +115,20 @@ simulate.sarimaTD <- function(
   ## Get to trajectories for originally observed time series ("orig") by
   ## adding seasonal lag of incidence and inverting the transformation
   orig_trajectory_samples <- raw_trajectory_samples
-  if(seasonal_difference) {
-    for(i in seq_len(nsim)) {
-      orig_trajectory_samples[i, ] <-
-        invert_seasonal_difference(
-          dy = raw_trajectory_samples[i, ],
-          y = transformed_y,
-          ts_frequency = ts_frequency)
-      orig_trajectory_samples[i, ] <-
-        invert_initial_transform(
-          y = orig_trajectory_samples[i, ],
-          transformation = transformation,
-          bc_params = est_bc_params)
-    }
-  } else {
-    for(i in seq_len(nsim)) {
-      orig_trajectory_samples[i, ] <-
-        invert_initial_transform(
-          y = raw_trajectory_samples[i, ],
-          transformation = transformation,
-          bc_params = est_bc_params)
-    }
+  for(i in seq_len(nsim)) {
+    orig_trajectory_samples[i, ] <-
+      invert_difference(
+        dy = raw_trajectory_samples[i, ],
+        y = transformed_y,
+        d = sarimaTD_d,
+        D = sarimaTD_D,
+        frequency = ts_frequency)
+    
+    orig_trajectory_samples[i, ] <-
+      invert_initial_transform(
+        y = orig_trajectory_samples[i, ],
+        transformation = transformation,
+        bc_params = est_bc_params)
   }
 
   attr(orig_trajectory_samples, "seed") <- seed

Diff for: R/utils.R

@@ -89,45 +89,88 @@ invert_bc_transform <- function(b, lambda, gamma) {
   return(z - gamma)
 }
 
-#' Do first-order seasonal differencing (go from original time series values to
-#' seasonally differenced time series values).
+#' Do first-order and seasonal differencing (go from original time series
+#' to differenced time series).
 #'
 #' @param y a univariate time series or numeric vector.
-#' @param ts_frequency frequency of time series.  Must be provided if y is not
-#'   of class "ts".  See the help for stats::ts for more.
+#' @param d order of first differencing
+#' @param D order of seasonal differencing
+#' @param frequency frequency of time series.  Must be provided if y is not
+#'   of class "ts" and D > 0.  See the help for stats::ts for more.
 #'
-#' @return a seasonally differenced time series object (of class 'ts'),
+#' @return a differenced time series object (of class 'ts'),
 #'   padded with leading NAs.
 #'
 #' @export
-do_seasonal_difference <- function(y, ts_frequency) {
-  differenced_y <- ts(c(rep(NA, ts_frequency),
-    y[seq(from = ts_frequency + 1, to = length(y))] -
-      y[seq(from = 1, to = length(y) - ts_frequency)]),
-    frequency = ts_frequency)
-  return(differenced_y)
+do_difference <- function(y, d = 0, D = 0, frequency = 1) {
+  # first differencing
+  for(i in seq_len(d)) {
+    y <- ts(
+      c(NA,
+        y[seq(from = 1 + 1, to = length(y))] -
+          y[seq(from = 1, to = length(y) - 1)]),
+    frequency = frequency)
+  }
+
+  # seasonal differencing
+  if(D > 0 && frequency < 2) {
+    stop("It doesn't make sense to do seasonal differencing with a time series frequency of 1.")
+  }
+  for(i in seq_len(D)) {
+    y <- ts(
+      c(rep(NA, frequency),
+        y[seq(from = frequency + 1, to = length(y))] -
+          y[seq(from = 1, to = length(y) - frequency)]),
+      frequency = frequency)
+  }
+  
+  return(y)
 }
 
-#' Invert first-order seasonal differencing (go from seasonally differenced time
-#' series values to original time series values).
+#' Invert first-order and seasonal differencing (go from seasonally differenced
+#' time series to original time series).
 #'
-#' @param dy a first-order seasonally differenced univariate time series with
-#'   values like y_{t} - y_{t - ts_frequency}
+#' @param dy a first-order and/or seasonally differenced univariate time series
+#'   with values like y_{t} - y_{t - ts_frequency}
 #' @param y a univariate time series or numeric vector with values like
 #'   y_{t - ts_frequency}.
-#' @param ts_frequency frequency of time series.  Must be provided if y is not
-#'   of class "ts".  See the help for stats::ts for more.
+#' @param d order of first differencing
+#' @param D order of seasonal differencing
+#' @param frequency frequency of time series.  Must be provided if y is not
+#'   of class "ts" and D > 0.  See the help for stats::ts for more.
 #'
 #' @details y may have longer length than dy.  It is assumed that dy "starts"
-#'   one time index after y "ends": that is, if y is of length T then
-#'   dy[1] = y[T + 1] - y[T + 1 - ts_frequency]
+#'   one time index after y "ends": that is, if y is of length T, d = 0, and
+#'   D = 1 then dy[1] = y[T + 1] - y[T + 1 - ts_frequency]
 #'
 #' @return a time series object (of class 'ts')
 #'
 #' @export
-invert_seasonal_difference <- function(dy, y, ts_frequency) {
-  return(ts(dy + y[length(y) + seq_along(dy) - ts_frequency],
-    freq = ts_frequency))
+invert_difference <- function(dy, y, d, D, frequency) {
+  for(i in seq_len(d)) {
+    y_dm1 <- do_difference(y, d = d-i, D = D, frequency = frequency)
+    dy_full <- c(y_dm1, dy)
+    for(t in seq_len(length(dy))) {
+      dy_full[length(y_dm1) + t] <- dy_full[length(y_dm1) + t - 1] + dy_full[length(y_dm1) + t]
+    }
+    dy <- dy_full[length(y_dm1) + seq_along(dy)]
+  }
+
+  for(i in seq_len(D)) {
+    y_dm1 <- do_difference(y, d = 0, D = D-i, frequency = frequency)
+    dy_full <- c(y_dm1, dy)
+    for(t in seq_len(length(dy))) {
+      dy_full[length(y_dm1) + t] <- dy_full[length(y_dm1) + t - frequency] + dy_full[length(y_dm1) + t]
+    }
+    dy <- dy_full[length(y_dm1) + seq_along(dy)]
+  }
+  
+#  for(i in seq_len(D)) {
+#    y_dm1 <- do_difference(y, d = 0, D = D-i, frequency = frequency)
+#    dy <- dy + y_dm1[length(y_dm1) + seq_along(dy) - frequency]
+#  }
+  
+  return(ts(dy, frequency = frequency))
 }
 
 #' Remove leading values that are infinite or missing, and replace all internal