diff --git a/R/forecast.R b/R/forecast.R
index 338177e8..526f77a3 100644
--- a/R/forecast.R
+++ b/R/forecast.R
@@ -162,8 +162,14 @@ forecast.mdl_ts <- function(object, new_data = NULL, h = NULL, bias_adjust = NUL
   # Compute forecasts
   if(simulate || bootstrap) {
     fc <- generate(object, new_data, bootstrap = bootstrap, times = times, ...)
-    fc <- unname(split(object$transformation[[1]](fc[[".sim"]]), fc[[index_var(fc)]]))
-    fc <- distributional::dist_sample(fc)
+    fc_idx <- fc[[index_var(fc)]]
+    fc <- if (length(resp_vars) > 1) {
+      do.call(cbind, fc[resp_vars])
+    } else {
+      fc[[".sim"]]
+    }
+    
+    fc <- distributional::dist_sample(vctrs::vec_split(fc, fc_idx)$val)
   } else {
     # Compute specials with new_data
     object$model$stage <- "forecast"
@@ -180,54 +186,46 @@ forecast.mdl_ts <- function(object, new_data = NULL, h = NULL, bias_adjust = NUL
     object$model$remove_data()
     object$model$stage <- NULL
     fc <- forecast(object$fit, new_data, specials = specials, times = times, ...)
-  }
-  
-  # Back-transform forecast distributions
-  bt <- map(object$transformation, function(x){
-    trans <- x%@%"inverse"
-    inv_trans <- `attributes<-`(x, NULL)
-    req_vars <- setdiff(all.vars(body(trans)), names(formals(trans)))
-    if(any(req_vars %in% names(new_data))) {
-      trans <- lapply(
-        vec_chop(new_data[req_vars]),
-        function(transform_data) {
-          set_env(trans, new_environment(transform_data, get_env(trans)))
-        }
-      )
-      attr(trans, "inverse") <- lapply(
-        vec_chop(new_data[req_vars]),
-        function(transform_data) {
-          set_env(inv_trans, new_environment(transform_data, get_env(inv_trans)))
-        }
-      )
-      trans
-    } else {
-      structure(list(trans), inverse = list(inv_trans))
+    
+    # Back-transform forecast distributions
+    bt <- map(object$transformation, function(x){
+      trans <- x%@%"inverse"
+      inv_trans <- `attributes<-`(x, NULL)
+      req_vars <- setdiff(all.vars(body(trans)), names(formals(trans)))
+      if(any(req_vars %in% names(new_data))) {
+        trans <- lapply(
+          vec_chop(new_data[req_vars]),
+          function(transform_data) {
+            set_env(trans, new_environment(transform_data, get_env(trans)))
+          }
+        )
+        attr(trans, "inverse") <- lapply(
+          vec_chop(new_data[req_vars]),
+          function(transform_data) {
+            set_env(inv_trans, new_environment(transform_data, get_env(inv_trans)))
+          }
+        )
+        trans
+      } else {
+        structure(list(trans), inverse = list(inv_trans))
+      }
+    })
+    
+    is_transformed <- vapply(bt, function(x) !is_symbol(body(x[[1]])), logical(1L))
+    if(length(bt) > 1) {
+      if(any(is_transformed)){
+        abort("Transformations of multivariate forecasts distributions are not supported, use simulate = TRUE or bootstrap = TRUE.")
+      }
     }
-#     exists_vars <- map_lgl(req_vars, exists, env)
-#     if(any(!exists_vars)){
-#       bt <- custom_error(bt, sprintf(
-# "Unable to find all required variables to back-transform the forecasts (missing %s).
-# These required variables can be provided by specifying `new_data`.",
-#         paste0("`", req_vars[!exists_vars], "`", collapse = ", ")
-#       ))
-#     }
-  })
-  
-  is_transformed <- vapply(bt, function(x) !is_symbol(body(x[[1]])), logical(1L))
-  if(length(bt) > 1) {
-    if(any(is_transformed)){
-      abort("Transformations of multivariate forecasts are not yet supported")
-    }
-  }
-  if(any(is_transformed)) {
-    if (identical(unique(dist_types(fc)), "dist_sample")) {
-      fc <- distributional::dist_sample(
-        .mapply(exec, list(bt[[1]], distributional::parameters(fc)$x), MoreArgs = NULL)
-      )
-    } else {
-      bt <- bt[[1]]
-      fc <- distributional::dist_transformed(fc, `attributes<-`(bt, NULL), bt%@%"inverse")
+    if(any(is_transformed)) {
+      if (identical(unique(dist_types(fc)), "dist_sample")) {
+        fc <- distributional::dist_sample(
+          .mapply(exec, list(bt[[1]], distributional::parameters(fc)$x), MoreArgs = NULL)
+        )
+      } else {
+        bt <- bt[[1]]
+        fc <- distributional::dist_transformed(fc, `attributes<-`(bt, NULL), bt%@%"inverse")
+      }
     }
   }
   
diff --git a/R/generate.R b/R/generate.R
index ea52fadc..7d8e6422 100644
--- a/R/generate.R
+++ b/R/generate.R
@@ -7,7 +7,7 @@
 #' Innovations are sampled by the model's assumed error distribution. 
 #' If `bootstrap` is `TRUE`, innovations will be sampled from the model's 
 #' residuals. If `new_data` contains the `.innov` column, those values will be
-#' treated as innovations for the simulated paths..
+#' treated as innovations for the simulated paths.
 #' 
 #' @param x A mable.
 #' @param new_data The data to be generated (time index and exogenous regressors)
@@ -77,20 +77,18 @@ generate.mdl_ts <- function(x, new_data = NULL, h = NULL, times = 1, seed = NULL
   }
   
   if(bootstrap) {
-    if(length(x$response) > 1) abort("Generating bootstrap paths from multivariate models is not yet supported.")
     res <- residuals(x$fit)
-    res <- stats::na.omit(res) - mean(res, na.rm = TRUE)
-    new_data$.innov <- if(bootstrap_block_size == 1) {
-      sample(res, nrow(new_data), replace = TRUE)
+    f_mean <- if(length(x$response) == 1) mean else colMeans
+    res <- stats::na.omit(res) - f_mean(res, na.rm = TRUE)
+    i <- if(bootstrap_block_size == 1) {
+      sample.int(NROW(res), nrow(new_data), replace = TRUE)
     } else {
       if(any(has_gaps(x$data)$.gaps)) abort("Residuals must be regularly spaced without gaps to use a block bootstrap method.")
       kr <- tsibble::key_rows(new_data)
-      # idx <- x$data[[index_var(x$data)]]
-      # new_idx <- new_data[[index_var(new_data)]]
-      # block_pos <- ((new_idx - min(idx))%%bootstrap_block_size)+1
-      innov <- lapply(lengths(kr), function(n) block_bootstrap(res, bootstrap_block_size, size = n))
-      vec_c(!!!innov)
+      ki <- lapply(lengths(kr), function(n) block_bootstrap(NROW(res), bootstrap_block_size, size = n))
+      vec_c(!!!ki)
     }
+    new_data$.innov <- if(length(x$response) == 1) res[i] else res[i,]
   }
   
   # Compute specials with new_data
@@ -124,13 +122,13 @@ Does your model require extra variables to produce simulations?", e$message))
   .sim
 }
 
-block_bootstrap <- function (x, window_size, size = length(x)) {
+block_bootstrap <- function (n, window_size, size = length(x)) {
   n_blocks <- size%/%window_size + 2
   bx <- numeric(n_blocks * window_size)
   for (i in seq_len(n_blocks)) {
-    block_pos <- sample(seq_len(length(x) - window_size + 1), 1)
-    bx[((i - 1) * window_size + 1):(i * window_size)] <- x[block_pos:(block_pos + window_size - 1)]
+    block_pos <- sample(seq_len(n - window_size + 1), 1)
+    bx[((i - 1) * window_size + 1):(i * window_size)] <- block_pos:(block_pos + window_size - 1)
   }
-  start_from <- sample(0:(window_size - 1), 1) + 1
+  start_from <- sample.int(window_size, 1)
   bx[seq(start_from, length.out = size)]
 }
\ No newline at end of file
diff --git a/R/irf.R b/R/irf.R
index 5b1f9d6f..7df78cdf 100644
--- a/R/irf.R
+++ b/R/irf.R
@@ -3,8 +3,10 @@
 #' This function calculates the impulse response function (IRF) of a time series model.
 #' The IRF describes how a model's variables react to external shocks over time.
 #'
+#' If `new_data` contains the `.impulse` column, those values will be
+#' treated as impulses for the calculated impulse responses.
+#'
 #' @param x A fitted model object, such as from a VAR or ARIMA model. This model is used to compute the impulse response.
-#' @param impulse A character string specifying the name of the variable that is shocked (the impulse variable).
 #' @param ... Additional arguments to be passed to lower-level functions.
 #' 
 #' @details
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 6f2db910..de8913c4 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -3,6 +3,7 @@ Blaskowitz
 CRPS
 Heeyoung
 Herwartz
+IRF
 JRSS
 MAAPE
 MASE
@@ -10,6 +11,7 @@ MatrixM
 MinT
 ORCID
 Sungil
+VAR
 Wickramasuriya
 backtransform
 dable
@@ -17,7 +19,6 @@ dables
 doi
 dplyr
 emperical
-env
 erroring
 etc
 forecast's
@@ -38,6 +39,7 @@ seasonalities
 superceded
 tibble
 tidyr
+tidyselect
 tidyverse
 tidyverts
 tsibble
diff --git a/tests/testthat/test-generate.R b/tests/testthat/test-generate.R
index ff773618..ce1ee00e 100644
--- a/tests/testthat/test-generate.R
+++ b/tests/testthat/test-generate.R
@@ -17,11 +17,6 @@ test_that("generate", {
   expect_equal(gen_complex$index, yearmonth("1979 Jan") + rep(0:23, 2*2*3))
   expect_equal(unique(gen_complex$key), c("fdeaths", "mdeaths"))
   expect_equal(unique(gen_complex$.model), c("ets", "lm"))
-  
-  expect_error(
-    mbl_mv %>% generate(),
-    "Generating paths from multivariate models is not yet supported"
-  )
 })
 
 test_that("generate seed setting", {