How can you draw samples from a multidimensional time series? I'm aware of MCMC methods such as Metropolis Hastings and friends, but these methods assume a stable posterior distribution. Is there a way to draw samples from a multidimensional timeseries? For example, the purchases on a debit card. There is correlation and seasonality in purchases: buy coffee in the morning, lunch in the afternoon, dinner at night. These purchases are heavily correlated - when you buy coffee (\$2) you buy it at a coffee shop, when you buy groceries (\$100) you buy it in a supermarket, and so on. 
Edit for specificity:
For example, given a dataset of N observations, like debit card transactions:
| time | merchant | purchase | cost |
-------------------------------------
|17:00 |  shop A  |   food   | $10  |
|09:30 |  shop B  |  coffee  | $2   |
| ...  |   ...    |    ...   | ...  | 

Time is a continuous variable, merchant and purchase are categorical and cost is discrete. 
Is there a method for drawing samples from these sorts of distributions, with or without making model assumptions (I'm sure there is)? I mention MCMC because it's the closest thing I know to this, but I don't know if it is the solution!
 A: You can calculate some statistic of a given multi-variate time series using boot::tsboot. According to the documentation for boot:

The replicate time series can be generated using fixed or random block lengths or can be model-based replicates.

library(tidyverse)
library(magrittr)
library(mvtnorm)
library(boot)

set.seed(123)
sigma <- diag(c(2,3))    
sigma[1,2] <- sigma[2,1] <- 0.8
n <- 200
eps <- mvtnorm::rmvnorm(n = n, mean = rep(0, 2), sigma = sigma)

## Simulate two correlated AR(1) process
xy <- data.frame(x = arima.sim(n = n, model = list(ar = c(0.7)), innov = 
  eps[,1]), y = arima.sim(n = n, model = list(ar = c(-0.4)), innov = eps[,2]))

## Examine the correlation
cor_xy <- function(xy){cor(xy$x, xy$y)}
xy %>% cor_xy()

## Examine the correlation using bootstrap
(b = boot::tsboot(xy, cor_xy, R = 1000, l = 20, sim = "fixed"))
plot(b)

BLOCK BOOTSTRAP FOR TIME SERIES

Fixed Block Length of 20 

Call:
boot::tsboot(tseries = xy, statistic = cor_xy, R = 1000, l = 20, 
    sim = "fixed")


Bootstrap Statistics :
      original       bias    std. error
t1* 0.07804724 -0.003418753  0.05543288

The bootstrap results can be visualized immediately using plot().

Discussion about the general bootstrap method can be found in How do you do bootstrapping with time-series data?. See Understanding the output of a bootstrap performed in R for detailed interpretation of output in boot::tsboot.
The way I simulate the multivariate time series is from How do you simulate two correlated AR(p) time series?.
For now, I don't know how to extract the samples from the bootstrap, which might be important in some cases. I am not sure if the algorithm is similar to that for univariate time series. Also, I don't know if the method applies to a time series and a categorical series.
