update models

.Rbuildignore

@@ -19,3 +19,4 @@ vignettes/figs/
 ^CRAN-SUBMISSION$
 ^sticker$
 ^development_packages\.R$
+^fit_models\.R$

DESCRIPTION

@@ -39,6 +39,7 @@ Imports:
     MASS,
     oeli (>= 0.3.0),
     padr,
+    pracma,
     progress,
     Rcpp,
     stats,

NAMESPACE

@@ -42,6 +42,7 @@ export(reorder_states)
 export(set_controls)
 export(simulate_hmm)
 export(viterbi)
+importFrom(MASS,ginv)
 importFrom(Rcpp,evalCpp)
 importFrom(checkmate,assert_integerish)
 importFrom(checkmate,assert_number)
@@ -78,6 +79,7 @@ importFrom(graphics,points)
 importFrom(graphics,text)
 importFrom(graphics,title)
 importFrom(padr,pad)
+importFrom(pracma,hessian)
 importFrom(stats,AIC)
 importFrom(stats,BIC)
 importFrom(stats,acf)

R/compute_ci.R

@@ -32,21 +32,21 @@ compute_ci <- function(x, alpha = 0.05) {
   }
 
   ### compute confidence intervals using the inverse Hessian approach
-  inv_fisher <- 1 / x$hessian_diagonal
-  sds <- suppressWarnings(sqrt(inv_fisher))
+  inverse_fisher <- x$inverse_fisher
+  sds <- suppressWarnings(sqrt(inverse_fisher))
   z_alpha <- stats::qnorm(p = 1 - alpha / 2)
   lower_limit <- x$estimate - z_alpha * sds
   upper_limit <- x$estimate + z_alpha * sds
 
   ### if negative variance, replace by NA_real_
