Fit the Discrete-Time Vector Autoregressive Model

Model

The measurement model is given by $\begin{equation} \mathbf{y}_{i, t} = \boldsymbol{\Lambda} \boldsymbol{\eta}_{i, t} + \boldsymbol{\varepsilon}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\varepsilon}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Theta} \right) \end{equation}$ where $\mathbf{y}_{i, t}$ , $\boldsymbol{\eta}_{i, t}$ , and $\boldsymbol{\varepsilon}_{i, t}$ are random variables and $\boldsymbol{\Lambda}$ , and $\boldsymbol{\Theta}$ are model parameters. $\mathbf{y}_{i, t}$ represents a vector of observed random variables, $\boldsymbol{\eta}_{i, t}$ a vector of latent random variables, and $\boldsymbol{\varepsilon}_{i, t}$ a vector of random measurement errors, at time $t$ and individual $i$ . $\boldsymbol{\Lambda}$ denotes a matrix of factor loadings, and $\boldsymbol{\Theta}$ the covariance matrix of $\boldsymbol{\varepsilon}$ . In this model, $\boldsymbol{\Lambda}$ is an identity matrix and $\boldsymbol{\Theta}$ is a diagonal matrix.

The dynamic structure is given by $\begin{equation} \boldsymbol{\eta}_{i, t} = \boldsymbol{\alpha} + \boldsymbol{\beta} \boldsymbol{\eta}_{i, t - 1} + \boldsymbol{\zeta}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\zeta}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Psi} \right) \end{equation}$ where $\boldsymbol{\eta}_{i, t}$ , $\boldsymbol{\eta}_{i, t - 1}$ , and $\boldsymbol{\zeta}_{i, t}$ are random variables, and $\boldsymbol{\alpha}$ , $\boldsymbol{\beta}$ , and $\boldsymbol{\Psi}$ are model parameters. Here, $\boldsymbol{\eta}_{i, t}$ is a vector of latent variables at time $t$ and individual $i$ , $\boldsymbol{\eta}_{i, t - 1}$ represents a vector of latent variables at time $t - 1$ and individual $i$ , and $\boldsymbol{\zeta}_{i, t}$ represents a vector of dynamic noise at time $t$ and individual $i$ . $\boldsymbol{\alpha}$ denotes a vector of intercepts, $\boldsymbol{\beta}$ a matrix of autoregression and cross regression coefficients, and $\boldsymbol{\Psi}$ the covariance matrix of $\boldsymbol{\zeta}_{i, t}$ . In this model, $\boldsymbol{\Psi}$ is a symmetric matrix.

Data Generation

The parameters used in this example were based on Bringmann et al. (2013).

Notation

Let $t = 1000$ be the number of time points and $n = 1000$ be the number of individuals.

Let the measurement model intecept vector $\boldsymbol{\nu}$ be given by

$\begin{equation} \boldsymbol{\nu} = \left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right) . \end{equation}$

Let the factor loadings matrix $\boldsymbol{\Lambda}$ be given by

$\begin{equation} \boldsymbol{\Lambda} = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right) . \end{equation}$

Let the measurement error covariance matrix $\boldsymbol{\Theta}$ be given by

$\begin{equation} \boldsymbol{\Theta} = \left( \begin{array}{cc} 0.2 & 0 \\ 0 & 0.2 \\ \end{array} \right) . \end{equation}$

Let the initial condition $\boldsymbol{\eta}_{0}$ be given by

$\begin{equation} \boldsymbol{\eta}_{0} \sim \mathcal{N} \left( \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0}, \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} \right) \end{equation}$

$\begin{equation} \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{c} 3.8143546 \\ 2.5763481 \\ \end{array} \right) \end{equation}$

$\begin{equation} \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{cc} 1.3978842 & 0.5782369 \\ 0.5782369 & 1.6636513 \\ \end{array} \right) . \end{equation}$

Let the transition matrix $\boldsymbol{\beta}$ be given by

$\begin{equation} \left( \begin{array}{cc} 0.28 & -0.048 \\ -0.035 & 0.26 \\ \end{array} \right) . \end{equation}$

Let the dynamic process noise $\boldsymbol{\Psi}$ be given by

$\begin{equation} \boldsymbol{\Psi} = \left( \begin{array}{cc} 1.3 & 0.5696315 \\ 0.5696315 & 1.56 \\ \end{array} \right) . \end{equation}$

R Function Arguments

n
#> [1] 1000
time
#> [1] 1000
mu0
#> [1] 3.814355 2.576348
sigma0
#>           [,1]      [,2]
#> [1,] 1.3978842 0.5782369
#> [2,] 0.5782369 1.6636513
sigma0_l # sigma0_l <- t(chol(sigma0))
#>           [,1]     [,2]
#> [1,] 1.1823215 0.000000
#> [2,] 0.4890691 1.193509
alpha
#> [1] 2.87 2.04
beta
#>        [,1]   [,2]
#> [1,]  0.280 -0.048
#> [2,] -0.035  0.260
psi
#>           [,1]      [,2]
#> [1,] 1.3000000 0.5696315
#> [2,] 0.5696315 1.5600000
psi_l # psi_l <- t(chol(psi))
#>           [,1]     [,2]
#> [1,] 1.1401754 0.000000
#> [2,] 0.4995998 1.144727
nu
#> [1] 0 0
lambda
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
theta
#>      [,1] [,2]
#> [1,]  0.2  0.0
#> [2,]  0.0  0.2
theta_l # theta_l <- t(chol(theta))
#>           [,1]      [,2]
#> [1,] 0.4472136 0.0000000
#> [2,] 0.0000000 0.4472136

Visualizing the Dynamics Without Process Noise (n = 5 with Different Initial Condition)

