diff --git a/DESCRIPTION b/DESCRIPTION
index e0cea7c7..c3848c2f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: mvgam
 Title: Multivariate (Dynamic) Generalized Additive Models
-Version: 1.0.7
+Version: 1.0.8
 Date: 2023-11-01
 Authors@R:
     person("Nicholas J", "Clark", , "nicholas.j.clark1214@gmail.com", role = c("aut", "cre"),
diff --git a/R/forecast.mvgam.R b/R/forecast.mvgam.R
index c9da361c..9054fb7d 100644
--- a/R/forecast.mvgam.R
+++ b/R/forecast.mvgam.R
@@ -744,6 +744,14 @@ forecast_draws = function(object,
                             time = data_test$time,
                             cap = suppressWarnings(linkfun(data_test$cap,
                                                            link = family$link)))
+
+          if(any(is.na(cap$cap)) | any(is.infinite(cap$cap))){
+            stop(paste0('Missing or infinite values found for some "cap" terms\n',
+                        'after transforming to the ',
+                        family$link, ' link scale'),
+                 call. = FALSE)
+          }
+
         } else {
           cap <- NULL
         }
diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R
index 3937bf6b..0e287842 100644
--- a/R/get_mvgam_priors.R
+++ b/R/get_mvgam_priors.R
@@ -732,27 +732,27 @@ get_mvgam_priors = function(formula,
 
     if(trend_model %in% c('PWlinear', 'PWlogistic')){
       # Need to fix this as a next priority
-      trend_df <- NULL
-      # trend_df <- data.frame(param_name = c('vector[n_series] k_trend;',
-      #                                       'vector[n_series] m_trend;'),
-      #                        param_length = length(unique(data_train$series)),
-      #                        param_info = c('base trend growth rates',
-      #                                       'trend offset parameters'),
-      #                        prior = c('k_trend ~ std_normal();',
-      #                                  'm_trend ~ student_t(3, 0, 2.5);'),
-      #                        example_change = c(paste0(
-      #                          'k ~ normal(',
-      #                          round(runif(min = -1, max = 1, n = 1), 2),
-      #                          ', ',
-      #                          round(runif(min = 0.1, max = 1, n = 1), 2),
-      #                          ');'),
-      #                          paste0(
-      #                            'm ~ normal(',
-      #                            round(runif(min = -1, max = 1, n = 1), 2),
-      #                            ', ',
-      #                            round(runif(min = 0.1, max = 1, n = 1), 2),
-      #                            ');')))
-      #
+      # trend_df <- NULL
+      trend_df <- data.frame(param_name = c('vector[n_series] k_trend;',
+                                            'vector[n_series] m_trend;'),
+                             param_length = length(unique(data_train$series)),
+                             param_info = c('base trend growth rates',
+                                            'trend offset parameters'),
+                             prior = c('k_trend ~ std_normal();',
+                                       'm_trend ~ student_t(3, 0, 2.5);'),
+                             example_change = c(paste0(
+                               'k ~ normal(',
+                               round(runif(min = -1, max = 1, n = 1), 2),
+                               ', ',
+                               round(runif(min = 0.1, max = 1, n = 1), 2),
+                               ');'),
+                               paste0(
+                                 'm ~ normal(',
+                                 round(runif(min = -1, max = 1, n = 1), 2),
+                                 ', ',
+                                 round(runif(min = 0.1, max = 1, n = 1), 2),
+                                 ');')))
+
 
     }
     if(trend_model == 'GP'){
diff --git a/R/index-mvgam.R b/R/index-mvgam.R
index 8e5e34c0..369d4d11 100644
--- a/R/index-mvgam.R
+++ b/R/index-mvgam.R
@@ -138,7 +138,7 @@ variables.mvgam = function(x, ...){
                        'ar1', 'ar2',
                        'ar3', 'A',
                        'Sigma', 'error', 'theta',
-                       'k_trend', 'delta', 'm_trend'), collapse = '|'),
+                       'k_trend', 'delta_trend', 'm_trend'), collapse = '|'),
                parnames) &
          !grepl('sigma_obs', parnames, fixed = TRUE) &
          !grepl('sigma_raw', parnames, fixed = TRUE))){
@@ -147,7 +147,7 @@ variables.mvgam = function(x, ...){
                                 'ar1', 'ar2',
                                 'ar3', 'A',
                                 'Sigma', 'error', 'theta',
-                                'k_trend', 'delta', 'm_trend'), collapse = '|'),
+                                'k_trend', 'delta_trend', 'm_trend'), collapse = '|'),
                         parnames) &
       !grepl('sigma_obs', parnames, fixed = TRUE) &
       !grepl('sigma_raw', parnames, fixed = TRUE)
diff --git a/R/mvgam.R b/R/mvgam.R
index 745e28df..fc9c3d34 100644
--- a/R/mvgam.R
+++ b/R/mvgam.R
@@ -78,10 +78,11 @@
 #'   \item `'AR2'` or `AR(p = 2)`
 #'   \item `'AR3'` or `AR(p = 3)`
 #'   \item `'VAR1'`  or `VAR()`(only available in \code{Stan})
+#'   \item `'PWlogistic`, `'PWlinear` or `PW()` (only available in \code{Stan})
 #'   \item `'GP'` (Gaussian Process with squared exponential kernel;
 #'only available in \code{Stan})}
 #'
