diff --git a/DESCRIPTION b/DESCRIPTION
index 6e74072e1..ba18f9b6f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,8 +2,8 @@ Package: brms
 Encoding: UTF-8
 Type: Package
 Title: Bayesian Regression Models using 'Stan'
-Version: 2.18.4
-Date: 2022-11-24
+Version: 2.18.5
+Date: 2022-12-13
 Authors@R:
     c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com",
              role = c("aut", "cre")),
diff --git a/NAMESPACE b/NAMESPACE
index f472ae4a6..40cea4c08 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -584,6 +584,7 @@ export(t2)
 export(theme_black)
 export(theme_default)
 export(threading)
+export(unstr)
 export(update_adterms)
 export(validate_newdata)
 export(validate_prior)
diff --git a/NEWS.md b/NEWS.md
index 81a705ac2..711273220 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,8 @@
 
 ### New Features
 
+* Model unstructured autocorrelation matrices via the `unstr` term
+thanks to the help of Sebastian Weber. (#1435)
 * Improve user control over model recompilation via argument `recompile`
 in post-processing methods that require a compiled Stan model.
 * Extend control over the `point_estimate` feature in `prepare_predictions`
diff --git a/R/brmsfit-helpers.R b/R/brmsfit-helpers.R
index ef07abc88..ff88f0344 100644
--- a/R/brmsfit-helpers.R
+++ b/R/brmsfit-helpers.R
@@ -214,8 +214,9 @@ get_cor_matrix <- function(cor, size = NULL, ndraws = NULL) {
 # compute covariance matrices of autocor structures
 # @param prep a brmsprep object
 # @param obs observations for which to compute the covariance matrix
+# @param Jtime vector indicating to which time points obs belong
 # @param latent compute covariance matrix for latent residuals?
-get_cov_matrix_ac <- function(prep, obs = NULL, latent = FALSE) {
+get_cov_matrix_ac <- function(prep, obs = NULL, Jtime = NULL, latent = FALSE) {
   if (is.null(obs)) {
     obs <- seq_len(prep$nobs)
   }
@@ -238,6 +239,9 @@ get_cov_matrix_ac <- function(prep, obs = NULL, latent = FALSE) {
   } else if (has_ac_class(acef, "cosy")) {
     cosy <- as.numeric(prep$ac$cosy)
     cor <- get_cor_matrix_cosy(cosy, nobs)
+  } else if (has_ac_class(acef, "unstr")) {
+    cortime <- prep$ac$cortime
+    cor <- get_cor_matrix_unstr(cortime, Jtime)
   } else if (has_ac_class(acef, "fcor")) {
     cor <- get_cor_matrix_fcor(prep$ac$Mfcor, ndraws)
   } else {
@@ -347,6 +351,17 @@ get_cor_matrix_cosy <- function(cosy, nobs) {
   out
 }
 
+# compute unstructured time correlation matrices
+# @param cortime time correlation draws
+# @param Jtime indictor of rows/cols to consider in cortime
+# @return a numeric 'ndraws' x 'nobs' x 'nobs' array
+#   where nobs = length(Jtime[Jtime > 0])
+get_cor_matrix_unstr <- function(cortime, Jtime) {
+  stopifnot(length(Jtime) > 0L)
+  Jtime <- Jtime[Jtime > 0]
+  get_cor_matrix(cortime)[, Jtime, Jtime, drop = FALSE]
+}
+
 # prepare a fixed correlation matrix
 # @param Mfcor correlation matrix to be prepared
 # @param ndraws number of posterior draws
diff --git a/R/brmsterms.R b/R/brmsterms.R
index f7fefbbe9..c86d8ef4b 100644
--- a/R/brmsterms.R
+++ b/R/brmsterms.R
@@ -903,7 +903,7 @@ get_effect.btl <- function(x, target = "fe", ...) {
 
 #' @export
 get_effect.btnl <- function(x, target = "fe", ...) {
-  NULL
+  x[[target]]
 }
 
 all_terms <- function(x) {
@@ -924,7 +924,7 @@ regex_sp <- function(type = "all") {
   funs <- c(
     sm = "(s|(t2)|(te)|(ti))",
     gp = "gp", cs = "cse?", mmc = "mmc",
-    ac = "((arma)|(ar)|(ma)|(cosy)|(sar)|(car)|(fcor))"
+    ac = "((arma)|(ar)|(ma)|(cosy)|(unstr)|(sar)|(car)|(fcor))"
   )
   funs[all_sp_types()] <- all_sp_types()
   if ("sp" %in% type) {
diff --git a/R/data-predictor.R b/R/data-predictor.R
index 37c0e12df..d5c1efd6c 100644
--- a/R/data-predictor.R
+++ b/R/data-predictor.R
@@ -631,14 +631,13 @@ data_gp <- function(bterms, data, internal = FALSE, basis = NULL, ...) {
 }
 
 # data for autocorrelation variables
-# @param locations optional original locations for CAR models
 data_ac <- function(bterms, data, data2, basis = NULL, ...) {
   out <- list()
   N <- nrow(data)
   acef <- tidy_acef(bterms)
   if (has_ac_subset(bterms, dim = "time")) {
-    gr <- subset2(acef, dim = "time")$gr
-    if (gr != "NA") {
+    gr <- get_ac_vars(acef, "gr", dim = "time")
+    if (isTRUE(nzchar(gr))) {
       tgroup <- as.numeric(factor(data[[gr]]))
     } else {
       tgroup <- rep(1, N)
@@ -662,12 +661,37 @@ data_ac <- function(bterms, data, data2, basis = NULL, ...) {
   }
   if (use_ac_cov_time(acef)) {
     # data for the 'covariance' versions of time-series structures
+    # TODO: change begin[i]:end[i] notation to slice[i]:(slice[i+1] - 1)
+    #   see comment on PR #1435
     out$N_tg <- length(unique(tgroup))
     out$begin_tg <- as.array(ulapply(unique(tgroup), match, tgroup))
     out$nobs_tg <- as.array(with(out,
       c(if (N_tg > 1L) begin_tg[2:N_tg], N + 1) - begin_tg
     ))
     out$end_tg <- with(out, begin_tg + nobs_tg - 1)
+    if (has_ac_class(acef, "unstr")) {
+      time <- get_ac_vars(bterms, "time", dim = "time")
+      time_data <- get(time, data)
+      new_times <- extract_levels(time_data)
+      if (length(basis)) {
+        times <- basis$times
+        # unstr estimates correlations only for given time points
+        invalid_times <- setdiff(new_times, times)
+        if (length(invalid_times)) {
+          stop2("Cannot handle new time points in UNSTR models.")
+        }
+      } else {
+        times <- new_times
+      }
+      out$n_unique_t <- length(times)
+      out$n_unique_cortime <- out$n_unique_t * (out$n_unique_t - 1) / 2
+      Jtime <- match(time_data, times)
+      out$Jtime_tg <- matrix(0L, out$N_tg, max(out$nobs_tg))
+      for (i in seq_len(out$N_tg)) {
+        out$Jtime_tg[i, seq_len(out$nobs_tg[i])] <-
+          Jtime[out$begin_tg[i]:out$end_tg[i]]
+      }
+    }
   }
   if (has_ac_class(acef, "sar")) {
     acef_sar <- subset2(acef, class = "sar")
diff --git a/R/exclude_pars.R b/R/exclude_pars.R
index 5899b320f..c47bc299b 100644
--- a/R/exclude_pars.R
+++ b/R/exclude_pars.R
@@ -74,9 +74,13 @@ exclude_pars.mvbrmsterms <- function(x, save_pars, ...) {
 exclude_pars.brmsterms <- function(x, data, save_pars, ...) {
   resp <- usc(combine_prefix(x))
   data <- subset_data(data, x)
-  out <- "Lncor"
+  par_classes <- c("Lncor", "Cortime")
+  out <- paste0(par_classes, resp)
   if (!save_pars$all) {
-    par_classes <- c("ordered_Intercept", "fixed_Intercept", "theta", "Llncor")
+    par_classes <- c(
+      "ordered_Intercept", "fixed_Intercept",
+      "theta", "Llncor", "Lcortime"
+    )
     c(out) <- paste0(par_classes, resp)
   }
   for (dp in names(x$dpars)) {
diff --git a/R/formula-ac.R b/R/formula-ac.R
index c4a6d5db6..e9289259b 100644
--- a/R/formula-ac.R
+++ b/R/formula-ac.R
@@ -2,11 +2,11 @@
 #'
 #' Specify autocorrelation terms in \pkg{brms} models. Currently supported terms
 #' are \code{\link{arma}}, \code{\link{ar}}, \code{\link{ma}},
-#' \code{\link{cosy}}, \code{\link{sar}}, \code{\link{car}}, and
-#' \code{\link{fcor}}. Terms can be directly specified within the formula, or
-#' passed to the \code{autocor} argument of \code{\link{brmsformula}} in the
-#' form of a one-sided formula. For deprecated ways of specifying
-#' autocorrelation terms, see \code{\link{cor_brms}}.
+#' \code{\link{cosy}}, \code{\link{unstr}}, \code{\link{sar}},
+#' \code{\link{car}}, and \code{\link{fcor}}. Terms can be directly specified
+#' within the formula, or passed to the \code{autocor} argument of
+#' \code{\link{brmsformula}} in the form of a one-sided formula. For deprecated
+#' ways of specifying autocorrelation terms, see \code{\link{cor_brms}}.
 #'
 #' @name autocor-terms
 #'
@@ -17,8 +17,8 @@
 #'
 #' @seealso \code{\link{brmsformula}}, \code{\link{acformula}},
 #'   \code{\link{arma}}, \code{\link{ar}}, \code{\link{ma}},
-#'   \code{\link{cosy}}, \code{\link{sar}}, \code{\link{car}},
-#'   \code{\link{fcor}}
+#'   \code{\link{cosy}}, \code{\link{unstr}}, \code{\link{sar}},
+#'   \code{\link{car}}, \code{\link{fcor}}
 #'
 #' @examples
 #' # specify autocor terms within the formula
@@ -188,7 +188,6 @@ ma <- function(time = NA, gr = NA, q = 1, cov = FALSE) {
 #' }
 #'
 #' @export
-#' @export
 cosy <- function(time = NA, gr = NA) {
   label <- deparse(match.call())
   time <- deparse(substitute(time))
@@ -200,6 +199,39 @@ cosy <- function(time = NA, gr = NA) {
   out
 }
 
+#' Set up UNSTR correlation structures
+#'
+#' Set up an unstructured (UNSTR) correlation term in \pkg{brms}. The function does
+#' not evaluate its arguments -- it exists purely to help set up a model with
+#' UNSTR terms.
+#'
+#' @inheritParams arma
+#'
+#' @return An object of class \code{'unstr_term'}, which is a list
+#'   of arguments to be interpreted by the formula
+#'   parsing functions of \pkg{brms}.
+#'
+#' @seealso \code{\link{autocor-terms}}
+#'
+#' @examples
+#' \dontrun{
+#' # add an unstructured correlation matrix for visits within the same patient
+#' fit <- brm(count ~ Trt + unstr(visit, patient), data = epilepsy)
+#' summary(fit)
+#' }
+#'
+#' @export
+unstr <- function(time, gr) {
+  label <- deparse(match.call())
+  time <- deparse(substitute(time))
+  time <- as_one_variable(time)
+  gr <- deparse(substitute(gr))
+  stopif_illegal_group(gr)
+  out <- nlist(time, gr, label)
+  class(out) <- c("unstr_term", "ac_term")
+  out
+}
+
 #' Spatial simultaneous autoregressive (SAR) structures
 #'
 #' Set up an spatial simultaneous autoregressive (SAR) term in \pkg{brms}. The
@@ -443,7 +475,7 @@ tidy_acef.brmsterms <- function(x, ...) {
 }
 
 #' @export
-tidy_acef.btl <- function(x, data = NULL, ...) {
+tidy_acef.btl <- function(x, ...) {
   form <- x[["ac"]]
   if (!is.formula(form)) {
     return(empty_acef())
@@ -476,6 +508,13 @@ tidy_acef.btl <- function(x, data = NULL, ...) {
       out$gr[i] <- ac$gr
       out$cov[i] <- TRUE
     }
+    if (is.unstr_term(ac)) {
+      out$class[i] <- "unstr"
+      out$dim[i] <- "time"
+      out$time[i] <- ac$time
+      out$gr[i] <- ac$gr
+      out$cov[i] <- TRUE
+    }
     if (is.sar_term(ac)) {
       out$class[i] <- "sar"
       out$dim[i] <- "space"
@@ -638,6 +677,7 @@ validate_fcor_matrix <- function(M) {
 
 # regex to extract all parameter names of autocorrelation structures
 regex_autocor_pars <- function() {
+  # cortime is ignored here to allow custom renaming in summary.brmsfit
   p <- c("ar", "ma", "sderr", "cosy", "lagsar", "errorsar",
          "car", "sdcar", "rhocar")
   p <- paste0("(", p, ")", collapse = "|")
@@ -656,6 +696,10 @@ is.cosy_term <- function(x) {
   inherits(x, "cosy_term")
 }
 
+is.unstr_term <- function(x) {
+  inherits(x, "unstr_term")
+}
+
 is.sar_term <- function(x) {
   inherits(x, "sar_term")
 }
diff --git a/R/formula-re.R b/R/formula-re.R
index df841c543..b563732f0 100644
--- a/R/formula-re.R
+++ b/R/formula-re.R
@@ -796,7 +796,10 @@ get_levels <- function(...) {
 
 extract_levels <- function(x) {
   # do not check for NAs according to #1355
-  levels(factor(x))
+  if (!is.factor(x)) {
+    x <- factor(x)
+  }
+  levels(x)
 }
 
 # extract names of group-level effects
diff --git a/R/log_lik.R b/R/log_lik.R
index 35a1def1c..32fe27329 100644
--- a/R/log_lik.R
+++ b/R/log_lik.R
@@ -245,9 +245,10 @@ log_lik_student_mv <- function(i, prep) {
 
 log_lik_gaussian_time <- function(i, prep) {
   obs <- with(prep$ac, begin_tg[i]:end_tg[i])
+  Jtime <- prep$ac$Jtime_tg[i, ]
   Y <- as.numeric(prep$data$Y[obs])
   mu <- as.matrix(get_dpar(prep, "mu", i = obs))
-  Sigma <- get_cov_matrix_ac(prep, obs)
+  Sigma <- get_cov_matrix_ac(prep, obs, Jtime = Jtime)
   .log_lik <- function(s) {
     C <- as.matrix(Sigma[s, , ])
     Cinv <- solve(C)
@@ -264,10 +265,11 @@ log_lik_gaussian_time <- function(i, prep) {
 
 log_lik_student_time <- function(i, prep) {
   obs <- with(prep$ac, begin_tg[i]:end_tg[i])
+  Jtime <- prep$ac$Jtime_tg[i, ]
   Y <- as.numeric(prep$data$Y[obs])
   nu <- as.matrix(get_dpar(prep, "nu", i = obs))
   mu <- as.matrix(get_dpar(prep, "mu", i = obs))
-  Sigma <- get_cov_matrix_ac(prep, obs)
+  Sigma <- get_cov_matrix_ac(prep, obs, Jtime = Jtime)
   .log_lik <- function(s) {
     df <- nu[s, ]
     C <- as.matrix(Sigma[s, , ])
diff --git a/R/make_standata.R b/R/make_standata.R
index 18b59becc..452559077 100644
--- a/R/make_standata.R
+++ b/R/make_standata.R
@@ -283,13 +283,16 @@ standata_basis_ac <- function(x, data, ...) {
   out <- list()
   if (has_ac_class(x, "car")) {
     gr <- get_ac_vars(x, "gr", class = "car")
-    stopifnot(length(gr) <= 1L)
     if (isTRUE(nzchar(gr))) {
       out$locations <- extract_levels(get(gr, data))
     } else {
       out$locations <- NA
     }
   }
+  if (has_ac_class(x, "unstr")) {
+    time <- get_ac_vars(x, "time", dim = "time")
+    out$times <- extract_levels(get(time, data))
+  }
   out
 }
 
diff --git a/R/posterior_predict.R b/R/posterior_predict.R
index 417f65447..e19be009c 100644
--- a/R/posterior_predict.R
+++ b/R/posterior_predict.R
@@ -389,8 +389,9 @@ posterior_predict_student_mv <- function(i, prep, ...) {
 
 posterior_predict_gaussian_time <- function(i, prep, ...) {
   obs <- with(prep$ac, begin_tg[i]:end_tg[i])
+  Jtime <- prep$ac$Jtime_tg[i, ]
   mu <- as.matrix(get_dpar(prep, "mu", i = obs))
-  Sigma <- get_cov_matrix_ac(prep, obs)
+  Sigma <- get_cov_matrix_ac(prep, obs, Jtime = Jtime)
   .predict <- function(s) {
     rmulti_normal(1, mu = mu[s, ], Sigma = Sigma[s, , ])
   }
@@ -399,9 +400,10 @@ posterior_predict_gaussian_time <- function(i, prep, ...) {
 
 posterior_predict_student_time <- function(i, prep, ...) {
   obs <- with(prep$ac, begin_tg[i]:end_tg[i])
+  Jtime <- prep$ac$Jtime_tg[i, ]
   nu <- as.matrix(get_dpar(prep, "nu", i = obs))
   mu <- as.matrix(get_dpar(prep, "mu", i = obs))
-  Sigma <- get_cov_matrix_ac(prep, obs)
+  Sigma <- get_cov_matrix_ac(prep, obs, Jtime = Jtime)
   .predict <- function(s) {
     rmulti_student_t(1, df = nu[s, ], mu = mu[s, ], Sigma = Sigma[s, , ])
   }
diff --git a/R/predictor.R b/R/predictor.R
index 9ae7e8d65..276abee4b 100644
--- a/R/predictor.R
+++ b/R/predictor.R
@@ -440,21 +440,19 @@ predictor_offset <- function(prep, i, nobs) {
 # @note eta has to be passed to this function in
 #   order for ARMA structures to work correctly
 predictor_ac <- function(eta, prep, i, fprep = NULL) {
-  if (has_ac_class(prep$ac$acef, "arma")) {
-    if (!is.null(prep$ac$err)) {
-      # ARMA correlations via latent residuals
-      eta <- eta + p(prep$ac$err, i, row = FALSE)
-    } else {
-      # ARMA correlations via explicit natural residuals
-      if (!is.null(i)) {
-        stop2("Pointwise evaluation is not possible for ARMA models.")
-      }
-      eta <- .predictor_arma(
-        eta, ar = prep$ac$ar, ma = prep$ac$ma,
-        Y = prep$ac$Y, J_lag = prep$ac$J_lag,
-        fprep = fprep
-      )
+  if (!is.null(prep$ac[["err"]])) {
+    # auto-correlations via latent residuals
+    eta <- eta + p(prep$ac$err, i, row = FALSE)
+  } else if (has_ac_class(prep$ac$acef, "arma")) {
+    # ARMA correlations via explicit natural residuals
+    if (!is.null(i)) {
+      stop2("Pointwise evaluation is not possible for ARMA models.")
     }
+    eta <- .predictor_arma(
+      eta, ar = prep$ac$ar, ma = prep$ac$ma,
+      Y = prep$ac$Y, J_lag = prep$ac$J_lag,
+      fprep = fprep
+    )
   }
   if (has_ac_class(prep$ac$acef, "car")) {
     eta <- eta + .predictor_re(Z = p(prep$ac$Zcar, i), r = prep$ac$rcar)
diff --git a/R/prepare_predictions.R b/R/prepare_predictions.R
index c71df5eda..dc1b9499b 100644
--- a/R/prepare_predictions.R
+++ b/R/prepare_predictions.R
@@ -736,6 +736,11 @@ prepare_predictions_ac <- function(bterms, draws, sdata, oos = NULL,
     cosy_regex <- paste0("^cosy", p, "$")
     out$cosy <- prepare_draws(draws, cosy_regex, regex = TRUE)
   }
+  if (has_ac_class(acef, "unstr")) {
+    cortime_regex <- paste0("^cortime", p, "__")
+    out$cortime <- prepare_draws(draws, cortime_regex, regex = TRUE)
+    out$Jtime_tg <- sdata[[paste0("Jtime_tg", p)]]
+  }
   if (use_ac_cov_time(acef)) {
     # prepare predictions for the covariance structures of time-series models
     out$begin_tg <- sdata[[paste0("begin_tg", p)]]
@@ -757,8 +762,9 @@ prepare_predictions_ac <- function(bterms, draws, sdata, oos = NULL,
       out$sderr <- prepare_draws(draws, sderr_regex, regex = TRUE)
       for (i in seq_len(out$N_tg)) {
         obs <- with(out, begin_tg[i]:end_tg[i])
+        Jtime <- out$Jtime_tg[i, ]
+        cov <- get_cov_matrix_ac(list(ac = out), obs, Jtime = Jtime, latent = TRUE)
         zeros <- rep(0, length(obs))
-        cov <- get_cov_matrix_ac(list(ac = out), obs, latent = TRUE)
         .err <- function(s) rmulti_normal(1, zeros, Sigma = cov[s, , ])
         out$err[, obs] <- rblapply(seq_rows(draws), .err)
       }
diff --git a/R/priors.R b/R/priors.R
index 12edfbb55..07e47b6d6 100644
--- a/R/priors.R
+++ b/R/priors.R
@@ -588,7 +588,8 @@ prior_predictor.brmsterms <- function(x, data, internal = FALSE, ...) {
       dp_prior <- prior_predictor(
         x$dpars[[dp]], data = data,
         def_scale_prior = def_scale_prior,
-        def_dprior = def_dprior
+        def_dprior = def_dprior,
+        internal = internal
       )
     } else if (!is.null(x$fdpars[[dp]])) {
       # parameter is fixed
@@ -608,7 +609,8 @@ prior_predictor.brmsterms <- function(x, data, internal = FALSE, ...) {
     nlp_prior <- prior_predictor(
       x$nlpars[[nlp]], data = data,
       def_scale_prior = def_scale_prior,
-      def_dprior = def_dprior
+      def_dprior = def_dprior,
+      internal = internal
     )
     prior <- prior + nlp_prior
   }
@@ -964,7 +966,7 @@ prior_sm <- function(bterms, data, def_scale_prior, ...) {
 }
 
 # priors for autocor parameters
-prior_ac <- function(bterms, def_scale_prior, ...) {
+prior_ac <- function(bterms, def_scale_prior, internal = FALSE, ...) {
   prior <- empty_prior()
   acef <- tidy_acef(bterms)
   if (!NROW(acef)) {
@@ -995,6 +997,15 @@ prior_ac <- function(bterms, def_scale_prior, ...) {
     prior <- prior +
       brmsprior(class = "cosy", ls = px, lb = "0", ub = "1")
   }
+  if (has_ac_class(acef, "unstr")) {
+    if (internal) {
+      prior <- prior +
+        brmsprior("lkj_corr_cholesky(1)", class = "Lcortime", ls = px)
+    } else {
+      prior <- prior +
+        brmsprior("lkj(1)", class = "cortime", ls = px)
+    }
+  }
   if (has_ac_latent_residuals(bterms)) {
     prior <- prior +
       brmsprior(def_scale_prior, class = "sderr", ls = px, lb = "0")
@@ -1195,8 +1206,8 @@ validate_prior <- function(prior, formula, data, family = gaussian(),
   prior <- prior[!no_checks, ]
   # check for duplicated priors
   prior$class <- rename(
-    prior$class, c("^cor$", "^rescor$", "^corme$", "^lncor"),
-    c("L", "Lrescor", "Lme", "Llncor"), fixed = FALSE
+    prior$class, c("^cor$", "^rescor$", "^corme$", "^lncor", "^cortime"),
+    c("L", "Lrescor", "Lme", "Llncor", "Lcortime"), fixed = FALSE
   )
   if (any(duplicated(prior))) {
     stop2("Duplicated prior specifications are not allowed.")
@@ -1307,7 +1318,10 @@ check_prior_content <- function(prior) {
   lb_priors_regex <- paste0("^(", paste0(lb_priors, collapse = "|"), ")")
   ulb_priors <- c("beta", "uniform", "von_mises", "beta_proportion")
   ulb_priors_regex <- paste0("^(", paste0(ulb_priors, collapse = "|"), ")")
-  cormat_pars <- c("cor", "L", "rescor", "Lrescor", "corme", "Lme", "lncor", "Llncor")
+  cormat_pars <- c(
+    "cor", "L", "rescor", "Lrescor", "corme", "Lme",
+    "lncor", "Llncor", "cortime", "Lcortime"
+  )
   cormat_regex <- "^((lkj)|(constant))"
   simplex_pars <- c("simo", "theta", "sbhaz")
   simplex_regex <- "^((dirichlet)|(constant))\\("
diff --git a/R/rename_pars.R b/R/rename_pars.R
index 6e2d3d2c3..3f7ddba7f 100644
--- a/R/rename_pars.R
+++ b/R/rename_pars.R
@@ -98,7 +98,8 @@ change_effects.btl <- function(x, ...) {
     change_sm(x, ...),
     change_cs(x, ...),
     change_sp(x, ...),
-    change_gp(x, ...))
+    change_gp(x, ...),
+    change_ac(x, ...))
 }
 
 # helps in renaming fixed effects parameters
@@ -335,6 +336,19 @@ change_sm <- function(bterms, data, pars, ...) {
   out
 }
 
+# helps in renaming autocorrelation parameters
+change_ac <- function(bterms, data, pars, ...) {
+  out <- list()
+  acef <- tidy_acef(bterms)
+  if (has_ac_class(acef, "unstr")) {
+    time <- get_ac_vars(acef, "time", dim = "time")
+    times <- extract_levels(get(time, data))
+    cortime_names <- get_cornames(times, type = "cortime", brackets = FALSE)
+    lc(out) <- clist(grepl("^cortime\\[", pars), cortime_names)
+  }
+  out
+}
+
 # helps in renaming group-level parameters
 # @param ranef: data.frame returned by 'tidy_ranef'
 change_re <- function(ranef, pars, ...) {
@@ -545,8 +559,8 @@ do_renaming <- function(x, change) {
 # @param x brmsfit object
 reorder_pars <- function(x) {
   all_classes <- unique(c(
-    "b", "bs", "bsp", "bcs", "ar", "ma", "sderr", "lagsar", "errorsar",
-    "car", "rhocar", "sdcar", "cosy", "sd", "cor", "df", "sds", "sdgp",
+    "b", "bs", "bsp", "bcs", "ar", "ma", "sderr", "lagsar", "errorsar", "car",
+    "rhocar", "sdcar", "cosy", "cortime", "sd", "cor", "df", "sds", "sdgp",
     "lscale", valid_dpars(x), "lncor", "Intercept", "tmp", "rescor",
     "delta", "lasso", "simo", "r", "s", "zgp", "rcar", "sbhaz",
     "R2D2", "Ymi", "Yl", "meanme", "sdme", "corme", "Xme", "prior",
diff --git a/R/stan-helpers.R b/R/stan-helpers.R
index f76197e01..7b7e9d705 100644
--- a/R/stan-helpers.R
+++ b/R/stan-helpers.R
@@ -3,6 +3,8 @@
 
 # define Stan functions or globally used transformed data
 # TODO: refactor to not require extraction of information from all model parts
+#   'expand_include_statements' removes duplicates which opens the door
+#   for adding Stan functions at better places rather than globally here
 stan_global_defs <- function(bterms, prior, ranef, threads) {
   families <- family_names(bterms)
   links <- family_info(bterms, "link")
@@ -77,16 +79,39 @@ stan_global_defs <- function(bterms, prior, ranef, threads) {
   acterms <- get_effect(bterms, "ac")
   acefs <- lapply(acterms, tidy_acef)
   if (any(ulapply(acefs, has_ac_subset, dim = "time", cov = TRUE))) {
-    # TODO: include functions selectively
     str_add(out$fun) <- glue(
-      "  #include 'fun_normal_time.stan'\n",
-      "  #include 'fun_student_t_time.stan'\n",
-      "  #include 'fun_scale_time_err.stan'\n",
-      "  #include 'fun_cholesky_cor_ar1.stan'\n",
-      "  #include 'fun_cholesky_cor_ma1.stan'\n",
-      "  #include 'fun_cholesky_cor_arma1.stan'\n",
-      "  #include 'fun_cholesky_cor_cosy.stan'\n"
+      "  #include 'fun_sequence.stan'\n",
+      "  #include 'fun_is_equal.stan'\n",
+      "  #include 'fun_stack_vectors.stan'\n"
     )
+    if ("gaussian" %in% families) {
+      str_add(out$fun) <- glue(
+        "  #include 'fun_normal_time.stan'\n",
+        "  #include 'fun_normal_time_se.stan'\n"
+      )
+    }
+    if ("student" %in% families) {
+      str_add(out$fun) <- glue(
+        "  #include 'fun_student_t_time.stan'\n",
+        "  #include 'fun_student_t_time_se.stan'\n"
+      )
+    }
+    # TODO: include selectively once we have the 'latent' indicator
+    str_add(out$fun) <- glue(
+      "  #include 'fun_scale_time_err.stan'\n"
+    )
+    if (any(ulapply(acefs, has_ac_class, "arma"))) {
+      str_add(out$fun) <- glue(
+        "  #include 'fun_cholesky_cor_ar1.stan'\n",
+        "  #include 'fun_cholesky_cor_ma1.stan'\n",
+        "  #include 'fun_cholesky_cor_arma1.stan'\n"
+      )
+    }
+    if (any(ulapply(acefs, has_ac_class, "cosy"))) {
+      str_add(out$fun) <- glue(
+        "  #include 'fun_cholesky_cor_cosy.stan'\n"
+      )
+    }
   }
   if (any(ulapply(acefs, has_ac_class, "sar"))) {
     if ("gaussian" %in% families) {
diff --git a/R/stan-likelihood.R b/R/stan-likelihood.R
index f28402bea..d2c3daea3 100644
--- a/R/stan-likelihood.R
+++ b/R/stan-likelihood.R
@@ -368,13 +368,23 @@ stan_log_lik_gaussian_time <- function(bterms, resp = "", mix = "", ...) {
   if (stan_log_lik_adj(bterms)) {
     stop2("Invalid addition arguments for this model.")
   }
+  has_se <- is.formula(bterms$adforms$se)
+  flex <- has_ac_class(tidy_acef(bterms), "unstr")
   p <- stan_log_lik_dpars(bterms, FALSE, resp, mix)
-  v <- c("chol_cor", "se2", "nobs_tg", "begin_tg", "end_tg")
+  v <- c("Lcortime", "nobs_tg", "begin_tg", "end_tg")
+  if (has_se) {
+    c(v) <- "se2"
+  }
+  if (flex) {
+    c(v) <- "Jtime_tg"
+  }
   p[v] <- as.list(paste0(v, resp))
   sfx <- str_if("sigma" %in% names(bterms$dpars), "het", "hom")
+  sfx <- str_if(has_se, paste0(sfx, "_se"), sfx)
+  sfx <- str_if(flex, paste0(sfx, "_flex"), sfx)
   sdist(glue("normal_time_{sfx}"),
-    p$mu, p$sigma, p$chol_cor, p$se2,
-    p$nobs_tg, p$begin_tg, p$end_tg
+    p$mu, p$sigma, p$se2, p$Lcortime,
+    p$nobs_tg, p$begin_tg, p$end_tg, p$Jtime_tg
   )
 }
 
@@ -427,13 +437,23 @@ stan_log_lik_student_time <- function(bterms, resp = "", mix = "", ...) {
   if (stan_log_lik_adj(bterms)) {
     stop2("Invalid addition arguments for this model.")
   }
+  has_se <- is.formula(bterms$adforms$se)
+  flex <- has_ac_class(tidy_acef(bterms), "unstr")
   p <- stan_log_lik_dpars(bterms, FALSE, resp, mix)
-  v <- c("chol_cor", "se2", "nobs_tg", "begin_tg", "end_tg")
+  v <- c("Lcortime", "nobs_tg", "begin_tg", "end_tg")
+  if (has_se) {
+    c(v) <- "se2"
+  }
+  if (flex) {
+    c(v) <- "Jtime_tg"
+  }
   p[v] <- as.list(paste0(v, resp))
   sfx <- str_if("sigma" %in% names(bterms$dpars), "het", "hom")
+  sfx <- str_if(has_se, paste0(sfx, "_se"), sfx)
+  sfx <- str_if(flex, paste0(sfx, "_flex"), sfx)
   sdist(glue("student_t_time_{sfx}"),
-    p$nu, p$mu, p$sigma, p$chol_cor, p$se2,
-    p$nobs_tg, p$begin_tg, p$end_tg
+    p$nu, p$mu, p$sigma, p$se2, p$Lcortime,
+    p$nobs_tg, p$begin_tg, p$end_tg, p$Jtime_tg
   )
 }
 
diff --git a/R/stan-predictor.R b/R/stan-predictor.R
index 7893d484d..99fca9dc4 100644
--- a/R/stan-predictor.R
+++ b/R/stan-predictor.R
@@ -1250,7 +1250,7 @@ stan_ac <- function(bterms, data, prior, threads, normalize, ...) {
   slice <- stan_slice(threads)
   has_natural_residuals <- has_natural_residuals(bterms)
   has_ac_latent_residuals <- has_ac_latent_residuals(bterms)
-  acef <- tidy_acef(bterms, data)
+  acef <- tidy_acef(bterms)
 
   if (has_ac_latent_residuals) {
     # families that do not have natural residuals require latent
@@ -1374,12 +1374,30 @@ stan_ac <- function(bterms, data, prior, threads, normalize, ...) {
     )
   }
 
+  acef_unstr <- subset2(acef, class = "unstr")
+  if (NROW(acef_unstr)) {
+    # unstructured correlation matrix
+    # most code is shared with ARMA covariance models
+    # define prior on the Cholesky scale to consistency across
+    # autocorrelation structures
+    str_add_list(out) <- stan_prior(
+      prior, class = "Lcortime", px = px, suffix = p,
+      type = glue("cholesky_factor_corr[n_unique_t{p}]"),
+      header_type = "matrix",
+      comment = "cholesky factor of unstructured autocorrelation matrix",
+      normalize = normalize
+    )
+  }
+
   acef_time_cov <- subset2(acef, dim = "time", cov = TRUE)
   if (NROW(acef_time_cov)) {
     # use correlation structures in covariance matrix parameterization
-    # optional for ARMA models and obligatory for COSY models
+    # optional for ARMA models and obligatory for COSY and UNSTR models
     # can only model one covariance structure at a time
     stopifnot(NROW(acef_time_cov) == 1)
+    if (use_threading(threads)) {
+      stop2("Threading is not supported for covariance-based autocorrelation models.")
+    }
     str_add(out$data) <- glue(
       "  // see the functions block for details\n",
       "  int<lower=1> N_tg{p};\n",
@@ -1394,44 +1412,68 @@ stan_ac <- function(bterms, data, prior, threads, normalize, ...) {
       "  int max_nobs_tg{p} = max(nobs_tg{p});",
       "  // maximum dimension of the autocorrelation matrix\n"
     )
-    if (!is.formula(bterms$adforms$se)) {
-      str_add(out$tdata_def) <- glue(
-        "  // no known standard errors specified by the user\n",
-        "  vector[N{resp}] se2{p} = rep_vector(0.0, N{resp});\n"
+    if (acef_time_cov$class == "unstr") {
+      # unstructured time-covariances require additional data and cannot
+      # be represented directly via Cholesky factors due to potentially
+      # different time subsets
+      str_add(out$data) <- glue(
+        "  int<lower=0> Jtime_tg{p}[N_tg{p}, max(nobs_tg{p})];\n",
+        "  int n_unique_t{p};  // total number of unique time points\n",
+        "  int n_unique_cortime{p};  // number of unique correlations\n"
       )
-      str_add(out$pll_args) <- glue(", data vector se2{p}")
-    }
-    str_add(out$tpar_def) <- glue(
-      "  // cholesky factor of the autocorrelation matrix\n",
-      "  matrix[max_nobs_tg{p}, max_nobs_tg{p}] chol_cor{p};\n"
-    )
-    str_add(out$pll_args) <- glue(", matrix chol_cor{p}")
-    if (acef_time_cov$class == "arma") {
-      if (acef_time_cov$p > 0 && acef_time_cov$q == 0) {
-        cor_fun <- "ar1"
-        cor_args <- glue("ar{p}[1]")
-      } else if (acef_time_cov$p == 0 && acef_time_cov$q > 0) {
-        cor_fun <- "ma1"
-        cor_args <- glue("ma{p}[1]")
-      } else {
-        cor_fun <- "arma1"
-        cor_args <- glue("ar{p}[1], ma{p}[1]")
+      str_add(out$pll_args) <- glue(", int[,] Jtime_tg{p}")
+      if (has_ac_latent_residuals) {
+        str_add(out$tpar_comp) <- glue(
+          "  // compute correlated time-series residuals\n",
+          "  err{p} = scale_time_err_flex(",
+          "zerr{p}, sderr{p}, Lcortime{p}, nobs_tg{p}, begin_tg{p}, end_tg{p}, Jtime_tg{p});\n"
+        )
       }
-    } else if (acef_time_cov$class == "cosy") {
-      cor_fun <- "cosy"
-      cor_args <- glue("cosy{p}")
-    }
-    str_add(out$tpar_comp) <- glue(
-      "  // compute residual covariance matrix\n",
-      "  chol_cor{p} = cholesky_cor_{cor_fun}({cor_args}, max_nobs_tg{p});\n"
-    )
-    if (has_ac_latent_residuals) {
+      str_add(out$gen_def) <- glue(
+        "  // compute group-level correlations\n",
+        "  corr_matrix[n_unique_t{p}] Cortime{p}",
+        " = multiply_lower_tri_self_transpose(Lcortime{p});\n",
+        "  vector<lower=-1,upper=1>[n_unique_cortime{p}] cortime{p};\n"
+      )
+      str_add(out$gen_comp) <- stan_cor_gen_comp(
+        glue("cortime{p}"), glue("n_unique_t{p}")
+      )
+    } else {
+      # all other time-covariance structures can be represented directly
+      # through Cholesky factors of the correlation matrix
+      if (acef_time_cov$class == "arma") {
+        if (acef_time_cov$p > 0 && acef_time_cov$q == 0) {
+          cor_fun <- "ar1"
+          cor_args <- glue("ar{p}[1]")
+        } else if (acef_time_cov$p == 0 && acef_time_cov$q > 0) {
+          cor_fun <- "ma1"
+          cor_args <- glue("ma{p}[1]")
+        } else {
+          cor_fun <- "arma1"
+          cor_args <- glue("ar{p}[1], ma{p}[1]")
+        }
+      } else if (acef_time_cov$class == "cosy") {
+        cor_fun <- "cosy"
+        cor_args <- glue("cosy{p}")
+      }
+      str_add(out$tpar_def) <- glue(
+        "  // cholesky factor of the autocorrelation matrix\n",
+        "  matrix[max_nobs_tg{p}, max_nobs_tg{p}] Lcortime{p};\n"
+      )
+      str_add(out$pll_args) <- glue(", matrix Lcortime{p}")
       str_add(out$tpar_comp) <- glue(
-        "  // compute correlated time-series residuals\n",
-        "  err{p} = scale_time_err(",
-        "zerr{p}, sderr{p}, chol_cor{p}, nobs_tg{p}, begin_tg{p}, end_tg{p});\n"
+        "  // compute residual covariance matrix\n",
+        "  Lcortime{p} = cholesky_cor_{cor_fun}({cor_args}, max_nobs_tg{p});\n"
       )
+      if (has_ac_latent_residuals) {
+        str_add(out$tpar_comp) <- glue(
+          "  // compute correlated time-series residuals\n",
+          "  err{p} = scale_time_err(",
+          "zerr{p}, sderr{p}, Lcortime{p}, nobs_tg{p}, begin_tg{p}, end_tg{p});\n"
+        )
+      }
     }
+
   }
 
   acef_sar <- subset2(acef, class = "sar")
@@ -2107,7 +2149,7 @@ stan_dpar_transform <- function(bterms, prior, threads, normalize, ...) {
   }
   if (any(families %in% "logistic_normal")) {
     stopifnot(length(families) == 1L)
-    predcats <- get_predcats(bterms$family)
+    predcats <- make_stan_names(get_predcats(bterms$family))
     sigma_dpars <- glue("sigma{predcats}")
     reqn <- sigma_dpars %in% names(bterms$dpars)
     n <- ifelse(reqn, "[n]", "")
diff --git a/R/summary.R b/R/summary.R
index f2f6b7aa0..883e23635 100644
--- a/R/summary.R
+++ b/R/summary.R
@@ -98,9 +98,9 @@ summary.brmsfit <- function(object, priors = FALSE, prob = 0.95,
 
   variables <- variables(object)
   incl_classes <- c(
-    "b", "bs", "bcs", "bsp", "bmo", "bme", "bmi", "bm", 
-    valid_dpars(object), "delta", "lncor", "rescor", "ar", "ma", 
-    "sderr", "cosy", "lagsar", "errorsar", "car", "sdcar", "rhocar", 
+    "b", "bs", "bcs", "bsp", "bmo", "bme", "bmi", "bm",
+    valid_dpars(object), "delta", "lncor", "rescor", "ar", "ma", "sderr",
+    "cosy", "cortime", "lagsar", "errorsar", "car", "sdcar", "rhocar",
     "sd", "cor", "df", "sds", "sdgp", "lscale", "simo"
   )
   incl_regex <- paste0("^", regex_or(incl_classes), "(_|$|\\[)")
@@ -159,6 +159,13 @@ summary.brmsfit <- function(object, priors = FALSE, prob = 0.95,
   cor_pars <- variables[grepl(regex_autocor_pars(), variables)]
   out$cor_pars <- full_summary[cor_pars, , drop = FALSE]
   rownames(out$cor_pars) <- cor_pars
+  cortime_pars <- variables[grepl("^cortime_", variables)]
+  if (length(cortime_pars)) {
+    tmp <- full_summary[cortime_pars, , drop = FALSE]
+    cortime_pars <- sub("__", ",", sub("__", "(", cortime_pars))
+    rownames(tmp) <- paste0(cortime_pars, ")")
+    out$cor_pars <- rbind(out$cor_pars, tmp)
+  }
 
   # summary of group-level effects
   for (g in out$group) {
diff --git a/inst/chunks/fun_is_equal.stan b/inst/chunks/fun_is_equal.stan
new file mode 100644
index 000000000..1fbff55ac
--- /dev/null
+++ b/inst/chunks/fun_is_equal.stan
@@ -0,0 +1,15 @@
+  // are two 1D integer arrays equal?
+  int is_equal(int[] a, int[] b) {
+    int n_a = size(a);
+    int n_b = size(b);
+    if (n_a != n_b) {
+      return 0;
+    }
+    for (i in 1:n_a) {
+      if (a[i] != b[i]) {
+        return 0;
+      }
+    }
+    return 1;
+  }
+
diff --git a/inst/chunks/fun_normal_time.stan b/inst/chunks/fun_normal_time.stan
index e0663c1aa..cadee5cdf 100644
--- a/inst/chunks/fun_normal_time.stan
+++ b/inst/chunks/fun_normal_time.stan
@@ -1,12 +1,10 @@
   /* multi-normal log-PDF for time-series covariance structures
-   * assuming homogoneous variances
+   * in Cholesky parameterization and assuming homogoneous variances
    * Args:
    *   y: response vector
    *   mu: mean parameter vector
    *   sigma: residual standard deviation
    *   chol_cor: cholesky factor of the correlation matrix
-   *   se2: square of user defined standard errors
-   *     should be set to zero if none are defined
    *   nobs: number of observations in each group
    *   begin: the first observation in each group
    *   end: the last observation in each group
@@ -14,58 +12,130 @@
    *   sum of the log-PDF values of all observations
    */
   real normal_time_hom_lpdf(vector y, vector mu, real sigma, matrix chol_cor,
-                            data vector se2, int[] nobs, int[] begin, int[] end) {
+                            int[] nobs, int[] begin, int[] end) {
     int I = size(nobs);
-    int has_se = max(se2) > 0;
     vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] L = sigma * chol_cor;
     for (i in 1:I) {
-      matrix[nobs[i], nobs[i]] L;
-      L = sigma * chol_cor[1:nobs[i], 1:nobs[i]];
-      if (has_se) {
-        // need to add 'se' to the correlation matrix itself
-        L = multiply_lower_tri_self_transpose(L);
-        L += diag_matrix(se2[begin[i]:end[i]]);
-        L = cholesky_decompose(L);
-      }
+      matrix[nobs[i], nobs[i]] L_i = L[1:nobs[i], 1:nobs[i]];
       lp[i] = multi_normal_cholesky_lpdf(
-        y[begin[i]:end[i]] | mu[begin[i]:end[i]], L
+        y[begin[i]:end[i]] | mu[begin[i]:end[i]], L_i
       );
     }
     return sum(lp);
   }
   /* multi-normal log-PDF for time-series covariance structures
-   * assuming heterogenous variances
-   * Args:
-   *   y: response vector
-   *   mu: mean parameter vector
+   * in Cholesky parameterization and assuming heterogenous variances
+   * Deviating Args:
    *   sigma: residual standard deviation vector
-   *   chol_cor: cholesky factor of the correlation matrix
-   *   se2: square of user defined standard errors
-   *     should be set to zero if none are defined
-   *   nobs: number of observations in each group
-   *   begin: the first observation in each group
-   *   end: the last observation in each group
    * Returns:
    *   sum of the log-PDF values of all observations
    */
   real normal_time_het_lpdf(vector y, vector mu, vector sigma, matrix chol_cor,
-                            data vector se2, int[] nobs, int[] begin, int[] end) {
+                            int[] nobs, int[] begin, int[] end) {
     int I = size(nobs);
-    int has_se = max(se2) > 0;
     vector[I] lp;
     for (i in 1:I) {
-      matrix[nobs[i], nobs[i]] L;
-      L = diag_pre_multiply(sigma[begin[i]:end[i]],
-                            chol_cor[1:nobs[i], 1:nobs[i]]);
-      if (has_se) {
-        // need to add 'se' to the correlation matrix itself
-        L = multiply_lower_tri_self_transpose(L);
-        L += diag_matrix(se2[begin[i]:end[i]]);
-        L = cholesky_decompose(L);
-      }
+      matrix[nobs[i], nobs[i]] L_i;
+      L_i = diag_pre_multiply(sigma[begin[i]:end[i]], chol_cor[1:nobs[i], 1:nobs[i]]);
       lp[i] = multi_normal_cholesky_lpdf(
-        y[begin[i]:end[i]] | mu[begin[i]:end[i]], L
+        y[begin[i]:end[i]] | mu[begin[i]:end[i]], L_i
+      );
+    }
+    return sum(lp);
+  }
+  /* multi-normal log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming homogoneous variances
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real normal_time_hom_flex_lpdf(vector y, vector mu, real sigma, matrix chol_cor,
+                                 int[] nobs, int[] begin, int[] end, int[,] Jtime) {
+    real lp = 0.0;
+    int I = size(nobs);
+    int has_lp[I] = rep_array(0, I);
+    int i = 1;
+    matrix[rows(chol_cor), cols(chol_cor)] L;
+    matrix[rows(chol_cor), cols(chol_cor)] Cov;
+    L = sigma * chol_cor;
+    Cov = multiply_lower_tri_self_transpose(L);
+    while (i <= I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      int lp_terms[I-i+1] = rep_array(0, I-i+1);
+      matrix[nobs[i], nobs[i]] L_i;
+      if (is_equal(iobs, sequence(1, rows(L)))) {
+        // all timepoints are present in this group
+        L_i = L;
+      } else {
+        // arbitrary subsets cannot be taken on L directly
+        L_i = cholesky_decompose(Cov[iobs, iobs]);
+      }
+      has_lp[i] = 1;
+      lp_terms[1] = 1;
+      // find all additional groups where we have the same timepoints
+      for (j in (i+1):I) {
+        if (has_lp[j] == 0 && is_equal(Jtime[j], Jtime[i]) == 1) {
+          has_lp[j] = 1;
+          lp_terms[j-i+1] = 1;
+        }
+      }
+      // vectorize the log likelihood by stacking the vectors
+      lp += multi_normal_cholesky_lpdf(
+        stack_vectors(y, nobs[i], lp_terms, begin[i:I], end[i:I]) |
+        stack_vectors(mu, nobs[i], lp_terms, begin[i:I], end[i:I]), L_i
       );
+      while (i <= I && has_lp[i] == 1) {
+        i += 1;
+      }
+    }
+    return lp;
+  }
+  /* multi-normal log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming heterogenous variances
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   sigma: residual standard deviation vectors
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real normal_time_het_flex_lpdf(vector y, vector mu, vector sigma, matrix chol_cor,
+                                 int[] nobs, int[] begin, int[] end, int[,] Jtime) {
+    int I = size(nobs);
+    vector[I] lp;
+    int has_lp[I] = rep_array(0, I);
+    int i = 1;
+    matrix[rows(chol_cor), cols(chol_cor)] Cor;
+    Cor = multiply_lower_tri_self_transpose(chol_cor);
+    while (i <= I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Lcor_i;
+      matrix[nobs[i], nobs[i]] L_i;
+      if (is_equal(iobs, sequence(1, rows(chol_cor)))) {
+        // all timepoints are present in this group
+        Lcor_i = chol_cor;
+      } else {
+        // arbitrary subsets cannot be taken on chol_cor directly
+        Lcor_i = cholesky_decompose(Cor[iobs, iobs]);
+      }
+      L_i = diag_pre_multiply(sigma[begin[i]:end[i]], Lcor_i);
+      lp[i] = multi_normal_cholesky_lpdf(y[begin[i]:end[i]] | mu[begin[i]:end[i]], L_i);
+      has_lp[i] = 1;
+      // find all additional groups where we have the same timepoints
+      for (j in (i+1):I) {
+        if (has_lp[j] == 0 && is_equal(Jtime[j], Jtime[i]) == 1) {
+          // group j may have different sigmas that group i
+          L_i = diag_pre_multiply(sigma[begin[j]:end[j]], Lcor_i);
+          lp[j] = multi_normal_cholesky_lpdf(y[begin[j]:end[j]] | mu[begin[j]:end[j]], L_i);
+          has_lp[j] = 1;
+        }
+      }
+      while (i <= I && has_lp[i] == 1) {
+        i += 1;
+      }
     }
     return sum(lp);
   }
diff --git a/inst/chunks/fun_normal_time_se.stan b/inst/chunks/fun_normal_time_se.stan
new file mode 100644
index 000000000..ec43de249
--- /dev/null
+++ b/inst/chunks/fun_normal_time_se.stan
@@ -0,0 +1,101 @@
+  /* multi-normal log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming homogoneous variances
+   * and known standard errors
+   * Args:
+   *   y: response vector
+   *   mu: mean parameter vector
+   *   sigma: residual standard deviation
+   *   se2: square of user defined standard errors
+   *   chol_cor: cholesky factor of the correlation matrix
+   *   nobs: number of observations in each group
+   *   begin: the first observation in each group
+   *   end: the last observation in each group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real normal_time_hom_se_lpdf(vector y, vector mu, real sigma, data vector se2,
+                               matrix chol_cor, int[] nobs, int[] begin, int[] end) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cov;
+    Cov = multiply_lower_tri_self_transpose(sigma * chol_cor);
+    for (i in 1:I) {
+      matrix[nobs[i], nobs[i]] Cov_i = Cov[1:nobs[i], 1:nobs[i]];
+      // need to add 'se' to the covariance matrix itself
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_normal_lpdf(y[begin[i]:end[i]] | mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
+  /* multi-normal log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming heterogenous variances
+   * and known standard errors
+   * Deviating Args:
+   *   sigma: residual standard deviation vector
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real normal_time_het_se_lpdf(vector y, vector mu, vector sigma, data vector se2,
+                               matrix chol_cor, int[] nobs, int[] begin, int[] end) {
+    int I = size(nobs);
+    vector[I] lp;
+    for (i in 1:I) {
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = diag_pre_multiply(sigma[begin[i]:end[i]], chol_cor[1:nobs[i], 1:nobs[i]]);
+      // need to add 'se' to the covariance matrix itself
+      Cov_i = multiply_lower_tri_self_transpose(Cov_i);
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_normal_lpdf(y[begin[i]:end[i]] | mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
+  /* multi-normal log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming homogoneous variances
+   * and known standard errors
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real normal_time_hom_se_flex_lpdf(vector y, vector mu, real sigma, data vector se2,
+                                    matrix chol_cor, int[] nobs, int[] begin,
+                                    int[] end, int[,] Jtime) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cov;
+    Cov = multiply_lower_tri_self_transpose(sigma * chol_cor);
+    for (i in 1:I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Cov_i = Cov[iobs, iobs];
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_normal_lpdf(y[begin[i]:end[i]] | mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
+  /* multi-normal log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming heterogenous variances
+   * and known standard errors
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   sigma: residual standard deviation vector
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real normal_time_het_se_flex_lpdf(vector y, vector mu, vector sigma, data vector se2,
+                                    matrix chol_cor, int[] nobs, int[] begin,
+                                    int[] end, int[,] Jtime) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cor;
+    Cor = multiply_lower_tri_self_transpose(chol_cor);
+    for (i in 1:I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = quad_form_diag(Cor[iobs, iobs], sigma[begin[i]:end[i]]);
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_normal_lpdf(y[begin[i]:end[i]] | mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
diff --git a/inst/chunks/fun_scale_time_err.stan b/inst/chunks/fun_scale_time_err.stan
index ecdcfac9a..ab97271d2 100644
--- a/inst/chunks/fun_scale_time_err.stan
+++ b/inst/chunks/fun_scale_time_err.stan
@@ -1,4 +1,5 @@
   /* scale and correlate time-series residuals
+   * using the Cholesky factor of the correlation matrix
    * Args:
    *   zerr: standardized and independent residuals
    *   sderr: standard deviation of the residuals
@@ -13,8 +14,51 @@
                          int[] nobs, int[] begin, int[] end) {
      vector[rows(zerr)] err;
      for (i in 1:size(nobs)) {
-       err[begin[i]:end[i]] =
-         sderr * chol_cor[1:nobs[i], 1:nobs[i]] * zerr[begin[i]:end[i]];
+       matrix[nobs[i], nobs[i]] L_i;
+       L_i = sderr * chol_cor[1:nobs[i], 1:nobs[i]];
+       err[begin[i]:end[i]] = L_i * zerr[begin[i]:end[i]];
      }
      return err;
    }
+  /* scale and correlate time-series residuals
+   * allowx for flexible correlation matrix subsets
+   * Deviating Args:
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   vector of scaled and correlated residuals
+   */
+   vector scale_time_err_flex(vector zerr, real sderr, matrix chol_cor,
+                              int[] nobs, int[] begin, int[] end, int[,] Jtime) {
+     vector[rows(zerr)] err;
+     int I = size(nobs);
+     int has_err[I] = rep_array(0, I);
+     int i = 1;
+     matrix[rows(chol_cor), cols(chol_cor)] L;
+     matrix[rows(chol_cor), cols(chol_cor)] Cov;
+     L = sderr * chol_cor;
+     Cov = multiply_lower_tri_self_transpose(L);
+     while (i <= I) {
+       int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+       matrix[nobs[i], nobs[i]] L_i;
+       if (is_equal(iobs, sequence(1, rows(L)))) {
+         // all timepoints are present in this group
+         L_i = L;
+       } else {
+         // arbitrary subsets cannot be taken on L directly
+         L_i = cholesky_decompose(Cov[iobs, iobs]);
+       }
+       err[begin[i]:end[i]] = L_i * zerr[begin[i]:end[i]];
+       has_err[i] = 1;
+       // find all additional groups where we have the same timepoints
+       for (j in (i+1):I) {
+         if (has_err[j] == 0 && is_equal(Jtime[j], Jtime[i]) == 1) {
+           err[begin[j]:end[j]] = L_i * zerr[begin[j]:end[j]];
+           has_err[j] = 1;
+         }
+       }
+       while (i <= I && has_err[i] == 1) {
+         i += 1;
+       }
+    }
+    return err;
+  }
diff --git a/inst/chunks/fun_stack_vectors.stan b/inst/chunks/fun_stack_vectors.stan
new file mode 100644
index 000000000..1b80de7d0
--- /dev/null
+++ b/inst/chunks/fun_stack_vectors.stan
@@ -0,0 +1,17 @@
+  /* grouped data stored linearly in "data" as indexed by begin and end
+   * is repacked to be stacked into an array of vectors.
+   */
+  vector[] stack_vectors(vector long_data, int n, int[] stack,
+                         int[] begin, int[] end) {
+    int S = sum(stack);
+    int G = size(stack);
+    vector[n] stacked[S];
+    int j = 1;
+    for (i in 1:G) {
+      if (stack[i] == 1) {
+        stacked[j] = long_data[begin[i]:end[i]];
+        j += 1;
+      }
+    }
+    return stacked;
+  }
diff --git a/inst/chunks/fun_student_t_time.stan b/inst/chunks/fun_student_t_time.stan
index d81589d50..089222fde 100644
--- a/inst/chunks/fun_student_t_time.stan
+++ b/inst/chunks/fun_student_t_time.stan
@@ -1,13 +1,11 @@
   /* multi-student-t log-PDF for time-series covariance structures
-   * assuming homogoneous variances
+   * in Cholesky parameterization and assuming homogoneous variances
    * Args:
    *   y: response vector
    *   nu: degrees of freedom parameter
    *   mu: mean parameter vector
    *   sigma: scale parameter
    *   chol_cor: cholesky factor of the correlation matrix
-   *   se2: square of user defined standard errors
-   *     should be set to zero if none are defined
    *   nobs: number of observations in each group
    *   begin: the first observation in each group
    *   end: the last observation in each group
@@ -15,56 +13,88 @@
    *   sum of the log-PDF values of all observations
    */
   real student_t_time_hom_lpdf(vector y, real nu, vector mu, real sigma,
-                               matrix chol_cor, data vector se2, int[] nobs,
-                               int[] begin, int[] end) {
+                               matrix chol_cor, int[] nobs, int[] begin,
+                               int[] end) {
     int I = size(nobs);
-    int has_se = max(se2) > 0;
     vector[I] lp;
     for (i in 1:I) {
-      matrix[nobs[i], nobs[i]] Cov;
-      Cov = sigma * chol_cor[1:nobs[i], 1:nobs[i]];
-      Cov = multiply_lower_tri_self_transpose(Cov);
-      if (has_se) {
-        Cov += diag_matrix(se2[begin[i]:end[i]]);
-      }
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = sigma * chol_cor[1:nobs[i], 1:nobs[i]];
+      Cov_i = multiply_lower_tri_self_transpose(Cov_i);
       lp[i] = multi_student_t_lpdf(
-        y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov
+        y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i
       );
     }
     return sum(lp);
   }
   /* multi-student-t log-PDF for time-series covariance structures
-   * assuming heterogenous variances
-   * Args:
-   *   y: response vector
-   *   nu: degrees of freedom parameter
-   *   mu: mean parameter vector
-   *   sigma: scale parameter vector
-   *   chol_cor: cholesky factor of the correlation matrix
-   *   se2: square of user defined standard errors
-   *     should be set to zero if none are defined
-   *   nobs: number of observations in each group
-   *   begin: the first observation in each group
-   *   end: the last observation in each group
+   * in Cholesky parameterization and assuming heterogenous variances
+   * Deviating Args:
+   *   sigma: residual scale vector
    * Returns:
    *   sum of the log-PDF values of all observations
    */
   real student_t_time_het_lpdf(vector y, real nu, vector mu, vector sigma,
-                               matrix chol_cor, data vector se2, int[] nobs,
-                               int[] begin, int[] end) {
+                               matrix chol_cor, int[] nobs, int[] begin,
+                               int[] end) {
+    int I = size(nobs);
+    vector[I] lp;
+    for (i in 1:I) {
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = diag_pre_multiply(sigma[begin[i]:end[i]], chol_cor[1:nobs[i], 1:nobs[i]]);
+      Cov_i = multiply_lower_tri_self_transpose(Cov_i);
+      lp[i] = multi_student_t_lpdf(
+        y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i
+      );
+    }
+    return sum(lp);
+  }
+  /* multi-student-t log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming homogoneous variances
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real student_t_time_hom_flex_lpdf(vector y, real nu, vector mu, real sigma,
+                                    matrix chol_cor, int[] nobs, int[] begin,
+                                    int[] end, int[,] Jtime) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cov;
+    Cov = multiply_lower_tri_self_transpose(sigma * chol_cor);
+    for (i in 1:I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Cov_i = Cov[iobs, iobs];
+      lp[i] = multi_student_t_lpdf(
+        y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i
+      );
+    }
+    return sum(lp);
+  }
+  /* multi-student-t log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming heterogenous variances
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   sigma: scale parameter vector
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real student_t_time_het_flex_lpdf(vector y, real nu, vector mu, vector sigma,
+                                    matrix chol_cor, int[] nobs, int[] begin,
+                                    int[] end, int[,] Jtime) {
     int I = size(nobs);
-    int has_se = max(se2) > 0;
     vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cor;
+    Cor = multiply_lower_tri_self_transpose(chol_cor);
     for (i in 1:I) {
-      matrix[nobs[i], nobs[i]] Cov;
-      Cov = diag_pre_multiply(sigma[begin[i]:end[i]],
-                              chol_cor[1:nobs[i], 1:nobs[i]]);
-      Cov = multiply_lower_tri_self_transpose(Cov);
-      if (has_se) {
-        Cov += diag_matrix(se2[begin[i]:end[i]]);
-      }
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = quad_form_diag(Cor[iobs, iobs], sigma[begin[i]:end[i]]);
       lp[i] = multi_student_t_lpdf(
-        y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov
+        y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i
       );
     }
     return sum(lp);
diff --git a/inst/chunks/fun_student_t_time_se.stan b/inst/chunks/fun_student_t_time_se.stan
new file mode 100644
index 000000000..6568e5265
--- /dev/null
+++ b/inst/chunks/fun_student_t_time_se.stan
@@ -0,0 +1,103 @@
+  /* multi-student-t log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming homogoneous variances
+   * and known standard errors
+   * Args:
+   *   y: response vector
+   *   nu: degrees of freedom parameter
+   *   mu: mean parameter vector
+   *   sigma: scale parameter
+   *   se2: square of user defined standard errors
+   *   chol_cor: cholesky factor of the correlation matrix
+   *   nobs: number of observations in each group
+   *   begin: the first observation in each group
+   *   end: the last observation in each group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real student_t_time_hom_se_lpdf(vector y, real nu, vector mu, real sigma,
+                                  data vector se2, matrix chol_cor, int[] nobs,
+                                  int[] begin, int[] end) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cov;
+    Cov = multiply_lower_tri_self_transpose(sigma * chol_cor);
+    for (i in 1:I) {
+      matrix[nobs[i], nobs[i]] Cov_i = Cov[1:nobs[i], 1:nobs[i]];
+      // need to add 'se' to the covariance matrix itself
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_student_t_lpdf(y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
+  /* multi-student-t log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming heterogenous variances
+   * and known standard errors
+   * Deviating Args:
+   *   sigma: scale parameter vector
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real student_t_time_het_se_lpdf(vector y, real nu, vector mu, vector sigma,
+                                  data vector se2, matrix chol_cor, int[] nobs,
+                                  int[] begin, int[] end) {
+    int I = size(nobs);
+    vector[I] lp;
+    for (i in 1:I) {
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = diag_pre_multiply(sigma[begin[i]:end[i]], chol_cor[1:nobs[i], 1:nobs[i]]);
+      Cov_i = multiply_lower_tri_self_transpose(Cov_i);
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_student_t_lpdf(y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
+  /* multi-student-t log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming homogoneous variances
+   * and known standard errors
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real student_t_time_hom_se_flex_lpdf(vector y, real nu, vector mu, real sigma,
+                                       data vector se2, matrix chol_cor, int[] nobs,
+                                       int[] begin, int[] end, int[,] Jtime) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cov;
+    Cov = multiply_lower_tri_self_transpose(sigma * chol_cor);
+    for (i in 1:I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Cov_i = Cov[iobs, iobs];
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_student_t_lpdf(y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
+  /* multi-student-t log-PDF for time-series covariance structures
+   * in Cholesky parameterization and assuming heterogenous variances
+   * and known standard errors
+   * allows for flexible correlation matrix subsets
+   * Deviating Args:
+   *   sigma: scale parameter vector
+   *   Jtime: array of time indices per group
+   * Returns:
+   *   sum of the log-PDF values of all observations
+   */
+  real student_t_time_het_se_flex_lpdf(vector y, real nu, vector mu, vector sigma,
+                                       data vector se2, matrix chol_cor, int[] nobs,
+                                       int[] begin, int[] end, int[,] Jtime) {
+    int I = size(nobs);
+    vector[I] lp;
+    matrix[rows(chol_cor), cols(chol_cor)] Cor;
+    Cor = multiply_lower_tri_self_transpose(chol_cor);
+    for (i in 1:I) {
+      int iobs[nobs[i]] = Jtime[i, 1:nobs[i]];
+      matrix[nobs[i], nobs[i]] Cov_i;
+      Cov_i = quad_form_diag(Cor[iobs, iobs], sigma[begin[i]:end[i]]);
+      Cov_i += diag_matrix(se2[begin[i]:end[i]]);
+      lp[i] = multi_student_t_lpdf(y[begin[i]:end[i]] | nu, mu[begin[i]:end[i]], Cov_i);
+    }
+    return sum(lp);
+  }
diff --git a/man/autocor-terms.Rd b/man/autocor-terms.Rd
index 2eeccea78..288613761 100644
--- a/man/autocor-terms.Rd
+++ b/man/autocor-terms.Rd
@@ -6,11 +6,11 @@
 \description{
 Specify autocorrelation terms in \pkg{brms} models. Currently supported terms
 are \code{\link{arma}}, \code{\link{ar}}, \code{\link{ma}},
-\code{\link{cosy}}, \code{\link{sar}}, \code{\link{car}}, and
-\code{\link{fcor}}. Terms can be directly specified within the formula, or
-passed to the \code{autocor} argument of \code{\link{brmsformula}} in the
-form of a one-sided formula. For deprecated ways of specifying
-autocorrelation terms, see \code{\link{cor_brms}}.
+\code{\link{cosy}}, \code{\link{unstr}}, \code{\link{sar}},
+\code{\link{car}}, and \code{\link{fcor}}. Terms can be directly specified
+within the formula, or passed to the \code{autocor} argument of
+\code{\link{brmsformula}} in the form of a one-sided formula. For deprecated
+ways of specifying autocorrelation terms, see \code{\link{cor_brms}}.
 }
 \details{
 The autocor term functions are almost solely useful when called in
@@ -31,6 +31,6 @@ bf(y ~ x) + acformula(~ arma(p = 1, q = 1) + car(M))
 \seealso{
 \code{\link{brmsformula}}, \code{\link{acformula}},
   \code{\link{arma}}, \code{\link{ar}}, \code{\link{ma}},
-  \code{\link{cosy}}, \code{\link{sar}}, \code{\link{car}},
-  \code{\link{fcor}}
+  \code{\link{cosy}}, \code{\link{unstr}}, \code{\link{sar}},
+  \code{\link{car}}, \code{\link{fcor}}
 }
diff --git a/man/unstr.Rd b/man/unstr.Rd
new file mode 100644
index 000000000..935255641
--- /dev/null
+++ b/man/unstr.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/formula-ac.R
+\name{unstr}
+\alias{unstr}
+\title{Set up UNSTR correlation structures}
+\usage{
+unstr(time, gr)
+}
+\arguments{
+\item{time}{An optional time variable specifying the time ordering
+of the observations. By default, the existing order of the observations
+in the data is used.}
+
+\item{gr}{An optional grouping variable. If specified, the correlation
+structure is assumed to apply only to observations within the same grouping
+level.}
+}
+\value{
+An object of class \code{'unstr_term'}, which is a list
+  of arguments to be interpreted by the formula
+  parsing functions of \pkg{brms}.
+}
+\description{
+Set up an unstructured (UNSTR) correlation term in \pkg{brms}. The function does
+not evaluate its arguments -- it exists purely to help set up a model with
+UNSTR terms.
+}
+\examples{
+\dontrun{
+# add an unstructured correlation matrix for visits within the same patient
+fit <- brm(count ~ Trt + unstr(visit, patient), data = epilepsy)
+summary(fit)
+}
+
+}
+\seealso{
+\code{\link{autocor-terms}}
+}
diff --git a/tests/local/tests.models-4.R b/tests/local/tests.models-4.R
index b364d8d24..7d044e427 100644
--- a/tests/local/tests.models-4.R
+++ b/tests/local/tests.models-4.R
@@ -1,5 +1,29 @@
 source("setup_tests_local.R")
 
+test_that("UNSTR models work correctly", {
+  epilepsy2 <- epilepsy
+  epilepsy2$visit <- as.numeric(epilepsy2$visit)
+  fit <- brm(count ~ Trt + unstr(visit, patient), data = epilepsy2)
+  print(fit)
+  expect_ggplot(pp_check(fit))
+
+  waic <- waic(fit)
+  expect_range(waic$estimates[3, 1], 1550, 1600)
+  # ensure the the correlation are actually included in the predictions
+  waic_without <- waic(fit, incl_autocor = FALSE)
+  expect_true(waic$estimates[3, 1] + 200 < waic_without$estimates[3, 1])
+
+  waic_new <- waic(fit, newdata = epilepsy2[1:100, ])
+  expect_range(waic_new $estimates[3, 1], 700, 760)
+
+  newdat <- epilepsy2[1:5, ]
+  newdat$visit[1] <- 5
+  expect_error(
+    waic(fit, newdata = newdat),
+    "Cannot handle new time points in UNSTR models"
+  )
+})
+
 test_that("SAR models work correctly", {
   data(oldcol, package = "spdep")
   fit_lagsar <- brm(CRIME ~ INC + HOVAL + sar(COL.nb),
diff --git a/tests/testthat/tests.make_stancode.R b/tests/testthat/tests.make_stancode.R
index 0f2daf593..d854a2dc6 100644
--- a/tests/testthat/tests.make_stancode.R
+++ b/tests/testthat/tests.make_stancode.R
@@ -761,7 +761,7 @@ test_that("Stan code for ARMA models is correct", {
 
   bform <- bf(y ~ x, sigma ~ x) + acformula(~arma(time, cov = TRUE))
   scode <- make_stancode(bform, dat, family = student)
-  expect_match2(scode, "student_t_time_het_lpdf(Y | nu, mu, sigma, chol_cor")
+  expect_match2(scode, "student_t_time_het_lpdf(Y | nu, mu, sigma, Lcortime")
 
   bform <- bf(y ~ exp(eta) - 1, eta ~ x, autocor = ~ar(time), nl = TRUE)
   scode <- make_stancode(bform, dat, family = student,
@@ -773,9 +773,9 @@ test_that("Stan code for ARMA models is correct", {
     y ~ x + ar(time, cov = TRUE), dat, family = poisson,
     prior = prior(cauchy(0, 10), class = sderr)
   )
-  expect_match2(scode, "chol_cor = cholesky_cor_ar1(ar[1], max_nobs_tg);")
+  expect_match2(scode, "Lcortime = cholesky_cor_ar1(ar[1], max_nobs_tg);")
   expect_match2(scode,
-    "err = scale_time_err(zerr, sderr, chol_cor, nobs_tg, begin_tg, end_tg);"
+    "err = scale_time_err(zerr, sderr, Lcortime, nobs_tg, begin_tg, end_tg);"
   )
   expect_match2(scode, "mu += Intercept + Xc * b + err;")
   expect_match2(scode, "lprior += cauchy_lpdf(sderr | 0, 10)")
@@ -798,14 +798,53 @@ test_that("Stan code for compound symmetry models is correct", {
     prior = prior(normal(0, 2), cosy)
   )
   expect_match2(scode, "real<lower=0,upper=1> cosy;")
-  expect_match2(scode, "chol_cor = cholesky_cor_cosy(cosy, max_nobs_tg);")
+  expect_match2(scode, "Lcortime = cholesky_cor_cosy(cosy, max_nobs_tg);")
   expect_match2(scode, "lprior += normal_lpdf(cosy | 0, 2)")
 
   scode <- make_stancode(bf(y ~ x + cosy(time), sigma ~ x), dat)
-  expect_match2(scode, "normal_time_het_lpdf(Y | mu, sigma, chol_cor")
+  expect_match2(scode, "normal_time_het_lpdf(Y | mu, sigma, Lcortime")
 
   scode <- make_stancode(y ~ x + cosy(time), dat, family = poisson)
-  expect_match2(scode, "chol_cor = cholesky_cor_cosy(cosy, max_nobs_tg);")
+  expect_match2(scode, "Lcortime = cholesky_cor_cosy(cosy, max_nobs_tg);")
+})
+
+test_that("Stan code for UNSTR covariance terms is correct", {
+  dat <- data.frame(y = 1:12, x = rnorm(12), tim = c(5:1, 1:5, c(0, 4)),
+                    g = c(rep(3:4, 5), rep(2, 2)))
+
+  scode <- make_stancode(y ~ x + unstr(tim, g), data = dat)
+  expect_match2(scode, "normal_time_hom_flex_lpdf(Y | mu, sigma, Lcortime, nobs_tg, begin_tg, end_tg, Jtime_tg);")
+  expect_match2(scode, "cortime[choose(k - 1, 2) + j] = Cortime[j, k];")
+  expect_match2(scode, "lprior += lkj_corr_cholesky_lpdf(Lcortime | 1);")
+
+  scode <- make_stancode(
+    y ~ x + unstr(tim, g), data = dat,
+    family = student(), prior = prior(lkj(4), cortime)
+  )
+  expect_match2(scode, "student_t_time_hom_flex_lpdf(Y | nu, mu, sigma, Lcortime, nobs_tg, begin_tg, end_tg, Jtime_tg);")
+  expect_match2(scode, "lprior += lkj_corr_cholesky_lpdf(Lcortime | 4);")
+
+  # test standard error
+  scode <- make_stancode(
+    y | se(1, sigma = TRUE) ~ x + unstr(tim, g),
+    data = dat, family = gaussian(),
+  )
+  expect_match2(scode, "normal_time_hom_se_flex_lpdf(Y | mu, sigma, se2, Lcortime, nobs_tg, begin_tg, end_tg, Jtime_tg);")
+
+  # test latent representation
+  scode <- make_stancode(
+    y ~ x + unstr(tim, g), data = dat,
+    family = poisson()
+  )
+  expect_match2(scode, "err = scale_time_err_flex(zerr, sderr, Lcortime, nobs_tg, begin_tg, end_tg,")
+  expect_match2(scode, "mu += Intercept + Xc * b + err;")
+
+  # non-linear model
+  scode <- make_stancode(
+    bf(y ~ a, a ~ x, autocor = ~ unstr(tim, g), nl = TRUE),
+    data = dat, family = student(), prior = prior(normal(0,1), nlpar = a)
+  )
+  expect_match2(scode, "student_t_time_hom_flex_lpdf(Y | nu, mu, sigma, Lcortime, nobs_tg, begin_tg, end_tg, Jtime_tg);")
 })
 
 test_that("Stan code for intercept only models is correct", {
diff --git a/tests/testthat/tests.make_standata.R b/tests/testthat/tests.make_standata.R
index 3915f07d7..feac574ae 100644
--- a/tests/testthat/tests.make_standata.R
+++ b/tests/testthat/tests.make_standata.R
@@ -287,6 +287,17 @@ test_that("make_standata returns correct data for ARMA terms", {
   sdata <- make_standata(bform, dat)
 })
 
+test_that("make_standata returns correct data for UNSTR covariance terms", {
+  dat <- data.frame(y = 1:12, x = rnorm(12), tim = c(5:1, 1:5, c(0, 4)),
+                    g = c(rep(3:4, 5), rep(2, 2)))
+
+  sdata <- make_standata(y ~ x + unstr(tim, g), data = dat)
+  expect_equal(sdata$n_unique_t, 6)
+  expect_equal(sdata$n_unique_cortime, 15)
+  Jtime <- rbind(c(1, 5, 0, 0, 0), 2:6, 2:6)
+  expect_equal(sdata$Jtime_tg, Jtime)
+})
+
 test_that("make_standata allows to retrieve the initial data order", {
   dat <- data.frame(y1 = rnorm(100), y2 = rnorm(100),
                           id = sample(1:10, 100, TRUE),