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'

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)
  • 1
    $\begingroup$ 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. 👍 $\endgroup$ Nov 20, 2019 at 16:19

1 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:




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:


  • 1
    $\begingroup$ Both sources assume that the factors evolve as a Guassian random walk. I'm having trouble altering the code of the sim_dfa function so that the factors evolve as an VAR(1) process, but I keep ending up with NaNs. $\endgroup$ Nov 18, 2019 at 7:54
  • 1
    $\begingroup$ Would something like this help: rdocumentation.org/packages/tsDyn/versions/0.8-1/topics/VAR.sim? $\endgroup$ Nov 18, 2019 at 16:26
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 18, 2019 at 20:41
  • 1
    $\begingroup$ Most definitely. Will add it once I obtain something that runs successfully. $\endgroup$ Nov 18, 2019 at 22:28
  • 1
    $\begingroup$ I fixed the problem. 1) I should just fix the model parameter values instead of simulating them according to the prior distributions I mentioned in the OP. 2) The documentation has a misprint in the arguments of the 'VAR.sim' function. $\endgroup$ Nov 20, 2019 at 3:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.