-#'For all types apart from `'GP'`, moving average and/or correlated
+#'For all types apart from `'GP'` and `PW()`, moving average and/or correlated
 #'process error terms can also be estimated (for example, `RW(cor = TRUE)` will set up a
 #'multivariate Random Walk if `n_series > 1`). See [mvgam_trends] for more details
 #'@param trend_map Optional `data.frame` specifying which series should depend on which latent
@@ -170,7 +171,12 @@
 #'\cr
 #'\cr
 #'*Formula syntax*: Details of the formula syntax used by \pkg{mvgam} can be found in
-#'\code{\link[mgcv]{formula.gam}}.
+#'\code{\link[mgcv]{formula.gam}}. Note that it is possible to supply an empty formula where
+#'there are no predictors or intercepts in the observation model (i.e. `y ~ 0` or `y ~ -1`).
+#'In this case, an intercept-only observation model will be set up but the intercept coefficient
+#'will be fixed at zero. This can be handy if you wish to fit pure State-Space models where
+#'the variation in the dynamic trend controls the average expectation, and/or where intercepts
+#'are non-identifiable (as in piecewise trends, see examples below)
 #'\cr
 #'\cr
 #'*Families and link functions*: Details of families supported by \pkg{mvgam} can be found in
@@ -354,7 +360,7 @@
 #'                        trend = c(1,1,2))
 #'
 #' # Fit the model using AR1 trends
-#' mod1 <- mvgam(y ~ s(season, bs = 'cc'),
+#' mod <- mvgam(y ~ s(season, bs = 'cc'),
 #'               trend_map = trend_map,
 #'               trend_model = 'AR1',
 #'               data = mod_data,
@@ -362,12 +368,12 @@
 #'
 #' # The mapping matrix is now supplied as data to the model in the 'Z' element
 #' mod1$model_data$Z
-#' code(mod1)
+#' code(mod)
 #'
 #' # The first two series share an identical latent trend; the third is different
-#' plot(mod1, type = 'trend', series = 1)
-#' plot(mod1, type = 'trend', series = 2)
-#' plot(mod1, type = 'trend', series = 3)
+#' plot(mod, type = 'trend', series = 1)
+#' plot(mod, type = 'trend', series = 2)
+#' plot(mod, type = 'trend', series = 3)
 #'
 #' # Example of how to use dynamic coefficients
 #' # Simulate a time-varying coefficient for the effect of temperature
@@ -418,7 +424,7 @@
 #'
 #' # Fit a model that includes the offset in the linear predictor as well as
 #' # hierarchical seasonal smooths
-#' mod1 <- mvgam(formula = y ~ offset(offset) +
+#' mod <- mvgam(formula = y ~ offset(offset) +
 #'               s(series, bs = 're') +
 #'               s(season, bs = 'cc') +
 #'               s(season, by = series, m = 1, k = 5),
@@ -428,25 +434,25 @@
 #'
 #' # Inspect the model file to see the modification to the linear predictor
 #' # (eta)
-#' mod1$model_file
+#' mod$model_file
 #'
 #' # Forecasts for the first two series will differ in magnitude
 #' layout(matrix(1:2, ncol = 2))
-#' plot(mod1, type = 'forecast', series = 1, newdata = dat$data_test,
+#' plot(mod, type = 'forecast', series = 1, newdata = dat$data_test,
 #'      ylim = c(0, 75))
-#' plot(mod1, type = 'forecast', series = 2, newdata = dat$data_test,
+#' plot(mod, type = 'forecast', series = 2, newdata = dat$data_test,
 #'      ylim = c(0, 75))
 #' layout(1)
 #'
 #' # Changing the offset for the testing data should lead to changes in
 #' # the forecast
 #' dat$data_test$offset <- dat$data_test$offset - 2
-#' plot(mod1, 'forecast', newdata = dat$data_test)
+#' plot(mod, 'forecast', newdata = dat$data_test)
 #'
 #' # Relative Risks can be computed by fixing the offset to the same value
 #' # for each series
 #' dat$data_test$offset <- rep(1, NROW(dat$data_test))
-#' preds_rr <- predict(mod1, type = 'link', newdata = dat$data_test)
+#' preds_rr <- predict(mod, type = 'link', newdata = dat$data_test)
 #' series1_inds <- which(dat$data_test$series == 'series_1')
 #' series2_inds <- which(dat$data_test$series == 'series_2')
 #'
@@ -467,6 +473,66 @@
 #'  }
 #' layout(1)
 #'
