diff --git a/NEON_manuscript/ixodes_hyp_tests.R b/NEON_manuscript/ixodes_hyp_tests.R
index 6a56889f..a646cdb2 100644
--- a/NEON_manuscript/ixodes_hyp_tests.R
+++ b/NEON_manuscript/ixodes_hyp_tests.R
@@ -53,7 +53,7 @@ hyp3 = y ~
 
 # Fit each hypothesis
 burnin = 5000
-n_samples = 2000
+n_samples = 1000
 
 fit_null <- fit_mvgam(data_train = all_data$data_train,
                   data_test = all_data$data_test,
diff --git a/NEON_manuscript/neon_utility_functions.R b/NEON_manuscript/neon_utility_functions.R
index 67666118..5e0ea403 100644
--- a/NEON_manuscript/neon_utility_functions.R
+++ b/NEON_manuscript/neon_utility_functions.R
@@ -3,7 +3,7 @@
 calculate_pit = function(out_gam_mod, series, data_test, data_train){
 
   # Pull out predictions and truths for the specific series
-  ends <- seq(0, dim(MCMCvis::MCMCchains(out_gam_mod$jags_output, 'ypred'))[2],
+  ends <- seq(0, dim(MCMCvis::MCMCchains(out_gam_mod$model_output, 'ypred'))[2],
               length.out = NCOL(out_gam_mod$ytimes) + 1)
   starts <- ends + 1
   starts <- c(1, starts[-c(1, (NCOL(out_gam_mod$ytimes)+1))])
@@ -24,7 +24,7 @@ calculate_pit = function(out_gam_mod, series, data_test, data_train){
     dplyr::arrange(time) %>%
     dplyr::pull(y)
 
-  preds <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'ypred')[,starts[series]:ends[series]]
+  preds <- MCMCvis::MCMCchains(out_gam_mod$model_output, 'ypred')[,starts[series]:ends[series]]
   preds <- t(preds[, (last_obs +1):NCOL(preds)])
   if(any(!is.na(truth))){
 
@@ -83,7 +83,7 @@ calculate_drps = function(out_gam_mod, pred_matrix = NULL, series, data_test, da
     scores
   }
 
-  ends <- seq(0, dim(MCMCvis::MCMCchains(out_gam_mod$jags_output, 'ypred'))[2],
+  ends <- seq(0, dim(MCMCvis::MCMCchains(out_gam_mod$model_output, 'ypred'))[2],
               length.out = NCOL(out_gam_mod$ytimes) + 1)
   starts <- ends + 1
   starts <- c(1, starts[-c(1, (NCOL(out_gam_mod$ytimes)+1))])
@@ -105,7 +105,7 @@ calculate_drps = function(out_gam_mod, pred_matrix = NULL, series, data_test, da
     dplyr::pull(y)
 
   if(is.null(pred_matrix)){
-    preds <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'ypred')[,starts[series]:ends[series]]
+    preds <- MCMCvis::MCMCchains(out_gam_mod$model_output, 'ypred')[,starts[series]:ends[series]]
     preds <- t(preds[, (last_obs +1):NCOL(preds)])
   } else {
     preds <- t(pred_matrix[, (last_obs +1):NCOL(pred_matrix)])
@@ -136,7 +136,7 @@ plot_mvgam_season = function(out_gam_mod, series, data_test, data_train,
                           siteID = as.character(unique(data_train$siteID[which(data_train$series ==
                                                                                  levels(data_train$series)[series])])))
   Xp <- predict(out_gam_mod$mgcv_model, newdata = pred_dat, type = 'lpmatrix')
-  betas <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'b')
+  betas <- MCMCvis::MCMCchains(out_gam_mod$model_output, 'b')
   plot((Xp %*% betas[1,]) ~ pred_dat$season, ylim = range((Xp %*% betas[1,]) + 2.5,
                                                                          (Xp %*% betas[1,]) - 2.5),
 
@@ -172,11 +172,11 @@ plot_mvgam_gdd = function(out_gam_mod, series, data_test, data_train,
                           siteID = as.character(unique(data_train$siteID[which(data_train$series ==
                                                                                  levels(data_train$series)[series])])))
   Xp <- predict(out_gam_mod$mgcv_model, newdata = pred_dat, type = 'lpmatrix')
-  betas <- MCMCvis::MCMCchains(out_gam_mod$jags_output, 'b')
+  betas <- MCMCvis::MCMCchains(out_gam_mod$model_output, 'b')
   preds <- matrix(NA, nrow = 1000, ncol = length(pred_dat$cum_gdd))
   for(i in 1:1000){
     preds[i,] <- rnbinom(length(pred_dat$cum_gdd), mu = exp((Xp %*% betas[i,])),
-                         size = MCMCvis::MCMCsummary(out_gam_mod$jags_output, 'r')$mean)
+                         size = MCMCvis::MCMCsummary(out_gam_mod$model_output, 'r')$mean)
   }
   int <- apply(preds,
                2, hpd, 0.95)
@@ -343,7 +343,7 @@ fit_mvgam = function(data_train,
 
   # Condition the model on the observed data
   if(missing(knots)){
-    out_gam_mod <- mvjagam(formula = formula,
+    out_gam_mod <- mvgam(formula = formula,
                            data_train = data_train,
                            data_test = data_test,
                            burnin = burnin,
@@ -355,7 +355,7 @@ fit_mvgam = function(data_train,
                            family = family,
                            phi_prior = phi_prior)
   } else {
-    out_gam_mod <- mvjagam(formula = formula,
+    out_gam_mod <- mvgam(formula = formula,
                            data_train = data_train,
                            data_test = data_test,
                            burnin = burnin,
@@ -373,7 +373,7 @@ fit_mvgam = function(data_train,
   #### If GAM component is LESS supported, we should see evidence in the form of: ####
   # 1. Poorer convergence of smoothing parameter estimates, suggesting the model
   # is more 'mis-specified' and harder to fit
-  rho_summary <- MCMCvis::MCMCsummary(out_gam_mod$jags_output,
+  rho_summary <- MCMCvis::MCMCsummary(out_gam_mod$model_output,
                                        c('rho'), HPD = TRUE)
 
   # 2. Stronger residual correlations, suggesting we are missing some site-level structure
diff --git a/NEON_manuscript/simulation_component.R b/NEON_manuscript/simulation_component.R
index 8148c6bf..db4551e3 100644
--- a/NEON_manuscript/simulation_component.R
+++ b/NEON_manuscript/simulation_component.R
@@ -31,7 +31,7 @@ sim_results <- lapply(seq_len(nrow(run_parameters)), function(x){
                         seasonality = 'hierarchical')
 
   # Fit a multivariate gam with no seasonality (null model)
-  nullmod <- mvjagam(data_train = sim_data$data_train,
+  nullmod <- mvgam(data_train = sim_data$data_train,
                   data_test = sim_data$data_test,
                   formula = y ~ s(series, bs = 're'),
                   family = 'nb',
@@ -80,7 +80,7 @@ sim_results <- lapply(seq_len(nrow(run_parameters)), function(x){
                                                            correlations$mean_correlations))
 
   # Fit a multivariate gam with hierarchical seasonality
-  hier_mod <- mvjagam(data_train = sim_data$data_train,
+  hier_mod <- mvgam(data_train = sim_data$data_train,
                       data_test = sim_data$data_test,
                       formula = y ~ s(series, bs = 're') +
                         s(season, k = 12, m = 2, bs = 'cc') +
@@ -117,7 +117,7 @@ sim_results <- lapply(seq_len(nrow(run_parameters)), function(x){
   Xp <- rbind(predict(mgcv_hier_mod, type = 'lpmatrix'),
               predict(mgcv_hier_mod, newdata = sim_data$data_test, type = 'lpmatrix'))
   vc <- vcov(mgcv_hier_mod)
-  sim <- MASS::mvrnorm(dim(MCMCvis::MCMCchains(hier_mod$jags_output, 'ypred'))[1],
+  sim <- MASS::mvrnorm(dim(MCMCvis::MCMCchains(hier_mod$model_output, 'ypred'))[1],
                        mu = coef(mgcv_hier_mod), Sigma = vc)
   dims_needed <- dim(exp(Xp %*% t(sim)))
   mus <- as.vector(exp(Xp %*% t(sim)))
diff --git a/NEON_manuscript/testing_tweedie.R b/NEON_manuscript/testing_tweedie.R
index 14cf1cd8..d94da003 100644
--- a/NEON_manuscript/testing_tweedie.R
+++ b/NEON_manuscript/testing_tweedie.R
@@ -1,42 +1,100 @@
 # Generate series with independent latent GP trends
 library(mvgam)
-sim_data <- sim_mvgam(T = 70, n_series = 2,
+sim_data <- sim_mvgam(T = 90, n_series = 2,
                       use_lv = F,
                       trend_model = 'GP',
                       trend_rel = 0.6,
                       family = 'poisson',
-                      mu_obs = c(7, 10))
-mod1stan <- mvgam(data_train = sim_data$data_train,
-                    data_test = sim_data$data_test,
-                    formula = y ~ s(season, k = 12, bs = 'cc'),
-                    knots = list(season = c(0.5, 12.5)),
-                    family = 'poisson',
-                    trend_model = 'GP',
-                    chains = 4,
-                    burnin = 1000,
-                    use_stan = T, run_model = T)
-mod1stan$model_file
-summary(mod1stan)
-plot(mod1stan, 'forecast', series = 2)
+                      mu_obs = c(7, 10),
+                      train_prop = 1)
+
+# To run outside of mvgam (maybe show how other elements could be added?)
+train_dat <- sim_data$data_train %>%
+  dplyr::filter(time <= 75)
+test_dat <- sim_data$data_train %>%
+  dplyr::filter(time > 75)
+# mod1stan <- mvgam(data_train = train_dat,
+#                   data_test = test_dat,
+#                   formula = y ~ s(season, k = 12, bs = 'cc'),
+#                   knots = list(season = c(0.5, 12.5)),
+#                   family = 'nb',
+#                   trend_model = 'GP',
+#                   use_stan = T,
+#                   run_model = F)
+# test <- rstan::stan(model_code = mod1stan$model_file,
+#                     data = mod1stan$model_data,
+#                     chains = 1, iter = 1000,
+#                     init = mod1stan$inits)
+
+mod1jags <- mvgam(data_train = train_dat,
+                  data_test = test_dat,
+                  formula = y ~ s(season, k = 12, bs = 'cc'),
+                  knots = list(season = c(0.5, 12.5)),
+                  family = 'poisson',
+                  trend_model = 'AR3')
+
+mod1stan <- mvgam(data_train = train_dat,
+                  data_test = test_dat,
+                  formula = y ~ s(season, k = 12, bs = 'cc'),
+                  knots = list(season = c(0.5, 12.5)),
+                  family = 'poisson',
+                  trend_model = 'GP',
+                  use_stan = T,
+                  burnin = 1000)
+
+# Compare the out of sample DRPS for each model
+plot(mod1jags, 'forecast', series = 1, data_test = test_dat)
+plot(mod1stan, 'forecast', series = 1, data_test = test_dat)
+
+plot(mod1jags, 'forecast', series = 2, data_test = test_dat)
+plot(mod1stan, 'forecast', series = 2, data_test = test_dat)
+
+compare_mvgams(mod1jags, mod1stan, fc_horizon = 10,
+               n_evaluations = 20, n_cores = 4)
+
+plot(mod1stan, 'residuals', series = 1)
+plot(mod1stan, 'residuals', series = 2)
+plot(mod1stan, 'uncertainty', data_test = test_dat, series = 1)
+plot(mod1stan, 'uncertainty', data_test = test_dat, series = 2)
+plot(mod1stan, 'forecast', series = 1, data_test = test_dat)
+plot(mod1stan, 'forecast', series = 2, data_test = test_dat)
+plot(mod1stan, 'smooths', series = 1)
 plot(mod1stan, 'smooths', series = 2)
-plot(mod1stan, 'trend')
+plot(mod1stan, 'trend', series = 2)
+plot(sim_data$true_trends[,2])
 
-pfilter_mvgam_init(object = mod1stan, data_assim = sim_data$data_test %>%
-                     dplyr::filter(time == 54),
-                   n_particles = 10000)
+# Need to fix the particle filter so it can initiate even if only one series has a next observation
+pfilter_mvgam_init(object = mod1stan, data_assim = test_dat,
+                   n_particles = 1000)
 
-pfilter_mvgam_online(data_assim = sim_data$data_test[1:10, ], n_cores = 3, kernel_lambda = 1)
+pfilter_mvgam_online(data_assim = test_dat[1:14, ],
+                     n_cores = 3, use_resampling = T,
+                     threshold = 0.6)
 
-fcs <- pfilter_mvgam_fc(data_test = sim_data$data_test, ylim = c(0, 17))
+fcs <- pfilter_mvgam_fc(data_test = test_dat, ylim = c(0, 20), n_cores = 3)
 
-layout(matrix(c(1:2), ncol = 2, nrow = 1))
-plot_mvgam_fc(mod1stan, series = 1, data_test = sim_data$data_test,
-     ylim = c(0, 17))
+layout(matrix(c(1:2), ncol = 1, nrow = 2))
+plot_mvgam_fc(mod1stan, series = 1, data_test = test_dat,
+     ylim = c(0, 20))
 fcs$series_1()
 points(c(sim_data$data_train %>%
            dplyr::filter(series == 'series_1')%>%
            dplyr::pull(y),
          sim_data$data_test %>%
            dplyr::filter(series == 'series_1') %>%
-           dplyr::pull(y)))
+           dplyr::pull(y)),
+       pch = 16, cex = 0.4)
+
 
+layout(matrix(c(1:2), ncol = 1, nrow = 2))
+plot_mvgam_fc(mod1stan, series = 2, data_test = test_dat,
+              ylim = c(0, 20))
+fcs$series_2()
+points(c(sim_data$data_train %>%
+           dplyr::filter(series == 'series_2')%>%
+           dplyr::pull(y),
+         sim_data$data_test %>%
+           dplyr::filter(series == 'series_2') %>%
+           dplyr::pull(y)),
+       pch = 16, cex = 0.4)
+unlink('pfilter', recursive = T)
diff --git a/R/add_base_dgam_lines.R b/R/add_base_dgam_lines.R
index bcde184e..f516578d 100644
--- a/R/add_base_dgam_lines.R
+++ b/R/add_base_dgam_lines.R
@@ -67,8 +67,12 @@ add_base_dgam_lines = function(use_lv, stan = FALSE){
     tau = sigma ^ -2;
 
     // posterior predictions
-    int<lower=0> ypred[total_obs];
-    ypred = poisson_log_rng(eta + to_vector(trend));
+    matrix[n, n_series] ypred;
+    for(i in 1:n){
+    for(s in 1:n_series){
+    ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]] + trend[i, s]);
+    }
+    }
     }
     "
 
diff --git a/R/add_stan_data.R b/R/add_stan_data.R
index 9afe3e8f..d3e26d46 100644
--- a/R/add_stan_data.R
+++ b/R/add_stan_data.R
@@ -41,22 +41,22 @@ add_stan_data = function(jags_file, stan_file, jags_data, family = 'poisson'){
                'r_vec[1:n,s] = rep_vector(inv(r_inv[s]), n);\n}\n')
 
       to_negbin <- gsub('poisson_log_rng', 'neg_binomial_2_rng',
-                        stan_file[grep('ypred = poisson_log_rng', stan_file, fixed = T)])
-      stan_file[grep('ypred = poisson_log_rng', stan_file, fixed = T)] <-
-        gsub(');', ', to_vector(r_vec));', to_negbin)
+                        stan_file[grep('ypred[i, s] = poisson_log_rng', stan_file, fixed = T)])
+      stan_file[grep('ypred[i, s] = poisson_log_rng', stan_file, fixed = T)] <-
+        gsub(');', ', r_vec[i, s]);', to_negbin)
 
       add_exp_open <- gsub('\\(eta', '(exp(eta',
-                           stan_file[grep('ypred = neg_binomial', stan_file, fixed = T)])
+                           stan_file[grep('ypred[i, s] = neg_binomial', stan_file, fixed = T)])
 
-      if(any(grepl('to_vector(trend)', stan_file, fixed = T))){
-        add_exp_cl <- gsub('to_vector(trend)', 'to_vector(trend))',
+      if(any(grepl('trend[i, s]', stan_file, fixed = T))){
+        add_exp_cl <- gsub('trend[i, s]', 'trend[i, s])',
                            add_exp_open, fixed = T)
       } else {
-        add_exp_cl <- gsub('eta)', 'eta,),',
+        add_exp_cl <- gsub('eta[i, s])', 'eta[i, s],),',
                            add_exp_open)
       }
 
-      stan_file[grep('ypred = neg_binomial', stan_file, fixed = T)] <-
+      stan_file[grep('ypred[i, s] = neg_binomial', stan_file, fixed = T)] <-
         add_exp_cl
 
     stan_file <- readLines(textConnection(stan_file), n = -1)
@@ -79,6 +79,14 @@ add_stan_data = function(jags_file, stan_file, jags_data, family = 'poisson'){
                                      '; // mgcv smooth penalty matrix ', snames[i])
   }
 
+  # Get parametric prior locations and precisions if necessary
+  if('p_taus' %in% names(jags_data)){
+    p_terms <- paste0('real p_taus[', length(jags_data$p_taus),']; // prior precisions for parametric coefficients\n',
+                      'real p_coefs[', length(jags_data$p_taus), ']; // prior locations for parametric coefficients\n')
+  } else {
+    p_terms <- NULL
+  }
+
   # Search for any non-contiguous indices that sometimes are used by mgcv
   if(any(grep('in c\\(', jags_file))){
     add_idxs <- TRUE
@@ -104,40 +112,42 @@ add_stan_data = function(jags_file, stan_file, jags_data, family = 'poisson'){
 
     # Update the Stan data block
     stan_file[grep('##insert data',
-                         stan_file)] <- paste0('//Stan code generated by package mvgam\n',
-                                                     'data {',
-                                                     '\n',
-                                                     paste0(idx_data, collapse = '\n'), '\n',
-                                                     'int<lower=0> total_obs; // total number of observations\n',
-                                                     'int<lower=0> n; // number of timepoints per series\n',
-                                                     'int<lower=0> n_sp; // number of smoothing parameters\n',
-                                                     'int<lower=0> n_series; // number of series\n',
-                                                     'int<lower=0> num_basis; // total number of basis coefficients\n',
-                                                     'vector[num_basis] zero; // priro locations for basis coefficients\n',
-                                                     'matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix\n',
-                                                     'int<lower=0> ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)\n',
-                                                     paste0(smooth_penalty_data, collapse = '\n'), '\n',
-                                                     'int<lower=0, upper=1> y_observed[n, n_series]; // indices of missing vs observed\n',
-                                                     'int<lower=-1> y[n, n_series]; // time-ordered observations, with -1 indicating missing\n',
-                                                     '}\n')
+                         stan_file)] <- paste0('//Stan model code generated by package mvgam\n',
+                                               'data {',
+                                               '\n',
+                                               paste0(idx_data, collapse = '\n'), '\n',
+                                               'int<lower=0> total_obs; // total number of observations\n',
+                                               'int<lower=0> n; // number of timepoints per series\n',
+                                               'int<lower=0> n_sp; // number of smoothing parameters\n',
+                                               'int<lower=0> n_series; // number of series\n',
+                                               'int<lower=0> num_basis; // total number of basis coefficients\n',
+                                               p_terms,
+                                               'vector[num_basis] zero; // prior locations for basis coefficients\n',
+                                               'matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix\n',
+                                               'int<lower=0> ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)\n',
+                                               paste0(smooth_penalty_data, collapse = '\n'), '\n',
+                                               'int<lower=0, upper=1> y_observed[n, n_series]; // indices of missing vs observed\n',
+                                               'int<lower=-1> y[n, n_series]; // time-ordered observations, with -1 indicating missing\n',
+                                               '}\n')
   } else {
     add_idxs <- FALSE
     stan_file[grep('##insert data',
-                         stan_file)] <- paste0('//Stan code generated by package mvgam\n',
-                                                     'data {',
-                                                     '\n',
-                                                     'int<lower=0> total_obs; // total number of observations\n',
-                                                     'int<lower=0> n; // number of timepoints per series\n',
-                                                     'int<lower=0> n_sp; // number of smoothing parameters\n',
-                                                     'int<lower=0> n_series; // number of series\n',
-                                                     'int<lower=0> num_basis; // total number of basis coefficients\n',
-                                                     'vector[num_basis] zero; // prior locations for basis coefficients\n',
-                                                     'matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix\n',
-                                                     'int<lower=0> ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)\n',
-                                                     paste0(smooth_penalty_data, collapse = '\n'), '\n',
-                                                     'int<lower=0, upper=1> y_observed[n, n_series]; // indices of missing vs observed\n',
-                                                     'int<lower=-1> y[n, n_series]; // time-ordered observations, with -1 indicating missing\n',
-                                                     '}\n')
+                         stan_file)] <- paste0('//Stan model code generated by package mvgam\n',
+                                               'data {',
+                                               '\n',
+                                               'int<lower=0> total_obs; // total number of observations\n',
+                                               'int<lower=0> n; // number of timepoints per series\n',
+                                               'int<lower=0> n_sp; // number of smoothing parameters\n',
+                                               'int<lower=0> n_series; // number of series\n',
+                                               'int<lower=0> num_basis; // total number of basis coefficients\n',
+                                               'vector[num_basis] zero; // prior locations for basis coefficients\n',
+                                               p_terms,
+                                               'matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix\n',
+                                               'int<lower=0> ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)\n',
+                                               paste0(smooth_penalty_data, collapse = '\n'), '\n',
+                                               'int<lower=0, upper=1> y_observed[n, n_series]; // indices of missing vs observed\n',
+                                               'int<lower=-1> y[n, n_series]; // time-ordered observations, with -1 indicating missing\n',
+                                               '}\n')
   }
   stan_file <- readLines(textConnection(stan_file), n = -1)
 
@@ -252,7 +262,7 @@ add_stan_data = function(jags_file, stan_file, jags_data, family = 'poisson'){
              ':', as.numeric(substr(sub(".*\\:", "",
                                         stan_file[grep('// parametric effect', stan_file) + 1]),
                                     1, 1)),
-             ') {\nb_raw[i] ~ normal(0, 1);\n}')
+             ') {\nb_raw[i] ~ normal(p_coefs[i], 1 / p_taus[i]);\n}')
     stan_file <- readLines(textConnection(stan_file), n = -1)
   }
   unlink('base_gam_stan.txt')
@@ -303,6 +313,12 @@ add_stan_data = function(jags_file, stan_file, jags_data, family = 'poisson'){
     stan_data <- append(stan_data, idx_vals)
   }
 
+  # Add parametric prior means and precisions if required
+  if('p_taus' %in% names(jags_data)){
+    stan_data$p_coefs <- array(jags_data$p_coefs, dim = length(jags_data$p_taus))
+    stan_data$p_taus <- array(jags_data$p_taus, dim = length(jags_data$p_taus))
+  }
+
   return(list(stan_file = stan_file,
               model_data = stan_data))
 }
diff --git a/R/add_trend_lines.R b/R/add_trend_lines.R
index 0d617636..e81f854e 100644
--- a/R/add_trend_lines.R
+++ b/R/add_trend_lines.R
@@ -27,8 +27,8 @@ add_trend_lines = function(model_file, stan = FALSE,
         'y[i, s] ~ poisson_log(eta[ytimes[i, s]]);'
 
 
-      model_file[grep('ypred =', model_file, fixed = T)] <-
-        "ypred = poisson_log_rng(eta );"
+      model_file[grep('ypred[i, s] =', model_file, fixed = T)] <-
+        "ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]]);"
 
       model_file <- model_file[-grep('tau =', model_file)]
 
@@ -68,7 +68,7 @@ add_trend_lines = function(model_file, stan = FALSE,
                  'real<lower=0> boundary;\n',
                  'boundary = (5.0/4) * max(times);\n',
                  'int<lower=1> num_gp_basis;\n',
-                 'num_gp_basis = max(30, n);\n',
+                 'num_gp_basis = max(40, n);\n',
                  'matrix[n, num_gp_basis] gp_phi;\n',
                  'for (m in 1:num_gp_basis){\n',
                  'gp_phi[,m] = phi_SE(boundary, m, times);\n',
@@ -76,7 +76,8 @@ add_trend_lines = function(model_file, stan = FALSE,
                  'parameters {')
 
         model_file[grep('##insert data', model_file)-1] <-
-          paste0('functions {\n',
+          paste0('//Stan GP functions generated by package mvgam\n',
+                 'functions {\n',
                  'real lambda_gp(real L, int m) {\n',
                  'real lam;\n',
                  'lam = ((m*pi())/(2*L))^2;\n',
diff --git a/R/eval_mvgam.R b/R/eval_mvgam.R
index 519f959d..b1a1583f 100644
--- a/R/eval_mvgam.R
+++ b/R/eval_mvgam.R
@@ -55,12 +55,12 @@ eval_mvgam = function(object,
     }
   }
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
+  if(eval_timepoint < 3){
+    stop('argument "eval_timepoint" must be >= 3',
+         call. = FALSE)
   }
 
+
   #### 1. Generate linear predictor matrix for covariates and extract trend estimates at timepoint
   data_train <- object$obs_data
   n_series <- NCOL(object$ytimes)
@@ -161,12 +161,37 @@ eval_mvgam = function(object,
     starts <- ends + 1
     starts <- c(1, starts[-c(1, (n_series + 1))])
     ends <- ends[-1]
-    trends <- lapply(seq_len(n_series), function(series){
-      trend_estimates <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
-      as.matrix(trend_estimates[,(eval_timepoint - 2):eval_timepoint])
-    })
 
-    taus <- MCMCvis::MCMCchains(object$model_output, 'tau')
+    if(object$trend_model == 'GP'){
+      # Pull out entire series of trend estimates up to the evaluation timepoint if this is a GP trend
+      trends <- lapply(seq_len(n_series), function(series){
+        MCMCvis::MCMCchains(object$model_output, 'trend')[,seq(series,
+                                                               dim(MCMCvis::MCMCchains(object$model_output,
+                                                                                       'trend'))[2],
+                                                               by = NCOL(object$ytimes))]
+      })
+
+      taus <- NULL
+    } else {
+      # Otherwise only keep the preceeding three trend estimates to save memory
+      if(object$fit_engine == 'stan'){
+        trends <- lapply(seq_len(n_series), function(series){
+          trend_estimates <- MCMCvis::MCMCchains(object$model_output, 'trend')[,seq(series,
+                                                                                    dim(MCMCvis::MCMCchains(object$model_output,
+                                                                                                            'trend'))[2],
+                                                                                    by = NCOL(object$ytimes))]
+          as.matrix(trend_estimates[,(eval_timepoint - 2):eval_timepoint])
+        })
+      } else {
+        trends <- lapply(seq_len(n_series), function(series){
+          trend_estimates <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+          as.matrix(trend_estimates[,(eval_timepoint - 2):eval_timepoint])
+        })
+      }
+
+      taus <- MCMCvis::MCMCchains(object$model_output, 'tau')
+    }
+
     }
 
     lvs <- NULL
@@ -198,7 +223,17 @@ eval_mvgam = function(object,
   }
 
   # AR term estimates
+  if(object$trend_model == 'GP'){
+    alpha_gps <- MCMCvis::MCMCchains(object$model_output, 'alpha_gp')
+    rho_gps <- MCMCvis::MCMCchains(object$model_output, 'rho_gp')
+    ar1s <- NULL
+    ar2s <- NULL
+    ar3s <- NULL
+  }
+
   if(object$trend_model %in% c('RW', 'None')){
+    alpha_gps <- NULL
+    rho_gps <- NULL
     if(object$use_lv){
       ar1s <- matrix(0, nrow = NROW(betas), ncol = object$n_lv)
       ar2s <- matrix(0, nrow = NROW(betas), ncol = object$n_lv)
@@ -212,6 +247,8 @@ eval_mvgam = function(object,
 
   if(object$trend_model == 'AR1'){
     ar1s <- MCMCvis::MCMCchains(object$model_output, 'ar1')
+    alpha_gps <- NULL
+    rho_gps <- NULL
 
     if(object$use_lv){
       ar2s <- matrix(0, nrow = NROW(betas), ncol = object$n_lv)
@@ -225,6 +262,8 @@ eval_mvgam = function(object,
   if(object$trend_model == 'AR2'){
     ar1s <- MCMCvis::MCMCchains(object$model_output, 'ar1')
     ar2s <- MCMCvis::MCMCchains(object$model_output, 'ar2')
+    alpha_gps <- NULL
+    rho_gps <- NULL
 
     if(object$use_lv){
       ar3s <- matrix(0, nrow = NROW(betas), ncol = object$n_lv)
@@ -234,6 +273,8 @@ eval_mvgam = function(object,
   }
 
   if(object$trend_model == 'AR3'){
+    alpha_gps <- NULL
+    rho_gps <- NULL
     ar1s <- MCMCvis::MCMCchains(object$model_output, 'ar1')
     ar2s <- MCMCvis::MCMCchains(object$model_output, 'ar2')
     ar3s <- MCMCvis::MCMCchains(object$model_output, 'ar3')
@@ -266,6 +307,7 @@ eval_mvgam = function(object,
   #### 2. Run trends forward fc_horizon steps to generate the forecast distribution ####
   use_lv <- object$use_lv
   upper_bounds <- object$upper_bounds
+  trend_model <- object$trend_model
 
   # Function to simulate trends / latent factors ahead using ar3 model
   sim_ar3 = function(phi, ar1, ar2, ar3, tau, state, h){
@@ -281,6 +323,21 @@ eval_mvgam = function(object,
     states[-c(1:3)]
   }
 
+  # Function to simulate trends ahead using squared exponential GP
+  sim_gp = function(alpha_gp, rho_gp, state, h){
+    t <- 1:length(state)
+    t_new <- 1:(length(state) + h)
+
+    Sigma_new <- alpha_gp^2 * exp(- outer(t, t_new, "-")^2 / (2 * rho_gp^2))
+    Sigma_star <- alpha_gp^2 * exp(- outer(t_new, t_new, "-")^2 / (2 * rho_gp^2))
+    Sigma <- alpha_gp^2 * exp(- outer(t, t, "-")^2 / (2 * rho_gp^2)) +
+      diag(1e-4, length(state))
+
+    tail(t(Sigma_new) %*% solve(Sigma, state), h) +
+      tail(MASS::mvrnorm(1, mu = rep(0, dim(Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new))[2]),
+                         Sigma = Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new)), h)
+  }
+
   # Run particles forward in time to generate their forecasts
   if(n_cores > 1){
 
@@ -298,6 +355,9 @@ eval_mvgam = function(object,
                         'ar1s',
                         'ar2s',
                         'ar3s',
+                        'alpha_gps',
+                        'rho_gps',
+                        'trend_model',
                         'lvs',
                         'lv_coefs',
                         'trends',
@@ -306,7 +366,8 @@ eval_mvgam = function(object,
                         'twdis',
                         'n_series',
                         'upper_bounds',
-                        'sim_ar3'),
+                        'sim_ar3',
+                        'sim_gp'),
                 envir = environment())
 
   pbapply::pboptions(type = "none")
@@ -398,14 +459,19 @@ eval_mvgam = function(object,
         trends[[trend]][samp_index, ]
       })
 
-      # Sample AR parameters
-      phi <- phis[samp_index, ]
-      ar1 <- ar1s[samp_index, ]
-      ar2 <- ar2s[samp_index, ]
-      ar3 <- ar3s[samp_index, ]
+      if(trend_model == 'GP'){
+        alpha_gp <- alpha_gps[samp_index, ]
+        rho_gp <- rho_gps[samp_index, ]
+      } else {
+        # Sample AR parameters
+        phi <- phis[samp_index, ]
+        ar1 <- ar1s[samp_index, ]
+        ar2 <- ar2s[samp_index, ]
+        ar3 <- ar3s[samp_index, ]
 
-      # Sample trend precisions
-      tau <- taus[samp_index,]
+        # Sample trend precisions
+        tau <- taus[samp_index,]
+      }
 
       # Family-specific parameters
       if(family == 'Negative Binomial'){
@@ -418,13 +484,21 @@ eval_mvgam = function(object,
       }
 
       series_fcs <- lapply(seq_len(n_series), function(series){
-        trend_preds <- sim_ar3(phi = phi[series],
-                               ar1 = ar1[series],
-                               ar2 = ar2[series],
-                               ar3 = ar3[series],
-                               tau = tau[series],
-                               state = last_trends[[series]],
-                               h = fc_horizon)
+        if(trend_model == 'GP'){
+          trend_preds <- sim_gp(alpha_gp = alpha_gp[series],
+                                 rho_gp = rho_gp[series],
+                                 state = last_trends[[series]],
+                                 h = fc_horizon)
+        } else {
+          trend_preds <- sim_ar3(phi = phi[series],
+                                 ar1 = ar1[series],
+                                 ar2 = ar2[series],
+                                 ar3 = ar3[series],
+                                 tau = tau[series],
+                                 state = last_trends[[series]],
+                                 h = fc_horizon)
+        }
+
         if(family == 'Negative Binomial'){
           fc <-  rnbinom(fc_horizon,
                          mu = exp(as.vector((Xp[which(as.numeric(data_assim$series) == series),] %*%
@@ -544,14 +618,19 @@ eval_mvgam = function(object,
           trends[[trend]][samp_index, ]
         })
 
-        # Sample AR parameters
-        phi <- phis[samp_index, ]
-        ar1 <- ar1s[samp_index, ]
-        ar2 <- ar2s[samp_index, ]
-        ar3 <- ar3s[samp_index, ]
-
-        # Sample trend precisions
-        tau <- taus[samp_index,]
+        if(trend_model == 'GP'){
+          alpha_gp <- alpha_gps[samp_index, ]
+          rho_gp <- rho_gps[samp_index, ]
+        } else {
+          # Sample AR parameters
+          phi <- phis[samp_index, ]
+          ar1 <- ar1s[samp_index, ]
+          ar2 <- ar2s[samp_index, ]
+          ar3 <- ar3s[samp_index, ]
+
+          # Sample trend precisions
+          tau <- taus[samp_index,]
+        }
 
         # Family-specific parameters
         if(family == 'Negative Binomial'){
@@ -564,13 +643,21 @@ eval_mvgam = function(object,
         }
 
         series_fcs <- lapply(seq_len(n_series), function(series){
-          trend_preds <- sim_ar3(phi = phi[series],
-                                 ar1 = ar1[series],
-                                 ar2 = ar2[series],
-                                 ar3 = ar3[series],
-                                 tau = tau[series],
-                                 state = last_trends[[series]],
-                                 h = fc_horizon)
+          if(trend_model == 'GP'){
+            trend_preds <- sim_gp(alpha_gp = alpha_gp[series],
+                                  rho_gp = rho_gp[series],
+                                  state = last_trends[[series]],
+                                  h = fc_horizon)
+          } else {
+            trend_preds <- sim_ar3(phi = phi[series],
+                                   ar1 = ar1[series],
+                                   ar2 = ar2[series],
+                                   ar3 = ar3[series],
+                                   tau = tau[series],
+                                   state = last_trends[[series]],
+                                   h = fc_horizon)
+          }
+
           if(family == 'Negative Binomial'){
             fc <-  rnbinom(fc_horizon,
                            mu = exp(as.vector((Xp[which(as.numeric(data_assim$series) == series),] %*%
@@ -613,7 +700,7 @@ eval_mvgam = function(object,
 
   # Evaluate against the truth
   series_truths <- lapply(seq_len(n_series), function(series){
-    if(class(object$obs_data) == 'list'){
+    if(class(object$obs_data)[1] == 'list'){
       data_assim[['y']][which(as.numeric(data_assim$series) == series)]
     } else {
       data_assim[which(as.numeric(data_assim$series) == series),'y']
@@ -662,5 +749,4 @@ eval_mvgam = function(object,
   names(series_drps) <- levels(data_assim$series)
 
   return(series_drps)
-
 }
diff --git a/R/get_mvgam_resids.R b/R/get_mvgam_resids.R
index cd4a9a9c..a397b506 100644
--- a/R/get_mvgam_resids.R
+++ b/R/get_mvgam_resids.R
@@ -9,12 +9,8 @@
 #'@return A \code{list} of residual distributions
 get_mvgam_resids = function(object, n_cores = 1){
 
-  # Convert stanfit objects to coda samples
-  if(object$fit_engine == 'stan'){
-    preds <- rstan::extract(object$model_output, 'ypred')[[1]]
-  } else {
-    preds <- MCMCvis::MCMCchains(object$model_output, 'ypred')
-  }
+
+preds <- MCMCvis::MCMCchains(object$model_output, 'ypred')
 
 # Functions for calculating randomised quantile (Dunn-Smyth) residuals
 ds_resids_nb = function(truth, fitted, draw, size){
@@ -125,7 +121,14 @@ series_resids <- pbapply::pblapply(seq_len(NCOL(object$ytimes)), function(series
       dplyr::filter(series == !!(levels(object$obs_data$series)[series])) %>%
       nrow()
   }
-  preds <- preds[,starts[series]:ends[series]]
+
+  if(object$fit_engine == 'stan'){
+    preds <- preds[,seq(series,
+                        dim(MCMCvis::MCMCchains(object$model_output, 'ypred'))[2],
+                        by = NCOL(object$ytimes))]
+  } else {
+    preds <- preds[,starts[series]:ends[series]]
+  }
 
   if(class(object$obs_data)[1] == 'list'){
     obj_dat <- data.frame(y = object$obs_data$y,
diff --git a/R/mvgam.R b/R/mvgam.R
index e60688db..a7a5e264 100644
--- a/R/mvgam.R
+++ b/R/mvgam.R
@@ -126,7 +126,17 @@
 #'If the fitted model is appropriate then Dunn-Smyth residuals will be standard normal in distribution and no
 #'autocorrelation will be evident. When a particular observation is missing, the residual is calculated by comparing independent
 #'draws from the model's posterior distribution
-#'
+#'\cr
+#'\cr
+#'*Using Stan*: An in-development feature of `mvgam` is the ability to use Hamiltonian Monte Carlo for parameter estimation
+#'via the software `Stan` (using the `rstan` interface). Note that the `rstan` library is currently required for this option to work,
+#'though support for other `Stan` interfaces will be added in future. Also note that currently there is no support for
+#'fitting `Tweedie` responses or dynamic factor models in `Stan`, though again these will be added in future.
+#'However there are great advantages when using `Stan`, which includes the option to estimate smooth latent trends
+#'via [Hilbert space approximate Gaussian Processes](https://arxiv.org/abs/2004.11408). This often makes sense for
+#'ecological series, which we expect to change smoothly. In `mvgam`, latent squared exponential GP trends are approximated using
+#'by default \code{40} basis functions, which saves computational costs compared to fitting full GPs while adequately estimating
+#'GP \code{alpha} and \code{rho} parameters
 #'@author Nicholas J Clark
 #'
 #'@seealso \code{\link[rjags]{jags.model}}, \code{\link[mcgv]{jagam}}, \code{\link[mcgv]{gam}}
@@ -718,7 +728,7 @@ mvgam = function(formula,
     parametric_tdata <- NULL
   }
 
-  model_file <- c('JAGS code generated by package mvgam',
+  model_file <- c('JAGS model code generated by package mvgam',
                        '\n',
                        'GAM formula:',
                   gsub('\"', '', paste(formula[2],formula[3],sep=' ~ ')),
@@ -765,6 +775,19 @@ mvgam = function(formula,
                                     jags_data = ss_jagam$jags.data,
                                     family = family)
 
+      # Sensible inits needed for the betas
+      if(trend_model != 'None'){
+        initf1 <- function() {
+          list(b_raw = runif(stan_objects$model_data$num_basis, -0.25, 0.25))
+        }
+      } else {
+        initf1 <- function() {
+          list(b_raw = runif(stan_objects$model_data$num_basis, -0.25, 0.25),
+               sigma = runif(stan_objects$model_data$n_series, 0.075, 1))
+        }
+      }
+      inits <- initf1
+
       output <- structure(list(call = formula,
                                family = dplyr::case_when(family == 'tw' ~ 'Tweedie',
                                                          family == 'poisson' ~ 'Poisson',
@@ -774,6 +797,7 @@ mvgam = function(formula,
                                pregam = ss_jagam$pregam,
                                model_file = stan_objects$stan_file,
                                model_data = stan_objects$model_data,
+                               inits = inits,
                                mgcv_model = ss_gam,
                                sp_names = names(ss_jagam$pregam$lsp0),
                                ytimes = ytimes,
@@ -784,6 +808,9 @@ mvgam = function(formula,
                                fit_engine = 'stan'),
                           class = 'mvgam_prefit')
     } else {
+      initlist <- replicate(chains, ss_jagam$jags.ini,
+                            simplify = FALSE)
+      inits <- initlist
       output <- structure(list(call = formula,
                                family = dplyr::case_when(family == 'tw' ~ 'Tweedie',
                                                          family == 'poisson' ~ 'Poisson',
@@ -793,6 +820,7 @@ mvgam = function(formula,
                                pregam = ss_jagam$pregam,
                                model_file = trimws(model_file),
                                model_data = ss_jagam$jags.data,
+                               inits = inits,
                                mgcv_model = ss_gam,
                                sp_names = names(ss_jagam$pregam$lsp0),
                                ytimes = ytimes,
@@ -834,7 +862,7 @@ mvgam = function(formula,
       require(rstan)
       options(mc.cores = parallel::detectCores())
 
-      # Sensible inits needed for the betas
+      # Sensible inits needed for the betas (and trend alphas??)
       if(trend_model != 'None'){
         initf1 <- function() {
           list(b_raw = runif(stan_objects$model_data$num_basis, -0.25, 0.25))
@@ -854,12 +882,14 @@ mvgam = function(formula,
         stan_control <- list(max_treedepth = 10)
       }
 
+      message("Compiling the Stan model")
+      message()
       fit1 <- stan(model_code = stan_objects$stan_file,
                    iter = burnin + 1000,
                    chains = chains,
                    data = stan_objects$model_data,
                    cores = min(c(chains, parallel::detectCores() - 1)),
-                   init = initf1,
+                   init = inits,
                    verbose = F,
                    control = stan_control,
                    pars = get_monitor_pars(family = family,
@@ -893,6 +923,10 @@ mvgam = function(formula,
       # is more crucial than little adaptation but long 'burnin' https://mmeredith.net/blog/2016/Adapt_or_burn.htm
       unlink('base_gam.txt')
       cat(model_file, file = 'base_gam.txt', sep = '\n', append = T)
+
+      message("Compiling the JAGS model")
+      message()
+
       if(parallel){
         cl <- parallel::makePSOCKcluster(min(c(chains, parallel::detectCores() - 1)))
         setDefaultCluster(cl)
@@ -908,7 +942,7 @@ mvgam = function(formula,
                                      jags = jags_path,
                                      thin = thin,
                                      method = "rjparallel",
-                                     monitor = c(param, 'deviance'),
+                                     monitor = param,
                                      silent.jags = TRUE,
                                      cl = cl)
         stopCluster(cl)
@@ -925,7 +959,7 @@ mvgam = function(formula,
                                      jags = jags_path,
                                      thin = thin,
                                      method = "rjags",
-                                     monitor = c(param, 'deviance'),
+                                     monitor = param,
                                      silent.jags = TRUE)
       }
       out_gam_mod <- coda::as.mcmc.list(gam_mod)
diff --git a/R/pfilter_mvgam_fc.R b/R/pfilter_mvgam_fc.R
index 4fd8e44e..bd654f4d 100644
--- a/R/pfilter_mvgam_fc.R
+++ b/R/pfilter_mvgam_fc.R
@@ -62,17 +62,14 @@ pfilter_mvgam_fc = function(file_path = 'pfilter',
     t <- 1:length(state)
     t_new <- 1:(length(state) + h)
 
-    # Evaluate on a fine a grid of points to preserve the
-    # infinite dimensionality
-    scale_down <- (100 / length(state)) * rho_gp
-    t <- t / scale_down
-    t_new <- t_new / scale_down
-
-    Sigma_new <- alpha_gp^2 * exp(-(rho_gp/scale_down) * outer(t, t_new, "-")^2)
-    Sigma <- alpha_gp^2 * exp(-(rho_gp/scale_down) * outer(t, t, "-")^2) +
+    Sigma_new <- alpha_gp^2 * exp(- outer(t, t_new, "-")^2 / (2 * rho_gp^2))
+    Sigma_star <- alpha_gp^2 * exp(- outer(t_new, t_new, "-")^2 / (2 * rho_gp^2))
+    Sigma <- alpha_gp^2 * exp(- outer(t, t, "-")^2 / (2 * rho_gp^2)) +
       diag(1e-4, length(state))
 
-    tail(t(Sigma_new) %*% solve(Sigma, state), h)
+    tail(t(Sigma_new) %*% solve(Sigma, state), h) +
+      tail(MASS::mvrnorm(1, mu = rep(0, dim(Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new))[2]),
+                         Sigma = Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new)), h)
   }
 
   if(missing(ylim)){
diff --git a/R/pfilter_mvgam_init.R b/R/pfilter_mvgam_init.R
index 028b4779..97894ed7 100644
--- a/R/pfilter_mvgam_init.R
+++ b/R/pfilter_mvgam_init.R
@@ -29,12 +29,6 @@ pfilter_mvgam_init = function(object,
     stop('argument "object" must be of class "mvgam"')
   }
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
-  }
-
 #### 1. Generate linear predictor matrix for the next timepoint and extract last trend estimates
 # (NOTE, all series must have observations for the next timepoint, even if they are NAs!!!!) ####
 data_train <- object$obs_data
@@ -124,7 +118,14 @@ if(object$use_lv){
   starts <- c(1, starts[-c(1, (n_series + 1))])
   ends <- ends[-1]
   trends <- lapply(seq_len(n_series), function(series){
-    trend_estimates <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+    if(object$fit_engine == 'stan'){
+      trend_estimates <- MCMCvis::MCMCchains(object$model_output, 'trend')[,seq(series,
+                                                                                dim(MCMCvis::MCMCchains(object$model_output,
+                                                                                                        'trend'))[2],
+                                                                                by = NCOL(object$ytimes))]
+    } else {
+      trend_estimates <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+    }
 
     # Need to only use estimates from the training period
     end_train <- object$obs_data %>%
@@ -435,17 +436,15 @@ particles <- pbapply::pblapply(sample_seq, function(x){
         t <- 1:length(last_trends[[trend]])
         t_new <- 1:(length(last_trends[[trend]]) + 1)
 
-        # Evaluate on a fine a grid of points to preserve the
-        # infinite dimensionality
-        scale_down <- (100 / length(last_trends[[trend]])) * rho_gp[trend]
-        t <- t / scale_down
-        t_new <- t_new / scale_down
 
-        Sigma_new <- alpha_gp[trend]^2 * exp(-(rho_gp[trend]/scale_down) * outer(t, t_new, "-")^2)
-        Sigma <- alpha_gp[trend]^2 * exp(-(rho_gp[trend]/scale_down) * outer(t, t, "-")^2) +
+        Sigma_new <- alpha_gp[trend]^2 * exp(- outer(t, t_new, "-")^2 / (2 * rho_gp[trend]^2))
+        Sigma_star <- alpha_gp[trend]^2 * exp(- outer(t_new, t_new, "-")^2 / (2 * rho_gp[trend]^2))
+        Sigma <- alpha_gp[trend]^2 * exp(- outer(t, t, "-")^2 / (2 * rho_gp[trend]^2)) +
           diag(1e-4, length(last_trends[[trend]]))
 
-        tail(t(Sigma_new) %*% solve(Sigma, last_trends[[trend]]), 1)
+        tail(t(Sigma_new) %*% solve(Sigma, last_trends[[trend]]), 1) +
+          tail(MASS::mvrnorm(1, mu = rep(0, dim(Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new))[2]),
+                             Sigma = Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new)), 1)
       }))
     }
 
@@ -502,7 +501,6 @@ particles <- pbapply::pblapply(sample_seq, function(x){
                      alpha_gp = as.numeric(alpha_gp),
                      rho_gp = as.numeric(rho_gp),
                      trend_states = lapply(last_trends, unname),
-                     trend_states = trends,
                      weight = weight,
                      liks = liks,
                      upper_bounds = upper_bounds,
@@ -542,7 +540,6 @@ particles <- pbapply::pblapply(sample_seq, function(x){
                      alpha_gp = as.numeric(alpha_gp),
                      rho_gp = as.numeric(rho_gp),
                      trend_states = lapply(last_trends, unname),
-                     trend_states = trends,
                      weight = weight,
                      liks = liks,
                      upper_bounds = upper_bounds,
@@ -584,7 +581,6 @@ particles <- pbapply::pblapply(sample_seq, function(x){
                      alpha_gp = as.numeric(alpha_gp),
                      rho_gp = as.numeric(rho_gp),
                      trend_states = lapply(last_trends, unname),
-                     trend_states = trends,
                      weight = weight,
                      liks = liks,
                      upper_bounds = upper_bounds,
diff --git a/R/pfilter_mvgam_online.R b/R/pfilter_mvgam_online.R
index 4fa13096..5a599117 100644
--- a/R/pfilter_mvgam_online.R
+++ b/R/pfilter_mvgam_online.R
@@ -17,7 +17,9 @@
 #'the specified \code{threshold}. Default for this option is \code{FALSE}, relying instead on kernel smoothing only
 #'to maintain particle diversity
 #'@param kernel_lambda \code{proportional numeric} specifying the strength of kernel smoothing to use when
-#'pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}
+#'pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}. Note
+#'however that kernel smoothing is not currently supported for Gaussian Process trend models. For these models,
+#'it is still preferable to use resampling with a reasonable \code{threshold} value (i.e. \code{0.5})
 #'@param file_path \code{character} string specifying the file path for locating the particles
 #'@param n_cores \code{integer} specifying number of cores for generating particle forecasts in parallel
 #'@return A \code{list} object of \code{length = n_particles} containing information on parameters and
diff --git a/R/pfilter_mvgam_smooth.R b/R/pfilter_mvgam_smooth.R
index c962585e..b7ab07c7 100644
--- a/R/pfilter_mvgam_smooth.R
+++ b/R/pfilter_mvgam_smooth.R
@@ -22,11 +22,13 @@
 #'the specified \code{threshold}. Note that resampling can result in loss of the original model's diversity of
 #'GAM beta coefficients, which may have undesirable consequences for the forecast distribution. If
 #'\code{use_resampling} is \code{TRUE}, some effort is made to remedy this by assigning randomly sampled draws of
-#'beta coefficients from the original model's distribution to each particle. This does not however guarantee that there
-#'will be no loss of diversity, especially if successive resamples take place. Default for this option is therefore
+#'GAM beta coefficients from the original model's distribution to each particle. This does not however guarantee there
+#'will be no loss of diversity, especially if successive resampling take place. Default for this option is therefore
 #'\code{FALSE}
 #'@param kernel_lambda \code{proportional numeric} specifying the strength of smoothing to use when
-#'pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}
+#'pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}. Note
+#'however that kernel smoothing is not currently supported for Gaussian Process trend models. For these models,
+#'it is still preferable to use resampling with a reasonable \code{threshold} value (i.e. \code{0.5})
 #'@param file_path \code{character} string specifying the file path for locating the particles
 #'@param n_cores \code{integer} specifying number of cores for generating particle forecasts in parallel
 #'@return A \code{list} object of \code{length = n_particles} containing information on parameters and
@@ -93,18 +95,17 @@ pfilter_mvgam_smooth = function(particles,
           t <- 1:length(particles[[x]]$trend_states[[series]])
           t_new <- 1:(length(particles[[x]]$trend_states[[series]])+1)
 
-          # Evaluate on a fine a grid of points to preserve the
-          # infinite dimensionality
-          scale_down <- (100 / length(particles[[x]]$trend_states[[series]])) * particles[[x]]$rho_gp[series]
-          t <- t / scale_down
-          t_new <- t_new / scale_down
-
-          Sigma_new <- particles[[x]]$alpha_gp[series]^2 * exp(-(particles[[x]]$rho_gp[series]/scale_down) * outer(t, t_new, "-")^2)
+          Sigma_new <- particles[[x]]$alpha_gp[series]^2 *
+            exp(- outer(t, t_new, "-")^2 / (2 * particles[[x]]$rho_gp[series]^2))
+          Sigma_star <- particles[[x]]$alpha_gp[series]^2 *
+            exp(- outer(t_new, t_new, "-")^2 / (2 * particles[[x]]$rho_gp[series]^2))
           Sigma <- particles[[x]]$alpha_gp[series]^2 *
-            exp(-(particles[[x]]$rho_gp[series]/scale_down) * outer(t, t, "-")^2) +
+            exp(- outer(t, t, "-")^2 / (2 * particles[[x]]$rho_gp[series]^2)) +
             diag(1e-4, length(particles[[x]]$trend_states[[series]]))
 
-          tail(t(Sigma_new) %*% solve(Sigma, particles[[x]]$trend_states[[series]]), 1)
+          tail(t(Sigma_new) %*% solve(Sigma, particles[[x]]$trend_states[[series]]), 1) +
+            tail(MASS::mvrnorm(1, mu = rep(0, dim(Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new))[2]),
+                               Sigma = Sigma_star - t(Sigma_new) %*% solve(Sigma, Sigma_new)), 1)
 
         }))
       } else {
@@ -491,6 +492,7 @@ pfilter_mvgam_smooth = function(particles,
     orig_betas <- NULL
   } else {
     # If resampling, must preserve the original diversity of GAM beta coefficients
+    cat('Resampling particles ...\n')
     next_update_seq <- index
     orig_betas <- do.call(rbind, purrr::map(particles, 'betas'))
   }
@@ -661,7 +663,11 @@ pfilter_mvgam_smooth = function(particles,
     }
 
     # Return the updated particle, preserving original betas
-    betas <- particles[[x]]$betas
+    if(use_resampling){
+      betas <- orig_betas[x,]
+    } else {
+      betas <- particles[[x]]$betas
+    }
 
     if(particles[[x]]$family == 'Negative Binomial'){
       if(particles[[x]]$trend_model == 'GP'){
diff --git a/R/plot_mvgam_fc.R b/R/plot_mvgam_fc.R
index f977359d..222db527 100644
--- a/R/plot_mvgam_fc.R
+++ b/R/plot_mvgam_fc.R
@@ -40,7 +40,11 @@ plot_mvgam_fc = function(object, series = 1, data_test, hide_xlabels = FALSE, yl
   ends <- ends[-1]
 
   if(object$fit_engine == 'stan'){
-    preds <- rstan::extract(object$model_output, 'ypred')[[1]][,starts[series]:ends[series]]
+
+    # For stan objects, ypred is stored as a vector in column-major order
+    preds <- MCMCvis::MCMCchains(object$model_output, 'ypred')[,seq(series,
+                                                                    dim(MCMCvis::MCMCchains(object$model_output, 'ypred'))[2],
+                                                                    by = NCOL(object$ytimes))]
   } else {
     preds <- MCMCvis::MCMCchains(object$model_output, 'ypred')[,starts[series]:ends[series]]
   }
@@ -130,6 +134,7 @@ plot_mvgam_fc = function(object, series = 1, data_test, hide_xlabels = FALSE, yl
                                time = data_test$time)
     }
 
+    # Plot training and testing points
     points(dplyr::bind_rows(data_train, data_test) %>%
              dplyr::filter(series == s_name) %>%
              dplyr::select(time, y) %>%
@@ -143,6 +148,56 @@ plot_mvgam_fc = function(object, series = 1, data_test, hide_xlabels = FALSE, yl
              dplyr::arrange(time) %>%
              dplyr::pull(y), pch = 16, col = "black", cex = 0.55)
     abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed')
+
+    # Calculate out of sample DRPS and print the score
+    drps_score <- function(truth, fc, interval_width = 0.9){
+      nsum <- 1000
+      Fy = ecdf(fc)
+      ysum <- 0:nsum
+      indicator <- ifelse(ysum - truth >= 0, 1, 0)
+      score <- sum((indicator - Fy(ysum))^2)
+
+      # Is value within empirical interval?
+      interval <- quantile(fc, probs = c((1-interval_width)/2, (interval_width + (1-interval_width)/2)))
+      in_interval <- ifelse(truth <= interval[2] & truth >= interval[1], 1, 0)
+      return(c(score, in_interval))
+    }
+
+    # Wrapper to operate on all observations in fc_horizon
+    drps_mcmc_object <- function(truth, fc, interval_width = 0.9){
+      indices_keep <- which(!is.na(truth))
+      if(length(indices_keep) == 0){
+        scores = data.frame('drps' = rep(NA, length(truth)),
+                            'interval' = rep(NA, length(truth)))
+      } else {
+        scores <- matrix(NA, nrow = length(truth), ncol = 2)
+        for(i in indices_keep){
+          scores[i,] <- drps_score(truth = as.vector(truth)[i],
+                                   fc = fc[,i], interval_width)
+        }
+      }
+      scores
+    }
+
+    truth <- as.matrix(data_test %>%
+                         dplyr::filter(series == s_name) %>%
+                         dplyr::select(time, y) %>%
+                         dplyr::distinct() %>%
+                         dplyr::arrange(time) %>%
+                         dplyr::pull(y))
+    last_train <- length(data_train %>%
+      dplyr::filter(series == s_name) %>%
+      dplyr::select(time, y) %>%
+      dplyr::distinct() %>%
+      dplyr::arrange(time) %>%
+      dplyr::pull(y))
+
+    fc <- preds[,(last_train+1):NCOL(preds)]
+
+    message('Out of sample DRPS:')
+    print(sum(drps_mcmc_object(as.vector(truth),
+                                fc)[,1]))
+    message()
   } else {
     if(class(data_train)[1] == 'list'){
       data_train <- data.frame(series = data_train$series,
diff --git a/R/plot_mvgam_resids.R b/R/plot_mvgam_resids.R
index 6851c6c4..285e440b 100644
--- a/R/plot_mvgam_resids.R
+++ b/R/plot_mvgam_resids.R
@@ -19,12 +19,6 @@ plot_mvgam_resids = function(object, series = 1, n_bins = 25){
     stop('argument "object" must be of class "mvgam"')
   }
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
-  }
-
   if(sign(series) != 1){
     stop('argument "series" must be a positive integer',
          call. = FALSE)
diff --git a/R/plot_mvgam_trace.R b/R/plot_mvgam_trace.R
index c2d2ad64..e74ccbb0 100644
--- a/R/plot_mvgam_trace.R
+++ b/R/plot_mvgam_trace.R
@@ -22,12 +22,6 @@ plot_mvgam_trace = function(object, param = 'rho', overlay_prior = TRUE){
     stop('argument "object" must be of class "mvgam"')
   }
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
-  }
-
   if(param == 'rho'){
     param_names <- object$sp_names
     prior_mat <- matrix(NA, nrow = dim(MCMCvis::MCMCchains(object$model_output, 'rho'))[1],
diff --git a/R/plot_mvgam_trend.R b/R/plot_mvgam_trend.R
index 1e927ae1..2bb30655 100644
--- a/R/plot_mvgam_trend.R
+++ b/R/plot_mvgam_trend.R
@@ -16,12 +16,6 @@ plot_mvgam_trend = function(object, series = 1, data_test,
     stop('argument "object" must be of class "mvgam"')
   }
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
-  }
-
   if(sign(series) != 1){
     stop('argument "series" must be a positive integer',
          call. = FALSE)
@@ -40,7 +34,14 @@ plot_mvgam_trend = function(object, series = 1, data_test,
   starts <- c(1, starts[-c(1, (NCOL(object$ytimes)+1))])
   ends <- ends[-1]
 
-  preds <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+  if(object$fit_engine == 'stan'){
+    preds <- MCMCvis::MCMCchains(object$model_output, 'trend')[,seq(series,
+                                                           dim(MCMCvis::MCMCchains(object$model_output,
+                                                                                   'trend'))[2],
+                                                           by = NCOL(object$ytimes))]
+  } else {
+    preds <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+  }
   preds_last <- preds[1,]
   pred_vals <- seq(1:length(preds_last))
 
diff --git a/R/plot_mvgam_uncertainty.R b/R/plot_mvgam_uncertainty.R
index 15b4807f..154aece5 100644
--- a/R/plot_mvgam_uncertainty.R
+++ b/R/plot_mvgam_uncertainty.R
@@ -17,12 +17,6 @@ plot_mvgam_uncertainty = function(object, series = 1, data_test, legend_position
     stop('argument "object" must be of class "mvgam"')
   }
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
-  }
-
   if(sign(series) != 1){
     stop('argument "series" must be a positive integer',
          call. = FALSE)
@@ -81,7 +75,14 @@ plot_mvgam_uncertainty = function(object, series = 1, data_test, legend_position
   betas <- MCMCvis::MCMCchains(object$model_output, 'b')
 
   # Extract current trend estimates
-  trend <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+  if(object$fit_engine == 'stan'){
+    trend <- MCMCvis::MCMCchains(object$model_output, 'trend')[,seq(series,
+                                                                    dim(MCMCvis::MCMCchains(object$model_output,
+                                                                                            'trend'))[2],
+                                                                    by = NCOL(object$ytimes))]
+  } else {
+    trend <- MCMCvis::MCMCchains(object$model_output, 'trend')[,starts[series]:ends[series]]
+  }
 
   if(length(unique(data_train$series)) == 1){
     trend <- matrix(trend[, NCOL(trend)])
@@ -94,7 +95,6 @@ plot_mvgam_uncertainty = function(object, series = 1, data_test, legend_position
     }
   }
 
-
   # Function to calculate intersection of two uncertainty distributions
   intersect_hist = function(fullpreds, gampreds){
     from <- min(min(fullpreds, na.rm = T),
diff --git a/R/sim_mvgam.R b/R/sim_mvgam.R
index 814ca9d6..4a9f70b6 100644
--- a/R/sim_mvgam.R
+++ b/R/sim_mvgam.R
@@ -141,6 +141,7 @@ sim_mvgam = function(T = 100,
     cat('trend_rhos = ', trend_rhos, '\n\n')
     # Generate latent GP trends
     trends <- do.call(cbind, lapply(seq_len(n_lv), function(lv){
+      set.seed(runif(1, 1, 100))
       sim_exp_gp(N = T, alpha = trend_alphas[lv], rho = trend_rhos[lv])
     }))
 
diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R
index a29e6b3e..c6e09d58 100644
--- a/R/summary.mvgam.R
+++ b/R/summary.mvgam.R
@@ -15,12 +15,6 @@ summary.mvgam = function(object){
   # calculate effective degrees of freedom for smooth terms
   jam = object$jam_model
 
-  # Convert stanfit objects to coda samples
-  if(class(object$model_output) == 'stanfit'){
-    object$model_output <- coda::mcmc.list(lapply(1:NCOL(object$model_output),
-                                                  function(x) coda::mcmc(as.array(object$model_output)[,x,])))
-  }
-
 #### Standard summary of formula and model argumements ####
 message("GAM formula:")
 print(object$call)
diff --git a/Results/sim_results.rda b/Results/sim_results.rda
deleted file mode 100644
index 1c8016bc..00000000
Binary files a/Results/sim_results.rda and /dev/null differ
diff --git a/docs/README-unnamed-chunk-10-1.png b/docs/README-unnamed-chunk-10-1.png
index c48d9b8d..bb0cf3d3 100644
Binary files a/docs/README-unnamed-chunk-10-1.png and b/docs/README-unnamed-chunk-10-1.png differ
diff --git a/docs/README-unnamed-chunk-11-1.png b/docs/README-unnamed-chunk-11-1.png
index 00d14d83..59371081 100644
Binary files a/docs/README-unnamed-chunk-11-1.png and b/docs/README-unnamed-chunk-11-1.png differ
diff --git a/docs/README-unnamed-chunk-13-1.png b/docs/README-unnamed-chunk-13-1.png
index b27a9d12..24f5bf4e 100644
Binary files a/docs/README-unnamed-chunk-13-1.png and b/docs/README-unnamed-chunk-13-1.png differ
diff --git a/docs/README-unnamed-chunk-14-1.png b/docs/README-unnamed-chunk-14-1.png
index 0dbaaf2e..f0604c21 100644
Binary files a/docs/README-unnamed-chunk-14-1.png and b/docs/README-unnamed-chunk-14-1.png differ
diff --git a/docs/README-unnamed-chunk-15-1.png b/docs/README-unnamed-chunk-15-1.png
index 583fb601..62539648 100644
Binary files a/docs/README-unnamed-chunk-15-1.png and b/docs/README-unnamed-chunk-15-1.png differ
diff --git a/docs/README-unnamed-chunk-16-1.png b/docs/README-unnamed-chunk-16-1.png
index a19664c7..c2a467ae 100644
Binary files a/docs/README-unnamed-chunk-16-1.png and b/docs/README-unnamed-chunk-16-1.png differ
diff --git a/docs/README-unnamed-chunk-17-1.png b/docs/README-unnamed-chunk-17-1.png
index 095214af..da948d03 100644
Binary files a/docs/README-unnamed-chunk-17-1.png and b/docs/README-unnamed-chunk-17-1.png differ
diff --git a/docs/README-unnamed-chunk-18-1.png b/docs/README-unnamed-chunk-18-1.png
index 0ceb3e3e..e3734269 100644
Binary files a/docs/README-unnamed-chunk-18-1.png and b/docs/README-unnamed-chunk-18-1.png differ
diff --git a/docs/README-unnamed-chunk-19-1.png b/docs/README-unnamed-chunk-19-1.png
index 6271306e..162e61ac 100644
Binary files a/docs/README-unnamed-chunk-19-1.png and b/docs/README-unnamed-chunk-19-1.png differ
diff --git a/docs/README-unnamed-chunk-20-1.png b/docs/README-unnamed-chunk-20-1.png
index fe3701b4..fde5aa6b 100644
Binary files a/docs/README-unnamed-chunk-20-1.png and b/docs/README-unnamed-chunk-20-1.png differ
diff --git a/docs/README-unnamed-chunk-21-1.png b/docs/README-unnamed-chunk-21-1.png
index 4b5227a0..4a014027 100644
Binary files a/docs/README-unnamed-chunk-21-1.png and b/docs/README-unnamed-chunk-21-1.png differ
diff --git a/docs/README-unnamed-chunk-22-1.png b/docs/README-unnamed-chunk-22-1.png
index 7d84b99f..0c07348a 100644
Binary files a/docs/README-unnamed-chunk-22-1.png and b/docs/README-unnamed-chunk-22-1.png differ
diff --git a/docs/README-unnamed-chunk-23-1.png b/docs/README-unnamed-chunk-23-1.png
index be91fd92..1c78e949 100644
Binary files a/docs/README-unnamed-chunk-23-1.png and b/docs/README-unnamed-chunk-23-1.png differ
diff --git a/docs/README-unnamed-chunk-24-1.png b/docs/README-unnamed-chunk-24-1.png
index fa8e5c0a..a63fcd8b 100644
Binary files a/docs/README-unnamed-chunk-24-1.png and b/docs/README-unnamed-chunk-24-1.png differ
diff --git a/docs/README-unnamed-chunk-7-1.png b/docs/README-unnamed-chunk-7-1.png
index 3a3cc5ea..644314d0 100644
Binary files a/docs/README-unnamed-chunk-7-1.png and b/docs/README-unnamed-chunk-7-1.png differ
diff --git a/docs/README-unnamed-chunk-8-1.png b/docs/README-unnamed-chunk-8-1.png
index 759fbe87..fe6c60ad 100644
Binary files a/docs/README-unnamed-chunk-8-1.png and b/docs/README-unnamed-chunk-8-1.png differ
diff --git a/docs/README-unnamed-chunk-9-1.png b/docs/README-unnamed-chunk-9-1.png
index af7d500c..28772e82 100644
Binary files a/docs/README-unnamed-chunk-9-1.png and b/docs/README-unnamed-chunk-9-1.png differ
diff --git a/docs/README.Rmd b/docs/README.Rmd
index 0f5281ab..6ff1783a 100644
--- a/docs/README.Rmd
+++ b/docs/README.Rmd
@@ -54,16 +54,16 @@ lynx_train = lynx_full[1:50, ]
 lynx_test = lynx_full[51:60, ]
 ```
 
-Now fit an `mvgam` model; it fits a GAM in which a cyclic smooth function for `season` is estimated jointly with a full time series model for the errors (in this case an `AR2` process), rather than relying on smoothing splines that do not incorporate a concept of the future. We assume the outcome follows a Poisson distribution and estimate the model in `JAGS` using MCMC sampling (installation links are found [here](https://sourceforge.net/projects/mcmc-jags/files/)). Note that for some models, it is now possible to estimate parameters using Hamiltonian Monte Carlo in the `Stan` software via a call to the `rstan` package, which can be found [here](https://mc-stan.org/users/interfaces/rstan)
+Now fit an `mvgam` model; it fits a GAM in which a cyclic smooth function for `season` is estimated jointly with a full time series model for the errors (in this case an `AR3` process), rather than relying on smoothing splines that do not incorporate a concept of the future. We assume the outcome follows a Poisson distribution and estimate the model in `Stan` using MCMC sampling (installation links are found [here](https://mc-stan.org/users/interfaces/rstan)). Note that for some models, it is now possible to estimate parameters using Hamiltonian Monte Carlo in the `Stan` software via a call to the `rstan` package, which can be found [here](https://mc-stan.org/users/interfaces/rstan)
 ```{r, message=FALSE, warning=FALSE}
-lynx_mvgam <- mvjagam(data_train = lynx_train,
+lynx_mvgam <- mvgam(data_train = lynx_train,
                data_test = lynx_test,
                formula = y ~ s(season, bs = 'cc', k = 19),
                knots = list(season = c(0.5, 19.5)),
                family = 'poisson',
-               trend_model = 'AR2',
-               drift = F,
-               burnin = 10000,
+               trend_model = 'AR3',
+               use_stan = T,
+               burnin = 1000,
                chains = 4)
 ```
 
@@ -160,7 +160,7 @@ plot(lynx_mvgam, type = 'residuals')
 
 Another useful utility of `mvgam` is the ability to use rolling window forecasts to evaluate competing models that may represent different hypotheses about the series dynamics. Here we will fit a poorly specified model to showcase how this evaluation works. In this model, we ignore the cyclic pattern of seasonality and force it to be fairly non-wiggly. We also use a random walk process for the trend
 ```{r, message=FALSE, warning=FALSE}
-lynx_mvgam_poor <- mvjagam(data_train = lynx_train,
+lynx_mvgam_poor <- mvgam(data_train = lynx_train,
                data_test = lynx_test,
                formula = y ~ s(season, bs = 'gp', k = 3),
                family = 'poisson',
diff --git a/docs/README.md b/docs/README.md
index 3dc6553c..15c6d479 100644
--- a/docs/README.md
+++ b/docs/README.md
@@ -91,26 +91,31 @@ lynx_test = lynx_full[51:60, ]
 
 Now fit an `mvgam` model; it fits a GAM in which a cyclic smooth
 function for `season` is estimated jointly with a full time series model
-for the errors (in this case an `AR2` process), rather than relying on
+for the errors (in this case an `AR3` process), rather than relying on
 smoothing splines that do not incorporate a concept of the future. We
 assume the outcome follows a Poisson distribution and estimate the model
-in `JAGS` using MCMC sampling (installation links are found
-[here](https://sourceforge.net/projects/mcmc-jags/files/)). Note that
-for some models, it is now possible to estimate parameters using
-Hamiltonian Monte Carlo in the `Stan` software via a call to the `rstan`
-package, which can be found
-[here](https://mc-stan.org/users/interfaces/rstan)
+in `Stan` using MCMC sampling (installation links are found
+[here](https://mc-stan.org/users/interfaces/rstan)). Note that for some
+models, it is now possible to estimate parameters using Hamiltonian
+Monte Carlo in the `Stan` software via a call to the `rstan` package,
+which can be found [here](https://mc-stan.org/users/interfaces/rstan)
 
 ``` r
-lynx_mvgam <- mvjagam(data_train = lynx_train,
+lynx_mvgam <- mvgam(data_train = lynx_train,
                data_test = lynx_test,
                formula = y ~ s(season, bs = 'cc', k = 19),
                knots = list(season = c(0.5, 19.5)),
                family = 'poisson',
-               trend_model = 'AR2',
-               drift = F,
-               burnin = 10000,
+               trend_model = 'AR3',
+               use_stan = T,
+               burnin = 1000,
                chains = 4)
+#> [1] "n_eff / iter looks reasonable for all parameters"
+#> [1] "Rhat looks reasonable for all parameters"
+#> [1] "0 of 4000 iterations ended with a divergence (0%)"
+#> [1] "69 of 4000 iterations saturated the maximum tree depth of 10 (1.725%)"
+#> [1] "  Run again with max_treedepth set to a larger value to avoid saturation"
+#> [1] "E-FMI indicated no pathological behavior"
 ```
 
 Perform a series of posterior predictive checks to see if the model is
@@ -193,7 +198,7 @@ summary(lynx_mvgam)
 #> log
 #> 
 #> Trend model:
-#> AR2
+#> AR3
 #> 
 #> N series:
 #> 1
@@ -202,42 +207,43 @@ summary(lynx_mvgam)
 #> 50
 #> 
 #> Status:
-#> Fitted using runjags::run.jags()
+#> Fitted using rstan::stan()
 #> 
 #> GAM smooth term estimated degrees of freedom:
 #>            edf df
-#> s(season) 5214 17
+#> s(season) 9011 17
 #> 
 #> GAM coefficient (beta) estimates:
-#>                    2.5%         50%       97.5% Rhat n.eff
-#> (Intercept)   6.6278432  6.76661720  6.89681964 1.03    87
-#> s(season).1  -1.1277700 -0.68473902 -0.22962380 1.34   133
-#> s(season).2  -0.2229482  0.28944913  0.86452547 1.16   103
-#> s(season).3   0.6519581  1.13340639  1.73032941 1.06    45
-#> s(season).4   1.2603487  1.74322315  2.34746633 1.14    53
-#> s(season).5   1.4442400  1.94342290  2.70066342 1.42    59
-#> s(season).6   0.4456740  1.02894020  1.61028936 1.14    51
-#> s(season).7  -0.8433434 -0.22754522  0.28449496 1.07    94
-#> s(season).8  -1.3271281 -0.72048653 -0.14820460 1.01   110
-#> s(season).9  -1.5807970 -0.79675279 -0.04589625 1.03   141
-#> s(season).10 -1.2656095 -0.44005281  0.41195850 1.05    82
-#> s(season).11 -0.3393578  0.33089693  0.99094619 1.07    75
-#> s(season).12  0.7295668  1.34967548  1.91774349 1.19    56
-#> s(season).13  0.7150178  1.46943619  2.20605335 1.47    47
-#> s(season).14  0.2769248  1.09164305  1.76202078 1.42    52
-#> s(season).15 -0.6892971 -0.07206227  0.46319419 1.07    88
-#> s(season).16 -1.3544063 -0.81086047 -0.26011986 1.02   144
-#> s(season).17 -1.5709996 -1.07117890 -0.54800861 1.09   170
+#>                     2.5%         50%       97.5% Rhat n.eff
+#> (Intercept)   6.79095318  6.80222670  6.81317305    1  6082
+#> s(season).1  -1.22012054 -0.70687316 -0.13550066    1  1533
+#> s(season).2  -0.34422058  0.23657666  0.84322640    1  1784
+#> s(season).3   0.44107251  1.10363234  1.72974673    1  1725
+#> s(season).4   0.95960079  1.74495643  2.40358999    1  1310
+#> s(season).5   1.15835797  2.00171701  2.70516103    1  1244
+#> s(season).6   0.42168311  1.17645507  1.82805556    1  1499
+#> s(season).7  -0.78622518 -0.11099717  0.58368607    1  1758
+#> s(season).8  -1.35633593 -0.67924442  0.05796520    1  1762
+#> s(season).9  -1.58002424 -0.81464343  0.07228273    1  1510
+#> s(season).10 -1.22512971 -0.45385323  0.48789845    1  1448
+#> s(season).11 -0.52903561  0.27259212  1.13203943    1  1948
+#> s(season).12  0.29771680  1.22695736  2.01184910    1  1826
+#> s(season).13  0.31480665  1.39836005  2.18521162    1  1223
+#> s(season).14  0.05866046  1.09845407  1.83787390    1  1041
+#> s(season).15 -0.77257545 -0.04440334  0.54493012    1  1373
+#> s(season).16 -1.37554751 -0.78647085 -0.21917617    1  2239
+#> s(season).17 -1.56201026 -1.06732513 -0.49428374    1  1756
 #> 
 #> GAM smoothing parameter (rho) estimates:
-#>              2.5%     50%    97.5% Rhat n.eff
-#> s(season) 3.08043 3.88746 4.563715 1.01  1930
+#>               2.5%      50%   97.5% Rhat n.eff
+#> s(season) 3.287832 4.159986 4.87422    1  2857
 #> 
 #> Latent trend parameter estimates:
-#>             2.5%        50%     97.5% Rhat n.eff
-#> ar1    0.4236962  0.7313141 1.0300245 1.03  1112
-#> ar2   -0.4069726 -0.1131876 0.1877084 1.03  1213
-#> sigma  0.3680085  0.4601895 0.5877884 1.00  1010
+#>                2.5%         50%     97.5% Rhat n.eff
+#> ar1[1]    0.5231943  0.82710459 0.9895406    1  2007
+#> ar2[1]   -0.5841628 -0.22174449 0.1588171    1  2723
+#> ar3[1]   -0.3304613  0.07239184 0.4249637    1  1279
+#> sigma[1]  0.3677681  0.46126809 0.5929270    1  2374
 #> 
 ```
 
@@ -294,6 +300,9 @@ the entire series (testing and training)
 
 ``` r
 plot(lynx_mvgam, type = 'forecast', data_test = lynx_test)
+#> Out of sample DRPS:
+#> [1] 683.5812
+#> 
 ```
 
 <img src="README-unnamed-chunk-17-1.png" style="display: block; margin: auto;" />
@@ -376,7 +385,7 @@ ignore the cyclic pattern of seasonality and force it to be fairly
 non-wiggly. We also use a random walk process for the trend
 
 ``` r
-lynx_mvgam_poor <- mvjagam(data_train = lynx_train,
+lynx_mvgam_poor <- mvgam(data_train = lynx_train,
                data_test = lynx_test,
                formula = y ~ s(season, bs = 'gp', k = 3),
                family = 'poisson',
@@ -404,11 +413,11 @@ markedly better (far lower DRPS) for this evaluation timepoint
 
 ``` r
 summary(mod1_eval$series1$drps)
-#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-#>   0.3335   8.4668  96.5577  94.4811 130.9125 232.9926
+#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+#>   1.781  11.172 106.428  97.483 136.918 227.490
 summary(mod2_eval$series1$drps)
 #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-#>   40.37   46.08  300.00  282.02  436.08  658.25
+#>   39.26   48.04  299.49  284.69  446.13  672.41
 ```
 
 Nominal coverages for both models’ 90% prediction intervals
diff --git a/man/mvgam.Rd b/man/mvgam.Rd
index 9315480f..4aa0f03d 100644
--- a/man/mvgam.Rd
+++ b/man/mvgam.Rd
@@ -182,14 +182,9 @@ The \code{p} parameter is therefore fixed at \code{1.5} (i.e. a so-called Geomet
 \emph{Factor regularisation}: When using a dynamic factor model for the trends, factor precisions are given
 regularized penalty priors to theoretically allow some factors to be dropped from the model by squeezing increasing
 factors' variances to zero. This is done to help protect against selecting too many latent factors than are needed to
-capture dependencies in the data, so it can often be advantageous to set `n_lv`` to a slightly larger number. However
-larger numbers of factors do come with additional computational costs so these should be balanced as well.
-\cr
-\cr
-\emph{Residuals}: For each series, randomized quantile (i.e. Dunn-Smyth) residuals are calculated for inspecting model diagnostics
-If the fitted model is appropriate then Dunn-Smyth residuals will be standard normal in distribution and no
-autocorrelation will be evident. When a particular observation is missing, the residual is calculated by comparing independent
-draws from the model's posterior distribution
+capture dependencies in the data, so it can often be advantageous to set \verb{n_lv`` to a slightly larger number. However larger numbers of factors do come with additional computational costs so these should be balanced as well. \\cr \\cr *Residuals*: For each series, randomized quantile (i.e. Dunn-Smyth) residuals are calculated for inspecting model diagnostics If the fitted model is appropriate then Dunn-Smyth residuals will be standard normal in distribution and no autocorrelation will be evident. When a particular observation is missing, the residual is calculated by comparing independent draws from the model's posterior distribution \\cr \\cr *Using Stan*: An in-development feature of }mvgam\verb{is the ability to use Hamiltonian Monte Carlo for parameter estimation via the software}Stan\verb{(using the}rstan\verb{interface). Note that the}rstan\verb{library is currently required for this option to work, though support for other}Stan\verb{interfaces will be added in future. Also note that currently there is no support for fitting}Tweedie\verb{responses or dynamic factor models in}Stan\verb{, though again these will be added in future. However there are great advantages when using }Stan\verb{, which includes the option to estimate smooth latent trends via [Hilbert space approximate Gaussian Processes](https://arxiv.org/abs/2004.11408). This often makes sense for ecological series, which we expect to change smoothly. In }mvgam`, latent squared exponential GP trends are approximated using
+by default \code{40} basis functions, which saves computational costs compared to fitting full GPs while adequately estimating
+GP \code{alpha} and \code{rho} parameters
 }
 \seealso{
 \code{\link[rjags]{jags.model}}, \code{\link[mcgv]{jagam}}, \code{\link[mcgv]{gam}}
diff --git a/man/pfilter_mvgam_online.Rd b/man/pfilter_mvgam_online.Rd
index e590584e..74e7184e 100644
--- a/man/pfilter_mvgam_online.Rd
+++ b/man/pfilter_mvgam_online.Rd
@@ -31,7 +31,9 @@ the specified \code{threshold}. Default for this option is \code{FALSE}, relying
 to maintain particle diversity}
 
 \item{kernel_lambda}{\code{proportional numeric} specifying the strength of kernel smoothing to use when
-pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}}
+pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}. Note
+however that kernel smoothing is not currently supported for Gaussian Process trend models. For these models,
+it is still preferable to use resampling with a reasonable \code{threshold} value (i.e. \code{0.5})}
 
 \item{n_cores}{\code{integer} specifying number of cores for generating particle forecasts in parallel}
 }
diff --git a/man/pfilter_mvgam_smooth.Rd b/man/pfilter_mvgam_smooth.Rd
index a95dc0c3..01d4418a 100644
--- a/man/pfilter_mvgam_smooth.Rd
+++ b/man/pfilter_mvgam_smooth.Rd
@@ -36,12 +36,14 @@ Should be between \code{0} and \code{1}}
 the specified \code{threshold}. Note that resampling can result in loss of the original model's diversity of
 GAM beta coefficients, which may have undesirable consequences for the forecast distribution. If
 \code{use_resampling} is \code{TRUE}, some effort is made to remedy this by assigning randomly sampled draws of
-beta coefficients from the original model's distribution to each particle. This does not however guarantee that there
-will be no loss of diversity, especially if successive resamples take place. Default for this option is therefore
+GAM beta coefficients from the original model's distribution to each particle. This does not however guarantee there
+will be no loss of diversity, especially if successive resampling take place. Default for this option is therefore
 \code{FALSE}}
 
 \item{kernel_lambda}{\code{proportional numeric} specifying the strength of smoothing to use when
-pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}}
+pulling low weight particles toward the high likelihood state space. Should be between \code{0} and \code{1}. Note
+however that kernel smoothing is not currently supported for Gaussian Process trend models. For these models,
+it is still preferable to use resampling with a reasonable \code{threshold} value (i.e. \code{0.5})}
 
 \item{file_path}{\code{character} string specifying the file path for locating the particles}
 }