-  bad_inv_fisher <- which(
+  bad_inverse_fisher <- which(
     !vapply(
-      inv_fisher, checkmate::test_number, logical(1), na.ok = FALSE, 
+      inverse_fisher, checkmate::test_number, logical(1), na.ok = FALSE, 
       finite = TRUE, lower = 0
     )
   )
-  lower_limit[bad_inv_fisher] <- NA_real_
-  upper_limit[bad_inv_fisher] <- NA_real_
+  lower_limit[bad_inverse_fisher] <- NA_real_
+  upper_limit[bad_inverse_fisher] <- NA_real_
 
   ### create and return output
   out <- lapply(

R/data_and_models.R

@@ -155,7 +155,7 @@
 #' fit    = list("runs" = 10, "gradtol" = 1e-6, "steptol" = 1e-6)
 #' )
 #' dax_data <- prepare_data(controls)
-#' dax_model_2n <- fit_model(dax_data)
+#' dax_model_2n <- fit_model(dax_data, seed = 1)
 #' dax_model_2n <- decode_states(dax_model_2n)
 #' dax_model_2n <- compute_residuals(dax_model_2n)
 #' summary(dax_model_2n)
@@ -178,22 +178,26 @@
 #' @details
 #' The model was estimated via:
 #' \preformatted{
-#' # TODO refit
 #' controls <- set_controls(
-#' states = 3,
-#' sdds   = "t",
-#' data   = list(
-#'   file        = dax,
-#'   date_column = "Date",
-#'   data_column = "Close",
-#'   logreturns  = TRUE,
-#'   from        = "2000-01-03",
-#'   to          = "2022-12-31"
-#' ),
-#' fit    = list("runs" = 100, "gradtol" = 1e-6, "steptol" = 1e-6)
+#'   states = 3,
+#'   sdds   = "t",
+#'   data   = list(
+#'     file        = dax,
+#'     date_column = "Date",
+#'     data_column = "Close",
+#'     logreturns  = TRUE,
+#'     from        = "2000-01-03",
+#'     to          = "2022-12-31"
+#'   ),
+#'   fit    = list(
+#'     runs        = 100, 
+#'     iterlim     = 300,
+#'     gradtol     = 1e-6,
+#'     steptol     = 1e-6
+#'   )
 #' )
 #' dax_data <- prepare_data(controls)
-#' dax_model_3t <- fit_model(dax_data)
+#' dax_model_3t <- fit_model(dax_data, seed = 1, ncluster = 10)
 #' dax_model_3t <- decode_states(dax_model_3t)
 #' dax_model_3t <- compute_residuals(dax_model_3t)
 #' summary(dax_model_3t)
@@ -217,7 +221,7 @@
 #' @details
 #' The model was estimated via:
 #' \preformatted{
-#' # TODO refit
+#' try({
 #' controls <- set_controls(
 #'   hierarchy = TRUE,
 #'   states    = c(2, 2),
@@ -230,11 +234,18 @@
 #'     logreturns = c(TRUE, TRUE)
 #'   ),
 #'   fit       = list(
-#'     runs = 200
+#'     runs       = 200, 
+#'     iterlim    = 300,
+#'     gradtol    = 1e-6,
+#'     steptol    = 1e-6
 #'   )
 #' )
 #' dax_vw_data <- prepare_data(controls)
-#' dax_vw_model <- fit_model(dax_vw_data)
+#' dax_vw_model <- fit_model(dax_vw_data, seed = 1, ncluster = 10)
+#' dax_vw_model <- decode_states(dax_vw_model)
+#' dax_vw_model <- compute_residuals(dax_vw_model)
+#' summary(dax_vw_model)
+#' })
 #' }
 #'
 #' @format An object of class \code{\link{fHMM_model}}.
@@ -254,7 +265,6 @@
 #' @details
 #' The model was estimated via:
 #' \preformatted{
-#' # TODO refit
 #' controls <- set_controls(
 #'   hierarchy = TRUE,
 #'   states    = c(3, 2),
@@ -274,7 +284,7 @@
 #'   )
 #' )
 #' unemp_spx_data <- prepare_data(controls)
-#' unemp_spx_model_3_2 <- fit_model(unemp_spx_data)
+#' unemp_spx_model_3_2 <- fit_model(unemp_spx_data, seed = 1)
 #' }
 #'
 #' @format An object of class \code{\link{fHMM_model}}.
@@ -294,10 +304,10 @@
 #' The model was estimated via:
 #' \preformatted{
 #' controls <- set_controls(
-#' states  = 2,
-#' sdds    = "gamma(mu = 1|2)",
-#' horizon = 200,
-#' fit     = list("runs" = 10)
+#'   states  = 2,
+#'   sdds    = "gamma(mu = 1|2)",
+#'   horizon = 200,
+#'   runs    = 10
 #' )
 #' pars <- fHMM_parameters(
 #'   controls = controls,
@@ -317,29 +327,3 @@
 #' @keywords model
 "sim_model_2gamma"
 
-#' Simulated 4-state HMM with log-normal distributions
-#'
-#' @description
-#' A pre-computed 4-state HMM with state-dependent log-normal distributions
-#' on \code{1000} simulated observations.
-#'
-#' @usage data("sim_model_4lnorm")
-#'
-#' @details
-#' The model was estimated via:
-#' \preformatted{
-#' # TODO refit
-#' controls <- set_controls(
-#'   states  = 4,
-#'   sdds    = "lognormal",
-#'   horizon = 1000,
-#'   fit     = list(runs = 50)
-#' )
-#' data_sim <- prepare_data(controls, seed = 1)
-#' sim_model_4lnorm <- fit_model(data_sim, seed = 1)
-#' }
-#'
-#' @format An object of class \code{\link{fHMM_model}}.
-#'
-#' @keywords model
-"sim_model_4lnorm"

R/fHMM-package.R

@@ -48,7 +48,9 @@
 #' @importFrom grDevices col2rgb
 #' @importFrom grDevices colorRampPalette
 #' @importFrom grDevices rgb
+#' @importFrom MASS ginv
 #' @importFrom padr pad
+#' @importFrom pracma hessian
 #' @importFrom Rcpp evalCpp
 #' @importFrom stats acf
 #' @importFrom stats AIC

R/fHMM_model.R

@@ -22,9 +22,8 @@
 #' A \code{numeric} vector, the model log-likelihoods in all optimization runs.
 #' @param gradient
 #' A \code{numeric} vector, the gradient at the optimum.
-#' @param hessian_diagonal
-#' A \code{numeric} vector, the diagonal elements of the approximated Hessian at 
-#' the optimum.
+#' @param inverse_fisher
+#' A \code{numeric} vector, the inverse Fisher information for each parameter.
 #' @param decoding
 #' A \code{numeric} vector, the decoded time series.
 #' @param alpha
@@ -41,7 +40,7 @@
 
 fHMM_model <- function(
     data, estimate, nlm_output, estimation_time, ll, lls, gradient, 
-    hessian_diagonal, decoding
+    inverse_fisher, decoding
 ) {
   structure(
     list(
@@ -52,8 +51,9 @@ fHMM_model <- function(
       "ll" = ll,
       "lls" = lls,
       "gradient" = gradient,
-      "hessian_diagonal" = hessian_diagonal,
-      "decoding" = NULL
+      "inverse_fisher" = inverse_fisher,
+      "decoding" = NULL,
+      "system_information" = oeli::system_information()
     ), 
     class = "fHMM_model"
   )

R/fit_model.R

@@ -200,14 +200,25 @@ fit_model <- function(
     )
   }
 
-  ### compute diagonal elements of Hessian
+  ### compute inverse Fisher information
   if (verbose) message("Approximating Hessian...")
-  hessian_diagonal <- suppressWarnings(pracma::hessdiag(
+  fisher <- pracma::hessdiag(
     f = target,
     x = mods[[which.max(lls)]][["estimate"]],
     observations = data[["data"]],
     controls = data[["controls"]]
-  ))
+  )
+  if (all(fisher > 0)) {
+    inverse_fisher <- 1 / fisher
+  } else {
+    hessian <- suppressWarnings(pracma::hessian(
+      f = target,
+      x0 = mods[[which.max(lls)]][["estimate"]],
+      observations = data[["data"]],
+      controls = data[["controls"]]
+    ))
+    inverse_fisher <- diag(MASS::ginv(hessian))
+  }
 
   ### final message
   if (verbose) message("Fitting completed!")
@@ -228,7 +239,7 @@ fit_model <- function(
     ll = ll,
     lls = lls,
     gradient = mod$gradient,
-    hessian_diagonal = hessian_diagonal,
+    inverse_fisher = inverse_fisher,
     decoding = NULL
   )
 

R/reorder_states.R

@@ -115,7 +115,7 @@ reorder_states <- function(x, state_order = "mean") {
   permut_all <- diag(length(x$estimate))[match, ]
   x$estimate <- parUncon
   x$gradient <- as.vector(permut_all %*% x$gradient)
-  x$hessian_diagonal <- as.vector(permut_all %*% x$hessian_diagonal)
+  x$inverse_fisher <- as.vector(permut_all %*% x$inverse_fisher)
   x$nlm_output$estimate <- x$estimate
   x$nlm_output$gradient <- x$gradient