see NEWS

Author: alexiosg

DESCRIPTION

@@ -1,13 +1,13 @@
 Package: tsissm
 Type: Package
 Title: Linear Innovations State Space Unobserved Components Model
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c(person("Alexios", "Galanos", role = c("aut", "cre"), comment = c(ORCID = "0009-0000-9308-0457"), email = "alexios@4dscape.com"))
 Maintainer: Alexios Galanos <alexios@4dscape.com>
 Description: Unobserved components time series model using the linear innovations state space representation (single source of error) with choice of error distributions and option for dynamic variance. Methods for estimation using automatic differentiation, automatic model selection and ensembling, prediction, filtering, simulation and backtesting. Based on the model described in Hyndman et al (2012) <doi:10.1198/jasa.2011.tm09771>.
 Depends: R (>= 4.1.0), Rcpp (>= 0.12.9), tsmethods (>= 1.0.0)
 LinkingTo: Rcpp (>= 0.12.9), TMB, RcppEigen
-Imports: TMB (>= 1.7.20), methods, tsaux, tsdistributions, zoo, xts, copula, flextable, data.table, sandwich, nloptr, RTMB, viridisLite, future, future.apply, progressr
+Imports: TMB (>= 1.7.20), methods, tsaux, tsdistributions, zoo, xts, copula, flextable, data.table, sandwich, nloptr, Rsolnp, RTMB, viridisLite, future, future.apply, progressr
 Suggests: 
     rmarkdown,
     tstests,

NAMESPACE

@@ -49,6 +49,7 @@ export(as_flextable)
 export(bread)
 export(estfun)
 export(hresiduals)
+export(issm_control)
 export(issm_modelspec)
 import(data.table)
 import(methods)
@@ -62,6 +63,7 @@ importFrom(RTMB,eigen)
 importFrom(RTMB,matrix)
 importFrom(Rcpp,evalCpp)
 importFrom(Rcpp,loadModule)
+importFrom(Rsolnp,csolnp)
 importFrom(TMB,MakeADFun)
 importFrom(copula,P2p)
 importFrom(copula,normalCopula)

NEWS.md

@@ -1,3 +1,10 @@
+# tsissm 1.0.2
+
+* The SLSQP solver (from nloptr) sometimes will return NaN for the parameters. Added a stop
+clause else it may crash the process (due to RTMB eigen crashing).
+* Added option for solnp solver (csolnp from Rsolnp)
+* Added a solver control function (issm_control)
+
 # tsissm 1.0.1
 
 * Initial CRAN submission based on complete re-write of model on github which was

R/backtest.R

@@ -20,6 +20,7 @@
 #' @param seed an value specifying if and how the random number generator should 
 #' be initialized (\sQuote{seeded}). Either NULL or an integer that will be used in a 
 #' call to set.seed before simulating the response vectors.
+#' @param solver either \dQuote{nlortr} or \dQuote{solnp}.
 #' @param trace whether to show the progress bar. The user is expected to have
 #' set up appropriate handlers for this using the \dQuote{progressr} package.
 #' @param ... not used.
@@ -56,7 +57,7 @@
 #'
 tsbacktest.tsissm.autospec <- function(object, start = floor(length(object$y)/2), end = length(object$y),
                                    h = 1, estimate_every = 1, rolling = FALSE, weights_scheme = c("AIC","BIC","U"), weights = NULL,
-                                   seed = NULL, trace = FALSE, ...)
+                                   seed = NULL, solver = "nloptr", trace = FALSE, ...)
 {
     if (object$top_n > 1) {
         top_n <- object$top_n
@@ -67,9 +68,10 @@ tsbacktest.tsissm.autospec <- function(object, start = floor(length(object$y)/2)
             weights_scheme <- "U"
         }
         out <- .backtest_ensemble(object, start = start, end = end, h = h, estimate_every = estimate_every, rolling = rolling, weights_scheme = weights_scheme, 
-                                  weights = weights, seed = seed, trace = trace)
+                                  weights = weights, seed = seed, solver = solver, trace = trace)
     } else {
-        out <- .backtest_top(object, start = start, end = end, h = h, estimate_every = estimate_every, rolling = rolling, seed = seed, trace = trace)
+        out <- .backtest_top(object, start = start, end = end, h = h, estimate_every = estimate_every, rolling = rolling, seed = seed, solver = solver, 
+                             trace = trace)
     }
     return(out)
 }