Using the `SimSSMFixed` Function from the `simStateSpace` Package to Simulate Data

library(simStateSpace)
sim <- SimSSMFixed(
  n = n,
  time = time,
  mu0 = mu0,
  sigma0_l = sigma0_l,
  alpha = alpha,
  beta = beta,
  psi_l = psi_l,
  nu = nu,
  lambda = lambda,
  theta_l = theta_l
)
data <- as.data.frame(sim)
head(data)
#>   id time       y1        y2
#> 1  1    0 3.455114 1.9223092
#> 2  1    1 5.696312 0.6452660
#> 3  1    2 4.364664 1.4954592
#> 4  1    3 4.900281 2.6730811
#> 5  1    4 2.325515 0.2139305
#> 6  1    5 3.340405 1.9930331
plot(sim)

Model Fitting

The FitDTVARMx function fits a DT-VAR model assuming fixed parameters.

library(fitDTVARMx)
fit <- FitDTVARMx(
  data = data,
  observed = paste0("y", seq_len(k)),
  id = "id",
  ncores = parallel::detectCores()
)
#> Running DTVAR with 16 parameters
#> 
#> Beginning initial fit attempt
#> Running DTVAR with 16 parameters
#> 
#>  Lowest minimum so far:  6560584.96448053
#> 
#> Solution found

#> 
#>  Solution found!  Final fit=6560585 (started at 112160520)  (1 attempt(s): 1 valid, 0 errors)
#>  Start values from best fit:
#> 0.238845204438278,0.385079342937541,0.192341329589063,-0.319501502562121,2.40734832712481,1.93273981750249,0.765165356271534,0.552642713133481,7.5802858677544e-14,0.730254135746907,1.89415411434006,1,1,7.94885699044023,4.61237397585988,2.14952390785434
summary(fit)
#> Summary of DTVAR 
#>  
#> free parameters:
#>         name         matrix  row  col      Estimate   Std.Error A
#> 1    beta_11   DTVAR_1.beta eta1 eta1  2.388452e-01          NA  
#> 2    beta_21   DTVAR_1.beta eta2 eta1  3.850793e-01          NA  
#> 3    beta_12   DTVAR_1.beta eta1 eta2  1.923413e-01          NA  
#> 4    beta_22   DTVAR_1.beta eta2 eta2 -3.195015e-01          NA  
#> 5    alpha_1  DTVAR_1.gamma eta1    1  2.407348e+00          NA  
#> 6    alpha_2  DTVAR_1.gamma eta2    1  1.932740e+00          NA  
#> 7     psi_11    DTVAR_1.psi eta1 eta1  7.651654e-01          NA !
#> 8     psi_21    DTVAR_1.psi eta1 eta2  5.526427e-01 0.001584757  
#> 9     psi_22    DTVAR_1.psi eta2 eta2  7.580286e-14          NA !
#> 10  theta_11  DTVAR_1.theta   y1   y1  7.302541e-01          NA  
#> 11  theta_22  DTVAR_1.theta   y2   y2  1.894154e+00          NA  
#> 12     mu0_1    DTVAR_1.mu0 eta1  mu0  1.000000e+00          NA !
#> 13     mu0_2    DTVAR_1.mu0 eta2  mu0  1.000000e+00 0.095680038 !
#> 14 sigma0_11 DTVAR_1.sigma0    1    1  7.948857e+00          NA !
#> 15 sigma0_21 DTVAR_1.sigma0    1    2  4.612374e+00 0.272433470 !
#> 16 sigma0_22 DTVAR_1.sigma0    2    2  2.149524e+00 0.304208784 !
#>                  lbound ubound
#> 1                             
#> 2                             
#> 3                             
#> 4                             
#> 5                             
#> 6                             
#> 7  2.2250738585072e-308       
#> 8                             
#> 9                    0!       
#> 10 2.2250738585072e-308       
#> 11 2.2250738585072e-308       
#> 12                   -1     1!
#> 13                   -1     1!
#> 14 2.2250738585072e-308       
#> 15                            
#> 16 2.2250738585072e-308       
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
#>        Model:             16                1999984               6560585
#>    Saturated:             NA                     NA                    NA
#> Independence:             NA                     NA                    NA
#> Number of observations/statistics: 1e+06/2e+06
#> 
#> Information Criteria: 
#>       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
#> AIC:        2560617                6560617                  6560617
#> BIC:      -21070215                6560806                  6560755
#> To get additional fit indices, see help(mxRefModels)
#> timestamp: 2025-08-03 20:34:21 
#> Wall clock time: 3704.049 secs 
#> optimizer:  SLSQP 
#> OpenMx version number: 2.22.7 
#> Need help?  See help(mxSummary)

References

Bringmann, L. F., Vissers, N., Wichers, M., Geschwind, N., Kuppens, P., Peeters, F., Borsboom, D., & Tuerlinckx, F. (2013). A network approach to psychopathology: New insights into clinical longitudinal data. PLoS ONE, 8(4). https://doi.org/10.1371/journal.pone.0060188

Hunter, M. D. (2017). State space modeling in an open source, modular, structural equation modeling environment. Structural Equation Modeling: A Multidisciplinary Journal, 25(2), 307–324. https://doi.org/10.1080/10705511.2017.1369354

Neale, M. C., Hunter, M. D., Pritikin, J. N., Zahery, M., Brick, T. R., Kirkpatrick, R. M., Estabrook, R., Bates, T. C., Maes, H. H., & Boker, S. M. (2015). OpenMx 2.0: Extended structural equation and statistical modeling. Psychometrika, 81(2), 535–549. https://doi.org/10.1007/s11336-014-9435-8

R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/

Ivan Jacob Agaloos Pesigan

2025-08-03