# How can I simulate observations from a dynamic factor model?

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)

• Cool! I will play with this in the next few days to build some intuition for what it does. Thank you so much for sharing your code. 👍 Nov 20, 2019 at 16:19

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:

• I think it does! I can simulate $f_t$ for $t = 1, 2, \ldots, T$ using 'VAR.sim' and do the same thing to simulate $e_t$ for $t = 1, 2, \ldots, T$. Then use the results to obtain the observations. Nov 18, 2019 at 20:41