How can I simulate observations from a dynamic factor model?

Question

Consider the following dynamic factor model where for $t = 1, 2, \ldots, T$ \begin{align} x_t &= \Lambda f_t + e_t \quad \text{ where } e_t = \Phi_1 e_{t-1} + \epsilon_t \text{ and } \epsilon_t \overset{\text{iid}}{\sim} \mathcal{N} (0, \Sigma_x) \\ f_t &= \Psi_1 f_{t-1} + \varepsilon_t \quad \text{ where } \varepsilon_t \overset{\text{iid}}{\sim} \mathcal{N} (0, \Sigma_f) \end{align}

The observations $x_t$ are vectors with length $N$. The latent dynamic factors $f_t$ are vectors with length $K$. The matrix of factor loadings $\Lambda$ is of size $N \times K$ and the first $K \times K$ square part is lower triangular with 1's along the diagonal. The matrices of autoregressive coefficients $\Phi_1$ and $\Psi$ are diagonal and of size $N \times N$ and $K \times K$, respectively. The covariance matrices $\Sigma_x$ and $\Sigma_f$ are diagonal and of size $N \times N$ and $K \times K$, respectively. More precisely, \begin{align} \Phi_1 &= \text{diag} (\phi_{11}, \dots, \phi_{1N}) \\ \Psi_1 &= \text{diag} (\psi_{11}, \dots, \psi_{1K}) \\ \Sigma_x &= \text{diag} (\sigma_{x1}^2, \ldots, \sigma_{xN}^2) \\ \Sigma_f &= \text{diag} (\sigma_{f1}^2, \ldots, \sigma_{fK}^2) \\ \Lambda &= \begin{bmatrix} 1 & 0 \\ \lambda_{21} & 1 \\ \lambda_{31} & \lambda_{32} \\ \vdots & \vdots \\ \lambda_{N1} & \lambda_{N2} \end{bmatrix} \text{ for example if } K = 2 \end{align}

The prior distributions are \begin{align} \phi_{1i} &\overset{iid}{\sim} \mathcal{N} (0,1) \text{ for } i = 1,\ldots, N \\ \psi_{1j} &\overset{iid}{\sim} \mathcal{N} (0,1) \text{ for } j = 1,\ldots, K \\ \lambda_{nk} &\overset{iid}{\sim} \mathcal{N} (0,1) \text{ for the lower triangular part of } \Lambda \\ \sigma_{xi}^2 &\overset{iid}{\sim} \text{inv-}\chi^2 (\nu = 4, \tau^2 = 0.01) \text{ for } i = 1,\ldots, N \\ \sigma_{fj}^2 &\overset{iid}{\sim} \text{inv-}\chi^2 (\nu = 4, \tau^2 = 0.01) \text{ for } j = 1,\ldots, K \end{align}

Question: How can I simulate/generate observations according to this model? That is, how can I obtain $x_1, x_2, \ldots, x_T$?

My thoughts: I am new to simulations in general. The process described below seems to be the most obvious. While I believe my thought process to be logical, I don't really have a way to make sure that it is correct as most texts involving models of this nature use real datasets.

1. Simulate model parameters according to their prior distributions.
2. Simulate $\varepsilon_t$ from $\mathcal{N} (0, \Sigma_f)$ for $t = 1, \ldots, T$.
3. Simulate $f_0$ from $\mathcal{N} (0, A_f)$ where $A_f$ comes from some fixed point equation (idea comes from Bayesian Core).
4. Simulate $f_t$ from $\mathcal{N} (\Psi_1 f_{t-1}, \Sigma_f)$ for $t = 1, \ldots, T$
5. Simulate $e_0$ from $\mathcal{N} (0, A_x)$ where $A_x$ comes from some fixed point equation.
6. Simulate $e_t$ from $\mathcal{N} (\Phi_1 e_{t-1}, \Sigma_f)$ for $t = 1, \ldots, T$.
7. Obtain $x_t$ using $\Lambda f_t + e_t$.

---

Based on the help from @IsabellaGhement, I came up with the following code that I think does what I am looking for.

library(tsDyn) #for the function 'VAR.sim'
library(lgarch) #for the function 'rmnorm' that is used in the function 'VAR.sim'

set.seed(3)
N <- 250 #number of variables
TT <- 150 #number of time series observations
K <- 5 #number of hidden factors

# specify all the model parameters
sigma.X <- phi <- 0.5 * diag(N)
sigma.F <- psi <- 0.3 * diag(K)
lambda <- matrix( rep( 0, times = N * K ), nrow = N, ncol = K )
lambda[1:K, 1:K] <- diag(K) #turns the first square into the identity matrix
lambda[lower.tri(lambda)] <- rnorm( K * (K - 1) / 2 + (N - K) * K ) #fills lower triangular parts with values from the normal distribution

#generate idiosyncratic components, e_t
e <- VAR.sim( B = phi, n = TT, include = "none", 
            innov = rmnorm(TT, mean = rep(0, TT), vcov = sigma.X), 
            varcov = sigma.X )

#generate the latent factors, f_t
f <- VAR.sim( B = psi, n = TT, include = "none", 
            innov = rmnorm(TT, mean = rep(0, TT), vcov = sigma.F), 
            varcov = sigma.F )

#generate observations, x_t
obs <- lambda %*% t(f) + t(e)

Answer 1

The package bayesdfa in R has a function called sim_dfa() which simulates data from a dynamic factor model, so you can look into its bowels to get some ideas for how you can simulate data for this type of model:

install.packages("bayesdfa")

library(bayesdfa)

set.seed(1)

sim_dat <- sim_dfa(
   num_trends = 2, 
   num_years = 20,
   num_ts = 4
)

Just use the command edit(sim_dfa) to access the function definition. Also, see the bayesdfa package vignette: https://cran.r-project.org/web/packages/bayesdfa/vignettes/bayesdfa.html.

Another detailed example of simulating data for this type of model can be found in this post on Dynamic Factor Analysis with the greta package for R:

https://mdscheuerell.github.io/gretaDFA/