+#' # Example of logistic growth with possible changepoints
+#' # Simple logistic growth model
+#' dNt = function(r, N, k){
+#'    r * N * (k - N)
+#' }
+#'
+#' # Iterate growth through time
+#' Nt = function(r, N, t, k) {
+#' for (i in 1:(t - 1)) {
+#'
+#'  # population at next time step is current population + growth,
+#'  # but we introduce several 'shocks' as changepoints
+#'  if(i %in% c(5, 15, 25, 41, 45, 60, 80)){
+#'    N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1),
+#'                                  N[i], k))
+#'    } else {
+#'    N[i + 1] <- max(1, N[i] + dNt(r, N[i], k))
+#'    }
+#'   }
+#'  N
+#' }
+#'
+#' # Simulate expected values
+#' set.seed(11)
+#' expected <- Nt(0.004, 2, 100, 30)
+#' plot(expected, xlab = 'Time')
+#'
+#' # Take Poisson draws
+#' y <- rpois(100, expected)
+#' plot(y, xlab = 'Time')
+#'
+#' # Assemble data into dataframe and model. We set a
+#' # fixed carrying capacity of 35 for this example, but note that
+#' # this value is not required to be fixed at each timepoint
+#' mod_data <- data.frame(y = y,
+#'                        time = 1:100,
+#'                        cap = 35,
+#'                        series = as.factor('series_1'))
+#' plot_mvgam_series(data = mod_data)
+#'
+#' # The intercept is nonidentifiable when using piecewise
+#' # trends because the trend functions have their own offset
+#' # parameters 'm'; it is recommended to always drop intercepts
+#' # when using these trend models
+#' mod <- mvgam(y ~ 0,
+#'              trend_model = PW(growth = 'logistic'),
+#'              family = poisson(),
+#'              data = mod_data)
+#' summary(mod)
+#'
+#' # Plot the posterior hindcast
+#' plot(mod, type = 'forecast')
+#'
+#' # View the changepoints with ggplot2 utilities
+#' library(ggplot2)
+#' mcmc_plot(mod, variable = 'delta_trend',
+#'           regex = TRUE) +
+#' scale_y_discrete(labels = mod$trend_model$changepoints) +
+#' labs(y = 'Potential changepoint',
+#'      x = 'Rate change')
 #' }
 #'@export
 
@@ -754,12 +820,7 @@ mvgam = function(formula,
                                   family = family,
                                   use_lv = use_lv,
                                   n_lv = n_lv,
-                                  trend_model =
-                                    if(trend_model %in% c('PWlinear', 'PWlogistic')){
-                                      'RW'
-                                    } else {
-                                      trend_model
-                                    },
+                                  trend_model = trend_model,
                                   trend_map = trend_map,
                                   drift = drift)
 
diff --git a/R/piecewise_trends.R b/R/piecewise_trends.R
index 84f162b3..4aad1cb9 100644
--- a/R/piecewise_trends.R
+++ b/R/piecewise_trends.R
@@ -15,11 +15,36 @@
 #'   Large values will allow many changepoints and a more flexible trend, while
 #'   small values will allow few changepoints. Default is `0.05`.
 #' @param growth Character string specifying either 'linear' or 'logistic' growth of
-#' the trend. If 'logistic', a variable labelled 'cap' MUST be in \code{data} to specify the
-#' maximum saturation point for the trend (see examples in \code{\link{mvgam}} for details).
+#' the trend. If 'logistic', a variable labelled `cap` MUST be in \code{data} to specify the
+#' maximum saturation point for the trend (see details and examples in \code{\link{mvgam}} for
+#' more information).
 #' Default is 'linear'.
 #' @return An object of class \code{mvgam_trend}, which contains a list of
 #' arguments to be interpreted by the parsing functions in \code{mvgam}
