fix bugs in trend_map with VAR1 models

Author: nicholasjclark

R/mvgam.R

@@ -712,6 +712,11 @@ mvgam = function(formula,
     drift <- FALSE
   }
 
+  if(use_lv & trend_model == 'VAR1' & missing(trend_map)){
+    stop('Cannot identify dynamic factor models that evolve as VAR processes',
+         call. = FALSE)
+  }
+
   # JAGS cannot support latent GP or VAR trends
   if(!use_stan & trend_model %in% c('GP', 'VAR1')){
     stop('gaussian process and VAR trends not supported for JAGS',

R/predict.mvgam.R

@@ -76,10 +76,17 @@ predict.mvgam = function(object, newdata,
     }
 
     if(object$trend_model %in% c('VAR1')){
-      family_pars <- list(sigma_obs =
-                            mcmc_chains(object$model_output, 'Sigma')[ ,
-                                                                       seq(1, NCOL(object$ytimes)^2,
-                                                                           by = NCOL(object$ytimes)+1)])
+      if(object$use_lv){
+        family_pars <- list(sigma_obs =
+                              mcmc_chains(object$model_output, 'Sigma')[ ,
+                                                                         seq(1, object$n_lv^2,
+                                                                             by = object$n_lv+1)])
+      } else {
+        family_pars <- list(sigma_obs =
+                              mcmc_chains(object$model_output, 'Sigma')[ ,
+                                                                         seq(1, NCOL(object$ytimes)^2,
+                                                                             by = NCOL(object$ytimes)+1)])
+      }
     }
 
     # Indicators of which trend to use for each observation

R/stan_utils.R

@@ -2320,6 +2320,14 @@ if(trend_model != 'VAR1'){
                     model_file, fixed = TRUE)] <-
       "matrix[n_lv, n_lv] A;"
 
+    model_file[grep("vector[n_series] mu[n - 1];",
+                    model_file, fixed = TRUE)] <-
+      "vector[n_lv] mu[n];"
+
+    model_file[grep("array[n] vector[n_series] mu;",
+                    model_file, fixed = TRUE)] <-
+      "array[n] vector[n_lv] mu;"
+
     model_file[grep("matrix[n_series, n_series] Sigma;",
                     model_file, fixed = TRUE)] <-
       "matrix[n_lv, n_lv] Sigma;"