diff --git a/R/add_MACor.R b/R/add_MACor.R
index 7899567d..984db65d 100644
--- a/R/add_MACor.R
+++ b/R/add_MACor.R
@@ -11,7 +11,7 @@ add_MaCor = function(model_file,
                      trend_model = 'VAR1',
                      drift = FALSE){
 
-  if(trend_model %in% c('RW', 'AR1', 'AR2')){
+  if(trend_model %in% c('RW', 'AR1', 'AR2', 'AR3')){
 
     # Update transformed data
       if(any(grepl('vector<lower=0>[n_lv] sigma;',
@@ -370,6 +370,60 @@ add_MaCor = function(model_file,
                    '}\n')
         }
 
+        if(trend_model == 'AR3'){
+          model_file[max(grep('= b_raw[',
+                              model_file, fixed = TRUE))] <-
+            paste0(model_file[max(grep('= b_raw[',
+                                       model_file, fixed = TRUE))],
+                   '\n// derived latent states\n',
+                   'trend_raw[1] = ',
+                   if(drift){ 'drift + '} else {NULL},
+                   'error[1];\n',
+                   if(add_ma){
+                     paste0('epsilon[1] = error[1];\n',
+                            'epsilon[2] = theta * error[1];\n',
+                            'epsilon[3] = theta * error[2];\n')
+                   } else {
+                     NULL
+                   },
+                   'trend_raw[2] = ',
+                   if(drift){ 'drift + '} else {NULL},
+                   'ar1 .* trend_raw[1] + ',
+                   if(add_ma){
+                     'epsilon[2] + error[2];\n'
+                   } else {
+                     'error[2];\n'
+                   },
+                   'trend_raw[3] = ',
+                   if(drift){ 'drift + '} else {NULL},
+                   'ar1 .* trend_raw[2] + ',
+                   'ar2 .* trend_raw[1] + ',
+                   if(add_ma){
+                     'epsilon[3] + error[3];\n'
+                   } else {
+                     'error[3];\n'
+                   },
+                   'for (i in 4:n) {\n',
+                   if(add_ma){
+                     paste0('// lagged error ma process\n',
+                            'epsilon[i] = theta * error[i - 1];\n',
+                            '// full ARMA process\n')
+                   } else {
+                     '// full AR process\n'
+                   },
+                   'trend_raw[i] = ',
+                   if(drift){ 'drift + '} else {NULL},
+                   'ar1 .* trend_raw[i - 1] + ',
+                   'ar2 .* trend_raw[i - 2] + ',
+                   'ar3 .* trend_raw[i - 3] + ',
+                   if(add_ma){
+                     'epsilon[i] + error[i];\n'
+                   } else {
+                     'error[i];\n'
+                   },
+                   '}\n')
+        }
+
       } else {
         if(trend_model %in% c('AR1', 'RW')){
           model_file[max(grep('= b_raw[',
@@ -419,6 +473,40 @@ add_MaCor = function(model_file,
                    'epsilon[i, j] + error[i, j];\n',
                    '}\n}')
         }
+
+        if(trend_model == 'AR3'){
+          model_file[max(grep('= b_raw[',
+                              model_file, fixed = TRUE))] <-
+            paste0(model_file[max(grep('= b_raw[',
+                                       model_file, fixed = TRUE))],
+                   'for(j in 1:n_series){\n',
+                   'trend[1, j] = ',
+                   if(drift){ 'drift[j] + '} else {NULL},
+                   'error[1, j];\n',
+                   'epsilon[1, j] = error[1, j];\n',
+                   'epsilon[2, j] = theta[j] * error[1, j];\n',
+                   'epsilon[3, j] = theta[j] * error[2, j];\n',
+                   'trend[2, j] = ',
+                   if(drift){ 'drift[j] + '} else {NULL},
+                   'ar1[j] * trend[1, j] + ',
+                   'epsilon[2, j] + error[2, j];\n',
+                   'trend[3, j] = ',
+                   if(drift){ 'drift[j] + '} else {NULL},
+                   'ar1[j] * trend[2, j] + ',
+                   'ar2[j] * trend[1, j] + ',
+                   'epsilon[2, j] + error[2, j];\n',
+                   'for(i in 4:n){\n',
+                   '// lagged error ma process\n',
+                   'epsilon[i, j] = theta[j] * error[i-1, j];\n',
+                   '// full ARMA process\n',
+                   'trend[i, j] = ',
+                   if(drift){ 'drift[j] + '} else {NULL},
+                   'ar1[j] * trend[i - 1, j] + ',
+                   'ar2[j] * trend[i - 2, j] + ',
+                   'ar3[j] * trend[i - 3, j] + ',
+                   'epsilon[i, j] + error[i, j];\n',
+                   '}\n}')
+        }
       }
 
       model_file <- readLines(textConnection(model_file), n = -1)
@@ -469,10 +557,10 @@ add_MaCor = function(model_file,
                    paste0('for(i in 1:n_lv){\n',
                           'for(j in 1:n_lv){\n',
                           'if (i != j)\n',
-                          'theta[i, j] ~ std_normal();\n',
+                          'theta[i, j] ~ normal(0, 0.2);\n',
                           '}\n}')
                  } else {
-                   'theta ~ std_normal();'
+                   'theta ~ normal(0, 0.2);'
                  })
         } else {
           NULL
@@ -505,10 +593,10 @@ add_MaCor = function(model_file,
                    paste0('for(i in 1:n_series){\n',
                           'for(j in 1:n_series){\n',
                           'if (i != j)\n',
-                          'theta[i, j] ~ std_normal();\n',
+                          'theta[i, j] ~ normal(0, 0.2);\n',
                           '}\n}')
                  } else {
-                   'theta ~ std_normal();'
+                   'theta ~ normal(0, 0.2);'
                  })
         } else {
           NULL
diff --git a/R/gp.R b/R/gp.R
index c629f0ac..7cfebf79 100644
--- a/R/gp.R
+++ b/R/gp.R
@@ -136,8 +136,13 @@ make_gp_additions = function(gp_details, data,
 
   # Add coefficient indices to attribute table and to Stan data
   for(covariate in seq_along(gp_att_table)){
-    coef_indices <- grep(gp_att_table[[covariate]]$name,
-                         names(coef(mgcv_model)), fixed = TRUE)
+    # coef_indices <- grep(gp_att_table[[covariate]]$name,
+    #                      names(coef(mgcv_model)), fixed = TRUE)
+    coef_indices <- which(grepl(gp_att_table[[covariate]]$name,
+                                names(coef(mgcv_model)), fixed = TRUE) &
+                            !grepl(paste0(gp_att_table[[covariate]]$name,':'),
+                                   names(coef(mgcv_model)), fixed = TRUE) == TRUE)
+
     gp_att_table[[covariate]]$first_coef <- min(coef_indices)
     gp_att_table[[covariate]]$last_coef <- max(coef_indices)
 
diff --git a/R/plot_mvgam_fc.R b/R/plot_mvgam_fc.R
index 15317bc9..cc46d0cd 100644
--- a/R/plot_mvgam_fc.R
+++ b/R/plot_mvgam_fc.R
@@ -248,12 +248,18 @@ plot_mvgam_fc = function(object, series = 1, newdata, data_test,
       dplyr::arrange(time) %>%
       dplyr::pull(y)
 
-    if(tolower(object$family) %in% c('beta', 'lognormal', 'gamma')){
+    if(tolower(object$family) %in% c('beta')){
       ylim <- c(min(cred, min(ytrain, na.rm = TRUE)),
                 max(cred, max(ytrain, na.rm = TRUE)))
       ymin <- max(0, ylim[1])
       ymax <- min(1, ylim[2])
       ylim <- c(ymin, ymax)
+    } else if(tolower(object$family) %in% c('lognormal', 'gamma')){
+      ylim <- c(min(cred, min(ytrain, na.rm = TRUE)),
+                max(cred, max(ytrain, na.rm = TRUE)))
+      ymin <- max(0, ylim[1])
+      ymax <- max(ylim)
+      ylim <- c(ymin, ymax)
     } else {
       ylim <- c(min(cred, min(ytrain, na.rm = TRUE)),
                 max(cred, max(ytrain, na.rm = TRUE)))
diff --git a/R/sim_mvgam.R b/R/sim_mvgam.R
index 028f1863..72cd7311 100644
--- a/R/sim_mvgam.R
+++ b/R/sim_mvgam.R
@@ -47,7 +47,7 @@
 #'be replicated `n_series` times. Defaults to small random values between `-0.5` and `0.5` on the link scale
 #'@param prop_missing \code{numeric} stating proportion of observations that are missing. Should be between
 #'\code{0} and \code{0.8}, inclusive
-#'@param prop_train \code{numeric} stating the proportion of data to use for training. Should be between \code{0.25} and \code{0.75}
+#'@param prop_train \code{numeric} stating the proportion of data to use for training. Should be between \code{0.2} and \code{1}
 #'@return A \code{list} object containing outputs needed for \code{\link{mvgam}}, including 'data_train' and 'data_test',
 #'as well as some additional information about the simulated seasonality and trend dependencies
 #'@examples
diff --git a/src/RcppExports.o b/src/RcppExports.o
index 3738bd53..849c29a4 100644
Binary files a/src/RcppExports.o and b/src/RcppExports.o differ
diff --git a/src/mvgam.dll b/src/mvgam.dll
index 8f336fda..0d00c5a1 100644
Binary files a/src/mvgam.dll and b/src/mvgam.dll differ
diff --git a/src/trend_funs.o b/src/trend_funs.o
index 371ea49f..faf098a5 100644
Binary files a/src/trend_funs.o and b/src/trend_funs.o differ
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index 072501ea..d4d55f4e 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