+#' @details
+#'*Offsets and intercepts*:
+#' For each of these trend models, an offset parameter is included in the trend
+#' estimation process. This parameter will be incredibly difficult to identify
+#' if you also include an intercept in the observation formula. For that reason,
+#' it is highly recommended that you drop the intercept from the formula
+#' (i.e. `y ~ x + 0` or `y ~ x - 1`, where `x` are your optional predictor terms).
+#'\cr
+#'\cr
+#'*Logistic growth and the cap variable*:
+#' When forecasting growth, there is often some maximum achievable point that
+#' a time series can reach. For example, total market size, total population size
+#' or carrying capacity in population dynamics. It can be advantageous for the forecast
+#' to saturate at or near this point so that predictions are more sensible.
+#' This function allows you to make forecasts using a logistic growth trend model,
+#' with a specified carrying capacity. Note that this capacity does not need to be static
+#' over time, it can vary with each series x timepoint combination if necessary. But you
+#' must supply a `cap` value for each observation in the data when using `growth = 'logistic'`.
+#' For observation families that use a non-identity link function, the `cap` value will
+#' be internally transformed to the link scale (i.e. your specified `cap` will be log
+#' transformed if you are using a `poisson()` or `nb()` family). It is therefore important
+#' that you specify the `cap` values on the scale of your outcome. Note also that
+#' no missing values are allowed in `cap`.
+#'
 #' @rdname piecewise_trends
 #' @export
 PW = function(n_changepoints = 10,
@@ -64,8 +89,13 @@ add_piecewise = function(model_file,
     if(!(exists('cap', where = data_train))) {
       stop('Capacities must be supplied as a variable named "cap" for logistic growth',
            call. = FALSE)
+    }
 
+    if(any(is.na(data_train$cap))){
+      stop('Missing values found for some "cap" terms',
+           call. = FALSE)
     }
+
     # Matrix of capacities per series (these must operate on the link scale)
     all_caps <- data.frame(series = as.numeric(data_train$series),
                            time = data_train$time,
@@ -73,6 +103,13 @@ add_piecewise = function(model_file,
                                                           link = family$link))) %>%
       dplyr::arrange(time, series)
 
+    if(any(is.na(all_caps$cap)) | any(is.infinite(all_caps$cap))){
+      stop(paste0('Missing or infinite values found for some "cap" terms\n',
+                  'after transforming to the ',
+                  family$link, ' link scale'),
+           call. = FALSE)
+    }
+
     if(!is.null(data_test)){
       if(!(exists('cap', where = data_test))) {
         stop('Capacities must also be supplied in "newdata" for logistic growth predictions',
@@ -85,6 +122,13 @@ add_piecewise = function(model_file,
                                    cap = suppressWarnings(linkfun(data_test$cap,
                                                                   link = family$link)))) %>%
         dplyr::arrange(time, series)
+
+      if(any(is.na(all_caps$cap)) | any(is.infinite(all_caps$cap))){
+        stop(paste0('Missing or infinite values found for some "cap" terms\n',
+                    'after transforming to the ',
+                    family$link, ' link scale'),
+             call. = FALSE)
+      }
     }
 
     cap <- matrix(NA, nrow = length(unique(all_caps$time)),
@@ -317,7 +361,7 @@ add_piecewise = function(model_file,
            "// trend offset parameters\n",
            "vector[n_series] m_trend;\n\n",
            "// trend rate adjustments per series\n",
-           "matrix[n_changepoints, n_series] delta;\n")
+           "matrix[n_changepoints, n_series] delta_trend;\n")
   model_file <- readLines(textConnection(model_file), n = -1)
 
   # Update the transformed parameters block
@@ -333,9 +377,9 @@ add_piecewise = function(model_file,
     '// trend estimates\n',
     'for (s in 1 : n_series) {\n',
     if(trend_model == 'PWlogistic'){
-      'trend[1 : n, s] = logistic_trend(k_trend[s], m_trend[s], to_vector(delta[,s]), time, to_vector(cap[,s]), A, t_change, n_changepoints);\n'
+      'trend[1 : n, s] = logistic_trend(k_trend[s], m_trend[s], to_vector(delta_trend[,s]), time, to_vector(cap[,s]), A, t_change, n_changepoints);\n'
     } else {
-      'trend[1 : n, s] = linear_trend(k_trend[s], m_trend[s], to_vector(delta[,s]), time, A, t_change);\n'
+      'trend[1 : n, s] = linear_trend(k_trend[s], m_trend[s], to_vector(delta_trend[,s]), time, A, t_change);\n'
     },
 
     '}\n')
@@ -353,13 +397,15 @@ add_piecewise = function(model_file,
     paste0('// trend parameter priors\n',
            'm_trend ~ student_t(3, 0, 2.5);\n',
            'k_trend ~ std_normal();\n',
-           'to_vector(delta) ~ double_exponential(0, changepoint_scale);\n',
+           'to_vector(delta_trend) ~ double_exponential(0, changepoint_scale);\n',
            model_file[grep("// likelihood functions", model_file, fixed = TRUE) - 1])
+  model_file <- readLines(textConnection(model_file), n = -1)
 
   # Update the generated quantities block
   model_file <- model_file[-grep("vector[n_series] tau;", model_file, fixed = TRUE)]
   tau_start <- grep("tau[s] = pow(sigma[s], -2.0);", model_file, fixed = TRUE) - 1
   model_file <- model_file[-c(tau_start:(tau_start + 2))]
+  model_file <- readLines(textConnection(model_file), n = -1)
 
   #### Return ####
   return(list(model_file = model_file,
diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R
index 67c6acc4..0fd677e9 100644
--- a/R/summary.mvgam.R
+++ b/R/summary.mvgam.R
@@ -480,6 +480,20 @@ if(object$use_lv){
 
 if(!object$use_lv){
   if(attr(object$model_data, 'trend_model') != 'None'){
+    if(attr(object$model_data, 'trend_model') %in% c('PWlinear', 'PWlogistic')){
+      cat("\nLatent trend growth rate estimates:\n")
+      print(suppressWarnings(mcmc_summary(object$model_output, c('k_trend'),
+                                          digits = digits,
+                                          variational = object$algorithm %in%
+                                            c('fullrank', 'meanfield')))[,c(3:7)])
+
+      cat("\nLatent trend offset estimates:\n")
+      print(suppressWarnings(mcmc_summary(object$model_output, c('m_trend'),
+                                          digits = digits,
+                                          variational = object$algorithm %in%
+                                            c('fullrank', 'meanfield')))[,c(3:7)])
+    }
+
     if(attr(object$model_data, 'trend_model') == 'RW'){
       if(object$drift){
         cat("\nLatent trend drift and sigma estimates:\n")
diff --git a/R/trends.R b/R/trends.R
index 63cfee29..30bb0f1d 100644
--- a/R/trends.R
+++ b/R/trends.R
@@ -8,6 +8,7 @@
 #'   \item `RW()`
 #'   \item `AR(p = 1, 2, or 3)`
 #'   \item `VAR()`(only available in \code{Stan})
+#'   \item `PW()` (piecewise linear or logistic trends; only available in \code{Stan})
 #'   \item `'GP'` (Gaussian Process with squared exponential kernel;
 #'only available in \code{Stan})}
 #'
@@ -40,6 +41,8 @@
 #'\item 'VARMA'
 #'\item 'VARMAcor'
 #'\item 'VARMA1,1cor'
+#'\item 'PWlinear'
+#'\item 'PWlogistic'
 #'\item 'GP'
 #'\item 'None'
 #'}
@@ -56,10 +59,15 @@
 #'Heaps (2022). Stationarity is not enforced when using `AR1`, `AR2` or `AR3` models,
 #'though this can be changed by the user by specifying lower and upper bounds on autoregressive
 #'parameters using functionality in [get_mvgam_priors] and the `priors` argument in
-#'[mvgam]
-#'@seealso \code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}}
+#'[mvgam]. Piecewise trends follow the formulation in the popular `prophet` package produced
+#'by `Facebook`, where users can allow for changepoints to control the potential flexibility
+#'of the trend. See Taylor and Letham (2018) for details
+#'@seealso \code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}}, \code{\link{PW}}
 #' @references Sarah E. Heaps (2022) Enforcing stationarity through the prior in Vector Autoregressions.
 #' Journal of Computational and Graphical Statistics. 32:1, 1-10.
+#'
+#' Sean J. Taylor and Benjamin Letham (2018) Forecasting at scale.
+#' The American Statistician 72.1, 37-45.
 #' @name mvgam_trends
 NULL
 
@@ -642,7 +650,7 @@ trend_par_names = function(trend_model,
     }
 
     if(trend_model %in% c('PWlinear', 'PWlogistic')){
-      param <- c('trend', 'delta', 'k_trend', 'm_trend')
+      param <- c('trend', 'delta_trend', 'k_trend', 'm_trend')
     }
 
   }
@@ -698,24 +706,24 @@ extract_trend_pars = function(object, keep_all_estimates = TRUE,
 
   # delta params for piecewise trends
   if(attr(object$model_data, 'trend_model') %in% c('PWlinear', 'PWlogistic')){
-    out$delta <- lapply(seq_along(levels(object$obs_data$series)), function(series){
+    out$delta_trend <- lapply(seq_along(levels(object$obs_data$series)), function(series){
       if(object$fit_engine == 'stan'){
-        delta_estimates <- mvgam:::mcmc_chains(object$model_output, 'delta')[,seq(series,
-                                                                          dim(mvgam:::mcmc_chains(object$model_output,
-                                                                                          'delta'))[2],
+        delta_estimates <- mvgam:::mcmc_chains(object$model_output, 'delta_trend')[,seq(series,
+                                                                          dim(mcmc_chains(object$model_output,
+                                                                                          'delta_trend'))[2],
                                                                           by = NCOL(object$ytimes))]
       } else {
-        delta_estimates <- mcmc_chains(object$model_output, 'delta')[,starts[series]:ends[series]]
+        delta_estimates <- mcmc_chains(object$model_output, 'delta_trend')[,starts[series]:ends[series]]
       }
     })
     if(attr(object$model_data, 'trend_model') == 'PWlogistic'){
       out$cap <- lapply(seq_along(levels(object$obs_data$series)), function(series){
-        t(replicate(NROW(out$delta[[1]]), object$trend_model$cap[,series]))
+        t(replicate(NROW(out$delta_trend[[1]]), object$trend_model$cap[,series]))
       })
     }
-    out$changepoints <- t(replicate(NROW(out$delta[[1]]), object$trend_model$changepoints))
-    out$change_freq <- replicate(NROW(out$delta[[1]]), object$trend_model$change_freq)
-    out$change_scale <- replicate(NROW(out$delta[[1]]), object$trend_model$changepoint_scale)
+    out$changepoints <- t(replicate(NROW(out$delta_trend[[1]]), object$trend_model$changepoints))
+    out$change_freq <- replicate(NROW(out$delta_trend[[1]]), object$trend_model$change_freq)
+    out$change_scale <- replicate(NROW(out$delta_trend[[1]]), object$trend_model$changepoint_scale)
   }
 
   # Latent trend loadings for dynamic factor models
@@ -934,10 +942,10 @@ extract_general_trend_pars = function(samp_index, trend_pars){
 
     if(names(trend_pars)[x] %in% c('last_lvs', 'lv_coefs', 'last_trends',
                                    'A', 'Sigma', 'theta', 'b_gp', 'error',
-                                   'delta', 'cap')){
+                                   'delta_trend', 'cap')){
 
       if(names(trend_pars)[x] %in% c('last_lvs', 'lv_coefs', 'last_trends',
-                                     'b_gp', 'delta', 'cap')){
+                                     'b_gp', 'delta_trend', 'cap')){
         out <- unname(lapply(trend_pars[[x]], `[`, samp_index, ))
       }
 
@@ -1381,7 +1389,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
     }
 
     if(trend_model == 'PWlinear'){
-      trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta), function(x){
+      trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta_trend), function(x){
 
         # Sample forecast horizon changepoints
         n_changes <- stats::rpois(1, (trend_pars$change_freq *
@@ -1397,7 +1405,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
         t_change_new <- sample(min(time):max(time), n_changes)
 
         # Combine with changepoints from the history
-        deltas <- c(trend_pars$delta[[x]], deltas_new)
+        deltas <- c(trend_pars$delta_trend[[x]], deltas_new)
         changepoint_ts <- c(trend_pars$changepoints, t_change_new)
 
         # Generate a trend draw
@@ -1413,7 +1421,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
     }
 
     if(trend_model == 'PWlogistic'){
-      trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta), function(x){
+      trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta_trend), function(x){
 
         # Sample forecast horizon changepoints
         n_changes <- stats::rpois(1, (trend_pars$change_freq *
@@ -1428,7 +1436,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
         t_change_new <- sample(min(time):max(time), n_changes)
 
         # Combine with changepoints from the history
-        deltas <- c(trend_pars$delta[[x]], deltas_new)
+        deltas <- c(trend_pars$delta_trend[[x]], deltas_new)
         changepoint_ts <- c(trend_pars$changepoints, t_change_new)
 
         # Get historical capacities
diff --git a/man/get_mvgam_priors.Rd b/man/get_mvgam_priors.Rd
index a23abb86..b56b8b55 100644
--- a/man/get_mvgam_priors.Rd
+++ b/man/get_mvgam_priors.Rd
@@ -78,10 +78,11 @@ and the observation process is the only source of error; similarly to what is es
 \item \code{'AR2'} or \code{AR(p = 2)}
 \item \code{'AR3'} or \code{AR(p = 3)}
 \item \code{'VAR1'}  or \code{VAR()}(only available in \code{Stan})
+\item \verb{'PWlogistic}, \verb{'PWlinear} or \code{PW()} (only available in \code{Stan})
 \item \code{'GP'} (Gaussian Process with squared exponential kernel;
 only available in \code{Stan})}
 
-For all types apart from \code{'GP'}, moving average and/or correlated
+For all types apart from \code{'GP'} and \code{PW()}, moving average and/or correlated
 process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a
 multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details}
 
diff --git a/man/mvgam.Rd b/man/mvgam.Rd
index f6735ee7..896c64bd 100644
--- a/man/mvgam.Rd
+++ b/man/mvgam.Rd
@@ -133,10 +133,11 @@ and the observation process is the only source of error; similarly to what is es
 \item \code{'AR2'} or \code{AR(p = 2)}
 \item \code{'AR3'} or \code{AR(p = 3)}
 \item \code{'VAR1'}  or \code{VAR()}(only available in \code{Stan})
+\item \verb{'PWlogistic}, \verb{'PWlinear} or \code{PW()} (only available in \code{Stan})
 \item \code{'GP'} (Gaussian Process with squared exponential kernel;
 only available in \code{Stan})}
 
-For all types apart from \code{'GP'}, moving average and/or correlated
+For all types apart from \code{'GP'} and \code{PW()}, moving average and/or correlated
 process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a
 multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details}
 
@@ -253,7 +254,12 @@ in a Bayesian framework using Markov Chain Monte Carlo.
 \cr
 \cr
 \emph{Formula syntax}: Details of the formula syntax used by \pkg{mvgam} can be found in
-\code{\link[mgcv]{formula.gam}}.
+\code{\link[mgcv]{formula.gam}}. Note that it is possible to supply an empty formula where
+there are no predictors or intercepts in the observation model (i.e. \code{y ~ 0} or \code{y ~ -1}).
+In this case, an intercept-only observation model will be set up but the intercept coefficient
+will be fixed at zero. This can be handy if you wish to fit pure State-Space models where
+the variation in the dynamic trend controls the average expectation, and/or where intercepts
+are non-identifiable (as in piecewise trends, see examples below)
 \cr
 \cr
 \emph{Families and link functions}: Details of families supported by \pkg{mvgam} can be found in
@@ -428,7 +434,7 @@ trend_map <- data.frame(series = unique(mod_data$series),
                        trend = c(1,1,2))
 
 # Fit the model using AR1 trends
-mod1 <- mvgam(y ~ s(season, bs = 'cc'),
+mod <- mvgam(y ~ s(season, bs = 'cc'),
               trend_map = trend_map,
               trend_model = 'AR1',
               data = mod_data,
@@ -436,12 +442,12 @@ mod1 <- mvgam(y ~ s(season, bs = 'cc'),
 
 # The mapping matrix is now supplied as data to the model in the 'Z' element
 mod1$model_data$Z
-code(mod1)
+code(mod)
 
 # The first two series share an identical latent trend; the third is different
-plot(mod1, type = 'trend', series = 1)
-plot(mod1, type = 'trend', series = 2)
-plot(mod1, type = 'trend', series = 3)
+plot(mod, type = 'trend', series = 1)
+plot(mod, type = 'trend', series = 2)
+plot(mod, type = 'trend', series = 3)
 
 # Example of how to use dynamic coefficients
 # Simulate a time-varying coefficient for the effect of temperature
@@ -492,7 +498,7 @@ dat$data_test$offset <- 0.5 * as.numeric(dat$data_test$series)
 
 # Fit a model that includes the offset in the linear predictor as well as
 # hierarchical seasonal smooths
-mod1 <- mvgam(formula = y ~ offset(offset) +
+mod <- mvgam(formula = y ~ offset(offset) +
               s(series, bs = 're') +
               s(season, bs = 'cc') +
               s(season, by = series, m = 1, k = 5),
@@ -502,25 +508,25 @@ mod1 <- mvgam(formula = y ~ offset(offset) +
 
 # Inspect the model file to see the modification to the linear predictor
 # (eta)
-mod1$model_file
+mod$model_file
 
 # Forecasts for the first two series will differ in magnitude
 layout(matrix(1:2, ncol = 2))
-plot(mod1, type = 'forecast', series = 1, newdata = dat$data_test,
+plot(mod, type = 'forecast', series = 1, newdata = dat$data_test,
      ylim = c(0, 75))
-plot(mod1, type = 'forecast', series = 2, newdata = dat$data_test,
+plot(mod, type = 'forecast', series = 2, newdata = dat$data_test,
      ylim = c(0, 75))
 layout(1)
 
 # Changing the offset for the testing data should lead to changes in
 # the forecast
 dat$data_test$offset <- dat$data_test$offset - 2
-plot(mod1, 'forecast', newdata = dat$data_test)
+plot(mod, 'forecast', newdata = dat$data_test)
 
 # Relative Risks can be computed by fixing the offset to the same value
 # for each series
 dat$data_test$offset <- rep(1, NROW(dat$data_test))
-preds_rr <- predict(mod1, type = 'link', newdata = dat$data_test)
+preds_rr <- predict(mod, type = 'link', newdata = dat$data_test)
 series1_inds <- which(dat$data_test$series == 'series_1')
 series2_inds <- which(dat$data_test$series == 'series_2')
 
@@ -541,6 +547,66 @@ for(i in 2:50){
  }
 layout(1)
 
+# Example of logistic growth with possible changepoints
+# Simple logistic growth model
+dNt = function(r, N, k){
+   r * N * (k - N)
+}
+
+# Iterate growth through time
+Nt = function(r, N, t, k) {
+for (i in 1:(t - 1)) {
+
+ # population at next time step is current population + growth,
+ # but we introduce several 'shocks' as changepoints
+ if(i \%in\% c(5, 15, 25, 41, 45, 60, 80)){
+   N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1),
+                                 N[i], k))
+   } else {
+   N[i + 1] <- max(1, N[i] + dNt(r, N[i], k))
+   }
+  }
+ N
+}
+
+# Simulate expected values
+set.seed(11)
+expected <- Nt(0.004, 2, 100, 30)
+plot(expected, xlab = 'Time')
+
+# Take Poisson draws
+y <- rpois(100, expected)
+plot(y, xlab = 'Time')
+
+# Assemble data into dataframe and model. We set a
+# fixed carrying capacity of 35 for this example, but note that
+# this value is not required to be fixed at each timepoint
+mod_data <- data.frame(y = y,
+                       time = 1:100,
+                       cap = 35,
+                       series = as.factor('series_1'))
+plot_mvgam_series(data = mod_data)
+
+# The intercept is nonidentifiable when using piecewise
+# trends because the trend functions have their own offset
+# parameters 'm'; it is recommended to always drop intercepts
+# when using these trend models
+mod <- mvgam(y ~ 0,
+             trend_model = PW(growth = 'logistic'),
+             family = poisson(),
+             data = mod_data)
+summary(mod)
+
+# Plot the posterior hindcast
+plot(mod, type = 'forecast')
+
+# View the changepoints with ggplot2 utilities
+library(ggplot2)
+mcmc_plot(mod, variable = 'delta_trend',
+          regex = TRUE) +
+scale_y_discrete(labels = mod$trend_model$changepoints) +
+labs(y = 'Potential changepoint',
+     x = 'Rate change')
 }
 }
 \references{
diff --git a/man/mvgam_trends.Rd b/man/mvgam_trends.Rd
index 35fc942a..85ee8e17 100644
--- a/man/mvgam_trends.Rd
+++ b/man/mvgam_trends.Rd
@@ -14,6 +14,7 @@ and the observation process is the only source of error; similarly to what is es
 \item \code{RW()}
 \item \verb{AR(p = 1, 2, or 3)}
 \item \code{VAR()}(only available in \code{Stan})
+\item \code{PW()} (piecewise linear or logistic trends; only available in \code{Stan})
 \item \code{'GP'} (Gaussian Process with squared exponential kernel;
 only available in \code{Stan})}
 
@@ -46,6 +47,8 @@ currently supported is:
 \item 'VARMA'
 \item 'VARMAcor'
 \item 'VARMA1,1cor'
+\item 'PWlinear'
+\item 'PWlogistic'
 \item 'GP'
 \item 'None'
 }
@@ -61,12 +64,17 @@ the latent process is enforced through the prior using the parameterisation give
 Heaps (2022). Stationarity is not enforced when using \code{AR1}, \code{AR2} or \code{AR3} models,
 though this can be changed by the user by specifying lower and upper bounds on autoregressive
 parameters using functionality in \link{get_mvgam_priors} and the \code{priors} argument in
-\link{mvgam}
+\link{mvgam}. Piecewise trends follow the formulation in the popular \code{prophet} package produced
+by \code{Facebook}, where users can allow for changepoints to control the potential flexibility
+of the trend. See Taylor and Letham (2018) for details
 }
 \references{
 Sarah E. Heaps (2022) Enforcing stationarity through the prior in Vector Autoregressions.
 Journal of Computational and Graphical Statistics. 32:1, 1-10.
+
+Sean J. Taylor and Benjamin Letham (2018) Forecasting at scale.
+The American Statistician 72.1, 37-45.
 }
 \seealso{
-\code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}}
+\code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}}, \code{\link{PW}}
 }
diff --git a/man/piecewise_trends.Rd b/man/piecewise_trends.Rd
index 6482644c..7789943a 100644
--- a/man/piecewise_trends.Rd
+++ b/man/piecewise_trends.Rd
@@ -26,8 +26,9 @@ Large values will allow many changepoints and a more flexible trend, while
 small values will allow few changepoints. Default is \code{0.05}.}
 
 \item{growth}{Character string specifying either 'linear' or 'logistic' growth of
-the trend. If 'logistic', a variable labelled 'cap' MUST be in \code{data} to specify the
-maximum saturation point for the trend (see examples in \code{\link{mvgam}} for details).
+the trend. If 'logistic', a variable labelled \code{cap} MUST be in \code{data} to specify the
+maximum saturation point for the trend (see details and examples in \code{\link{mvgam}} for
+more information).
 Default is 'linear'.}
 }
 \value{
@@ -40,3 +41,27 @@ in \code{mvgam}. These functions do not evaluate their arguments –
 they exist purely to help set up a model with particular piecewise
 trend models.
 }
+\details{
+\emph{Offsets and intercepts}:
+For each of these trend models, an offset parameter is included in the trend
+estimation process. This parameter will be incredibly difficult to identify
+if you also include an intercept in the observation formula. For that reason,
+it is highly recommended that you drop the intercept from the formula
+(i.e. \code{y ~ x + 0} or \code{y ~ x - 1}, where \code{x} are your optional predictor terms).
+\cr
+\cr
+\emph{Logistic growth and the cap variable}:
+When forecasting growth, there is often some maximum achievable point that
+a time series can reach. For example, total market size, total population size
+or carrying capacity in population dynamics. It can be advantageous for the forecast
+to saturate at or near this point so that predictions are more sensible.
+This function allows you to make forecasts using a logistic growth trend model,
+with a specified carrying capacity. Note that this capacity does not need to be static
+over time, it can vary with each series x timepoint combination if necessary. But you
+must supply a \code{cap} value for each observation in the data when using \code{growth = 'logistic'}.
+For observation families that use a non-identity link function, the \code{cap} value will
+be internally transformed to the link scale (i.e. your specified \code{cap} will be log
+transformed if you are using a \code{poisson()} or \code{nb()} family). It is therefore important
+that you specify the \code{cap} values on the scale of your outcome. Note also that
+no missing values are allowed in \code{cap}.
+}
diff --git a/man/update.mvgam.Rd b/man/update.mvgam.Rd
index 532918da..171a1a14 100644
--- a/man/update.mvgam.Rd
+++ b/man/update.mvgam.Rd
@@ -57,10 +57,11 @@ and the observation process is the only source of error; similarly to what is es
 \item \code{'AR2'} or \code{AR(p = 2)}
 \item \code{'AR3'} or \code{AR(p = 3)}
 \item \code{'VAR1'}  or \code{VAR()}(only available in \code{Stan})
+\item \verb{'PWlogistic}, \verb{'PWlinear} or \code{PW()} (only available in \code{Stan})
 \item \code{'GP'} (Gaussian Process with squared exponential kernel;
 only available in \code{Stan})}
 
-For all types apart from \code{'GP'}, moving average and/or correlated
+For all types apart from \code{'GP'} and \code{PW()}, moving average and/or correlated
 process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a
 multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details}
 
diff --git a/src/mvgam.dll b/src/mvgam.dll
index ee1d68fd..18c35555 100644
Binary files a/src/mvgam.dll and b/src/mvgam.dll differ
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index 8e286ed8..ca3e39fd 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ
diff --git a/tests/testthat/test-piecewise.R b/tests/testthat/test-piecewise.R
new file mode 100644
index 00000000..d061912d
--- /dev/null
+++ b/tests/testthat/test-piecewise.R
@@ -0,0 +1,67 @@
+context("piecewise")
+# Simulate data from a piecewise logistic trend
+ts <- mvgam:::piecewise_logistic(t = 1:100,
+                                 cap = 8.5,
+                                 deltas = extraDistr::rlaplace(10, 0, 0.025),
+                                 k = 0.075,
+                                 m = 0,
+                                 changepoint_ts = sample(1:100, 10))
+y <- rnorm(100, ts, 0.75)
+
+# Don't put 'cap' variable in dataframe
+df <- data.frame(y = y,
+                 time = 1:100,
+                 series = as.factor('series1'),
+                 #cap = 8.75,
+                 fake = rnorm(100))
+
+test_that("logistic should error if cap is missing", {
+  expect_error(mvgam(formula = y ~ 0,
+                     data = df,
+                     trend_model = PW(growth = 'logistic',
+                                      n_changepoints = 10),
+                     priors = prior(normal(2, 5), class = k_trend),
+                     family = gaussian(),
+                     run_model = TRUE,
+                     return_model_data = TRUE),
+               'Capacities must be supplied as a variable named "cap" for logistic growth')
+})
+
+# Now include some missing values in 'cap'
+df <- data.frame(y = y,
+                 time = 1:100,
+                 series = as.factor('series1'),
+                 cap = sample(c(8.75, NA), 100, TRUE),
+                 fake = rnorm(100))
+
+test_that("logistic should error if cap has NAs", {
+  expect_error(mvgam(formula = y ~ 0,
+                     data = df,
+                     trend_model = PW(growth = 'logistic',
+                                      n_changepoints = 10),
+                     priors = prior(normal(2, 5), class = k_trend),
+                     family = gaussian(),
+                     run_model = TRUE,
+                     return_model_data = TRUE),
+               'Missing values found for some "cap" terms')
+})
+
+# Missing values can also happen when transforming to the link scale
+y <- rpois(100, ts + 5)
+df <- data.frame(y = y,
+                 time = 1:100,
+                 series = as.factor('series1'),
+                 cap = -1,
+                 fake = rnorm(100))
+
+test_that("logistic should error if cap has NAs after link transformation", {
+  expect_error(mvgam(formula = y ~ 0,
+                     data = df,
+                     trend_model = PW(growth = 'logistic',
+                                      n_changepoints = 10),
+                     family = poisson(),
+                     run_model = TRUE,
+                     return_model_data = TRUE),
+               paste0('Missing or infinite values found for some "cap" terms\n',
+               'after transforming to the log link scale'))
+})