Merge pull request #440 from tidyverts/439-update-ets-to-allow-missing-values

Updated ETS to allow missing values

Author: mitchelloharawild

DESCRIPTION

@@ -68,4 +68,4 @@ Encoding: UTF-8
 Language: en-GB
 Roxygen: list(markdown = TRUE, roclets=c('rd', 'collate',
     'namespace'))
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3

NAMESPACE

@@ -54,6 +54,7 @@ S3method(glance,model_mean)
 S3method(hfitted,ARIMA)
 S3method(hfitted,ETS)
 S3method(interpolate,ARIMA)
+S3method(interpolate,ETS)
 S3method(interpolate,TSLM)
 S3method(interpolate,model_mean)
 S3method(model_sum,AR)

NEWS.md

@@ -2,6 +2,7 @@
 
 ## Improvements
 
+* `ETS()` now supports missing values.
 * Documentation improvements.
 
 ## Bug fixes

R/ets.R

@@ -17,10 +17,6 @@ train_ets <- function(.data, specials, opt_crit,
   y <- unclass(.data)[[measured_vars(.data)]]
   idx <- unclass(.data)[[index_var(.data)]]
 
-  if (any(is.na(y))) {
-    abort("ETS does not support missing values.")
-  }
-
   # Build possible models
   model_opts <- expand.grid(
     errortype = ets_spec$error$method,
@@ -33,7 +29,7 @@ train_ets <- function(.data, specials, opt_crit,
 
   # Remove bad models
   if (NROW(model_opts) > 1) {
-    if (min(y) <= 0) {
+    if (min(y, na.rm = TRUE) <= 0) {
       model_opts <- model_opts[model_opts$errortype != "M", ]
     }
     if (restrict) {
@@ -720,4 +716,26 @@ initial_ets_states <- function(object) {
   )
   colnames(states_init) <- unsplit(states_names, states_type)
   states_init
-}
+}
+
+
+#' @inherit interpolate.model_mean
+#'
+#' @examples
+#' library(tsibbledata)
+#'
+#' olympic_running |>
+#'   model(mean = ETS(Time)) %>%
+#'   interpolate(olympic_running)
+#' @export
+interpolate.ETS <- function(object, new_data, specials, ...) {
+  # Get missing values
+  y <- unclass(new_data)[[measured_vars(new_data)]]
+  miss_val <- which(is.na(y))
+  # Forward fitted values
+  forward_fits <- object$est[[".fitted"]][miss_val]
+  # Ideally, we would also apply the model to the reversed time series
+  # and get the backward fitted values, then combine the two.
+  new_data[[measured_vars(new_data)]][miss_val] <- forward_fits
+  new_data
+}

R/etsmodel.R

@@ -140,16 +140,11 @@ estimate_ets <- function(y, m, init.state, errortype, trendtype, seasontype,
     alpha = unname(alpha), beta = unname(beta),
     gamma = unname(gamma), phi = unname(phi), init.state
   )
-  if (errortype == "A") {
-    fits <- y - e$e
-  } else {
-    fits <- y / (1 + e$e)
-  }
 
   return(list(
     loglik = -0.5 * e$lik, aic = aic, bic = bic, aicc = aicc,
     mse = mse, amse = amse, mae = mae,
-    residuals = e$e, fitted = fits,
+    residuals = e$e, fitted = e$fits,
     states = states, par = fit.par
   ))
 }
@@ -429,7 +424,7 @@ pegelsresid.C <- function(y, m, init.state, errortype, trendtype, seasontype, da
   p <- length(init.state)
   x <- numeric(p * (n + 1))
   x[1:p] <- init.state
-  e <- numeric(n)
+  e <- fits <- numeric(n)
   lik <- 0
   if (!damped) {
     phi <- 1
@@ -457,18 +452,15 @@ pegelsresid.C <- function(y, m, init.state, errortype, trendtype, seasontype, da
     as.double(gamma),
     as.double(phi),
     as.double(e),
+    as.double(fits),
     as.double(lik),
     as.double(amse),
     as.integer(nmse),
+    NAOK = TRUE,
     PACKAGE = "fable"
   )
-  if (!is.na(Cout[[13]])) {
-    if (abs(Cout[[13]] + 99999) < 1e-7) {
-      Cout[[13]] <- NA
-    }
-  }
 
-  return(list(lik = Cout[[13]], amse = Cout[[14]], e = Cout[[12]], states = matrix(Cout[[3]], nrow = n + 1, ncol = p, byrow = TRUE)))
+  return(list(lik = Cout[[14]], amse = Cout[[15]], e = Cout[[12]], fits = Cout[[13]], states = matrix(Cout[[3]], nrow = n + 1, ncol = p, byrow = TRUE)))
 }
 
 admissible <- function(alpha, beta, gamma, phi, m) {

man/interpolate.ETS.Rd

@@ -69,6 +69,7 @@ void EtsTargetFunction::init(std::vector<double> & p_y, int p_nstate, int p_erro
 	//	for(int i=0; i < n; i++) this->e.push_back(0);
 	this->amse.resize(30, 0);
 	this->e.resize(n, 0);
+	this->fits.resize(n, 0);
 
 }
 
@@ -165,8 +166,7 @@ void EtsTargetFunction::eval(const double* p_par, int p_par_length) {
 	for(int i=0; i <= p*this->y.size(); i++) state.push_back(0);
 
 	etscalc(&this->y[0], &this->n, &this->state[0], &this->m, &this->errortype, &this->trendtype, &this->seasontype,
-			&this->alpha, &this->beta, &this->gamma, &this->phi, &this->e[0], &this->lik, &this->amse[0], &this->nmse);
-
+			&this->alpha, &this->beta, &this->gamma, &this->phi, &this->e[0], &this->fits[0], &this->lik, &this->amse[0], &this->nmse);
 
 	// Avoid perfect fits
 	if (this->lik < -1e10) this->lik = -1e10;

src/etsTargetFunction.cpp

@@ -5,12 +5,9 @@
 extern "C" {
 
 void etscalc(double *, int *, double *, int *, int *, int *, int *,
-		double *, double *, double *, double *, double *, double *, double *, int *);
+		double *, double *, double *, double *, double *, double *, double *, double *, int *);
 
-void cpolyroot(double *opr, double *opi, int *degree,
-			double *zeror, double *zeroi, Rboolean *fail);
 }
-
 class EtsTargetFunction {
 
 public:
@@ -53,6 +50,7 @@ class EtsTargetFunction {
 	double alpha, beta, gamma, phi;
 
 	std::vector<double> e;
+	std::vector<double> fits;
 	std::vector<double> amse;
 
 	double lik, objval;