@@ -78,7 +80,7 @@ tsbacktest.tsissm.autospec <- function(object, start = floor(length(object$y)/2)
 #' @rdname tsbacktest
 #' @export
 tsbacktest.tsissm.spec <- function(object, start = floor(length(object$target$y_orig)/2), end = length(object$target$y_orig),
-                                   h = 1, estimate_every = 1, rolling = FALSE, seed = NULL, trace = FALSE, ...)
+                                   h = 1, estimate_every = 1, rolling = FALSE, seed = NULL, solver = "nloptr", trace = FALSE, ...)
 {
     parameter <- b <- forecast_dates <- NULL
     data <- xts(object$target$y_orig, object$target$index)
@@ -150,7 +152,7 @@ tsbacktest.tsissm.spec <- function(object, start = floor(length(object$target$y_
                                sampling = object$target$sampling, sample_n = object$variance$sample_n, 
                                init_garch =  object$variance$init_variance, garch_order = object$variance$order, 
                                variance = object$variance$type, distribution = object$distribution$type)
-        mod <- try(estimate(spec, scores = FALSE), silent = TRUE)
+        mod <- try(estimate(spec, solver = solver, scores = FALSE), silent = TRUE)
         model_coef <- coef(mod)
         log_lik <- as.numeric(logLik(mod))
         aic <- as.numeric(AIC(mod))
@@ -240,7 +242,8 @@ tsbacktest.tsissm.spec <- function(object, start = floor(length(object$target$y_
 }
 
 .backtest_top <- function(object, start = floor(length(object$y)/2), end = length(object$y), h = 1,
-                          estimate_every = 1, rolling = FALSE, seed = NULL, trace = FALSE, ...)
+                          estimate_every = 1, rolling = FALSE, seed = NULL, solver = "nloptr", 
+                          trace = FALSE, ...)
 {
     parameter <- b <- forecast_dates <- NULL
     data <- object$y
@@ -302,7 +305,7 @@ tsbacktest.tsissm.spec <- function(object, start = floor(length(object$target$y_
                                sampling = object$sampling, sample_n = object$sample_n, 
                                init_garch =  object$init_garch, garch_order = object$garch_order, 
                                variance = object$variance, distribution = object$distribution, top_n = 1)
-        mod <- try(estimate(spec, trace = FALSE), silent = TRUE)
+        mod <- try(estimate(spec, solver = solver, trace = FALSE), silent = TRUE)
         model_coef <- coef(mod)
         log_lik <- as.numeric(logLik(mod))
         aic <- as.numeric(AIC(mod))
@@ -391,7 +394,7 @@ tsbacktest.tsissm.spec <- function(object, start = floor(length(object$target$y_
 
 .backtest_ensemble <- function(object, start = floor(length(object$y)/2), end = length(object$y),
                                h = 1, estimate_every = 1, rolling = FALSE, weights_scheme = c("U","AIC","BIC"), weights = NULL,
-                               seed = NULL, trace = FALSE, ...)
+                               seed = NULL, solver = "nloptr", trace = FALSE, ...)
 {
     if (weights_scheme == "U") {
         wfun <- function(x) {
@@ -467,7 +470,7 @@ tsbacktest.tsissm.spec <- function(object, start = floor(length(object$target$y_
                                sampling = object$sampling, sample_n = object$sample_n, 
                                init_garch =  object$init_garch, garch_order = object$garch_order, 
                                variance = object$variance, distribution = object$distribution, top_n = top_n)
-        mod <- try(estimate(spec, trace = FALSE), silent = TRUE)
+        mod <- try(estimate(spec, solver = solver, trace = FALSE), silent = TRUE)
         w <- wfun(mod)
         L <- index_table[[i]]
         M <- split(L, by = "filter_date")

R/constraint.R

@@ -287,6 +287,7 @@ make_constraint <- function(spec, fun, issmenv) {
         constraint_fun <- NULL
     } else {
         constraint_fun <- function(pars, fun, issmenv) {
+            if (any(is.nan(pars))) stop("\nsolver returned NaN values. Try a different algorithm.")
             con <- unlist(sapply(C, do.call, list(pars, fun, issmenv)))
             jac <- do.call(rbind, lapply(J, function(x) x(pars, fun, issmenv)))
             return(list("constraints" = con, "jacobian" = jac))
@@ -295,6 +296,151 @@ make_constraint <- function(spec, fun, issmenv) {
     return(constraint_fun)
 }
 
+
+make_constraint_solnp <- function(spec, fun, issmenv) {
+    # case 1 only level
+    # case 2 only level and GARCH
+    # case 3 only level and arma
+    # case 4 only level and arma and GARCH
+    # case 5 level +
+    # case 6 level + and GARCH
+    # case 7 level + and arma
+    # case 8 level + and arma and GARCH
+    case_id <- model_case(spec)
+    use_case <- c(0, 0, 0)
+    not_null <- 0
+    if (case_id[2] == 0) {
+        issm_c <- function(pars, fun, issmenv) {
+            NULL
+        }
+        issm_j <- function(pars, fun, issmenv) {
+            NULL
+        }
+    } else {
+        not_null <- not_null + 1
+        issm_c <- function(pars, fun, issmenv) {
+            issmenv$constraint$fn(pars)
+        }
+        issm_j <- function(pars, fun, issmenv) {
+            if (all(abs(issmenv$constraint$report(pars)$tmp - 1) <= 1e-12)) {
+                J <- matrix(1, ncol = length(pars), nrow = 1)
+            } else {
+                J <- issmenv$constraint$gr(pars)
+            }
+            return(J)
+        }
+        use_case[1] <- 1
+    }
+    if (case_id[3] > 0) {
+        arma_order <- spec$arma$order
+        if (arma_order[1] > 1 & arma_order[2] > 1) {
+            not_null <- not_null + 1
+            arma_c <- function(pars, fun, issmenv) {
+                c(issmenv$arcons$fn(pars),issmenv$macons$fn(pars))
+            }
+            arma_j <- function(pars, fun, issmenv) {
+                rbind(issmenv$arcons$gr(pars), issmenv$macons$gr(pars))
+            }
+        } else if (arma_order[1] > 1 & arma_order[2] <= 1) {
+            not_null <- not_null + 1
+            arma_c <- function(pars, fun, issmenv) {
+                issmenv$arcons$fn(pars)
+            }
+            arma_j <- function(pars, fun, issmenv) {
+                issmenv$arcons$gr(pars)
+            }
+        } else if (arma_order[1] <= 1 & arma_order[2] > 1) {
+            not_null <- not_null + 1
+            arma_c <- function(pars, fun, issmenv) {
+                issmenv$macons$fn(pars)
+            }
+            arma_j <- function(pars, fun, issmenv) {
+                issmenv$macons$gr(pars)
+            }
+        } else {
+            arma_c <- function(pars, fun, issmenv) {
+                NULL
+            }
+            arma_j <- function(pars, fun, issmenv) {
+                NULL
+            }
+        }
+        use_case[2] <- 1
+    } else {
+        arma_c <- function(pars, fun, issmenv) {
+            NULL
+        }
+        arma_j <- function(pars, fun, issmenv) {
+            NULL
+        }
+    }
+    if (case_id[4] > 0) {
+        not_null <- not_null + 1
+        garch_c <- function(pars, fun, issmenv) {
+            garch_constraint(pars, fun, issmenv)
+        }
+        garch_j <- function(pars, fun, issmenv){
+            garch_jacobian(pars, fun, issmenv)
+        }
+        use_case[3] <- 1
+    } else {
+        garch_c <- function(pars, fun, issmenv) {
+            NULL
+        }
+        garch_j <- function(pars, fun, issmenv){
+            NULL
+        }
+    }
+    
+    C <- list()
+    J <- list()
+    k <- 0
+    if (use_case[1] == 1) {
+        k <- k + 1
+        C[[k]] <- issm_c
+        J[[k]] <- issm_j
+    }
+    if (use_case[2] == 1) {
+        k <- k + 1
+        C[[k]] <- arma_c
+        J[[k]] <- arma_j
+    }
+    if (use_case[3] == 1) {
+        k <- k + 1
+        C[[k]] <- garch_c
+        J[[k]] <- garch_j
+    }
+    
+    if (sum(use_case) == 0) {
+        constraint_fun <- NULL
+        constraint_jac <- NULL
+        ineq_lower <- NULL
+        ineq_upper <- NULL
+    } else {
+        if (not_null > 0) {
+            n_c <- length(unlist(sapply(C, do.call, list(fun$par, fun, issmenv))))
+            constraint_fun <- function(pars, fun, issmenv) {
+                if (any(is.nan(pars))) stop("\nsolver returned NaN values. Try a different algorithm.")
+                con <- unlist(sapply(C, do.call, list(pars, fun, issmenv)))
+                return(con)
+            }
+            constraint_jac <- function(pars, fun, issmenv) {
+                if (any(is.nan(pars))) stop("\nsolver returned NaN values. Try a different algorithm.")
+                jac <- do.call(rbind, lapply(J, function(x) x(pars, fun, issmenv)))
+                return(jac)
+            }
+            ineq_lower <- rep(-1e8, n_c)
+            ineq_upper <- rep(0, n_c)
+        } else {
+            constraint_fun <- NULL
+            constraint_jac <- NULL
+            ineq_lower <- NULL
+            ineq_upper <- NULL
+        }
+    }
+    return(list(ineq_fn = constraint_fun, ineq_jac = constraint_jac, ineq_lower = ineq_lower, ineq_upper = ineq_upper))
+}
+
 companion_matrix <- function(params) {
     n <- length(params)
     if (params[n] == 0) stop("The last parameter must not be zero to form a monic polynomial.")