Fit the Continuous-Time Vector Autoregressive Model Using the dynr Package
Ivan Jacob Agaloos Pesigan
2026-01-09
Source:vignettes/fit-ct-var-dynr.Rmd
fit-ct-var-dynr.RmdData Generation Using the SimSSMOUFixed Function from
the simStateSpace Package
n
#> [1] 100
time
#> [1] 1000
delta_t
#> [1] 0.1
mu0
#> [1] 0 0 0
sigma0
#> [,1] [,2] [,3]
#> [1,] 1.0 0.2 0.2
#> [2,] 0.2 1.0 0.2
#> [3,] 0.2 0.2 1.0
sigma0_l # sigma0_l <- t(chol(sigma0))
#> [,1] [,2] [,3]
#> [1,] 1.0 0.0000000 0.0000000
#> [2,] 0.2 0.9797959 0.0000000
#> [3,] 0.2 0.1632993 0.9660918
mu
#> [1] 0 0 0
phi
#> [,1] [,2] [,3]
#> [1,] -0.357 0.000 0.000
#> [2,] 0.771 -0.511 0.000
#> [3,] -0.450 0.729 -0.693
sigma
#> [,1] [,2] [,3]
#> [1,] 0.24455556 0.02201587 -0.05004762
#> [2,] 0.02201587 0.07067800 0.01539456
#> [3,] -0.05004762 0.01539456 0.07553061
sigma_l # sigma_l <- t(chol(sigma))
#> [,1] [,2] [,3]
#> [1,] 0.49452559 0.0000000 0.000000
#> [2,] 0.04451917 0.2620993 0.000000
#> [3,] -0.10120330 0.0759256 0.243975
nu
#> [1] 0 0 0
lambda
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
theta
#> [,1] [,2] [,3]
#> [1,] 0.2 0.0 0.0
#> [2,] 0.0 0.2 0.0
#> [3,] 0.0 0.0 0.2
theta_l # theta_l <- t(chol(theta))
#> [,1] [,2] [,3]
#> [1,] 0.4472136 0.0000000 0.0000000
#> [2,] 0.0000000 0.4472136 0.0000000
#> [3,] 0.0000000 0.0000000 0.4472136
library(simStateSpace)
sim <- SimSSMOUFixed(
n = n,
time = time,
delta_t = delta_t,
mu0 = mu0,
sigma0_l = sigma0_l,
mu = mu,
phi = phi,
sigma_l = sigma_l,
nu = nu,
lambda = lambda,
theta_l = theta_l,
type = 0
)
data <- as.data.frame(sim)
colnames(data) <- c("id", "time", "x", "m", "y")
head(data)
#> id time x m y
#> 1 1 0.0 -0.3504435 0.41877429 2.611996
#> 2 1 0.1 -0.5920330 1.07433208 1.669272
#> 3 1 0.2 -0.7619855 1.21483834 2.369837
#> 4 1 0.3 -1.6964652 0.21209722 2.128531
#> 5 1 0.4 -1.2282686 0.09950326 1.891140
#> 6 1 0.5 0.1433985 0.66784226 2.036033Model Fitting
We use the dynr package to fit the continuous-time
vector autoregressive model.
Prepare the Initial Condition
dynr_initial <- prep.initial(
values.inistate = c(
0, 0, 0
),
params.inistate = c(
"mu0_1_1", "mu0_2_1", "mu0_3_1"
),
values.inicov = matrix(
data = c(
1.0, 0.2, 0.2,
0.2, 1.0, 0.2,
0.2, 0.2, 1.0
),
nrow = 3
),
params.inicov = matrix(
data = c(
"sigma0_1_1", "sigma0_2_1", "sigma0_3_1",
"sigma0_2_1", "sigma0_2_2", "sigma0_3_2",
"sigma0_3_1", "sigma0_3_2", "sigma0_3_3"
),
nrow = 3
)
)Prepare the Measurement Model
dynr_measurement <- prep.measurement(
values.load = diag(3),
params.load = matrix(
data = "fixed",
nrow = 3,
ncol = 3
),
state.names = c("eta_x", "eta_m", "eta_y"),
obs.names = c("x", "m", "y")
)Prepare the Dynamic Model
dynr_dynamics <- prep.formulaDynamics(
formula = list(
eta_x ~ phi_1_1 * eta_x + phi_1_2 * eta_m + phi_1_3 * eta_y,
eta_m ~ phi_2_1 * eta_x + phi_2_2 * eta_m + phi_2_3 * eta_y,
eta_y ~ phi_3_1 * eta_x + phi_3_2 * eta_m + phi_3_3 * eta_y
),
startval = c(
phi_1_1 = -0.2, phi_2_1 = 0.0, phi_3_1 = 0.0,
phi_1_2 = 0.0, phi_2_2 = -0.2, phi_3_2 = 0.0,
phi_1_3 = 0.0, phi_2_3 = 0.0, phi_3_3 = -0.2
),
isContinuousTime = TRUE
)Prepare the Noise Matrices
dynr_noise <- prep.noise(
values.latent = matrix(
data = 0.2 * diag(3),
nrow = 3
),
params.latent = matrix(
data = c(
"sigma_1_1", "sigma_2_1", "sigma_3_1",
"sigma_2_1", "sigma_2_2", "sigma_3_2",
"sigma_3_1", "sigma_3_2", "sigma_3_3"
),
nrow = 3
),
values.observed = 0.2 * diag(3),
params.observed = matrix(
data = c(
"theta_1_1", "fixed", "fixed",
"fixed", "theta_2_2", "fixed",
"fixed", "fixed", "theta_3_3"
),
nrow = 3
)
)Prepare the Model
In this example, we increase the maximum evaluations in the optimization process and constrain the lower and upper bounds of the parameters for the drift matrix.
dynr_model <- dynr.model(
data = dynr_data,
initial = dynr_initial,
measurement = dynr_measurement,
dynamics = dynr_dynamics,
noise = dynr_noise,
outfile = "fit-ct-var-dynr.c"
)
dynr_model@options$maxeval <- 100000
lb <- ub <- rep(NA, times = length(dynr_model$xstart))
names(ub) <- names(lb) <- names(dynr_model$xstart)
lb[
c(
"phi_1_1", "phi_2_1", "phi_3_1",
"phi_1_2", "phi_2_2", "phi_3_2",
"phi_1_3", "phi_2_3", "phi_3_3"
)
] <- -1.5
ub[
c(
"phi_1_1", "phi_2_1", "phi_3_1",
"phi_1_2", "phi_2_2", "phi_3_2",
"phi_1_3", "phi_2_3", "phi_3_3"
)
] <- 1.5
dynr_model$lb <- lb
dynr_model$ub <- ubFit the Model
fit <- dynr.cook(
dynr_model,
verbose = FALSE
)
#> [1] "Get ready!!!!"
#> using C compiler: ‘gcc (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0’
#> Optimization function called.
#> Starting Hessian calculation ...
#> Finished Hessian calculation.
#> Original exit flag: 3
#> Modified exit flag: 3
#> Optimization terminated successfully: ftol_rel or ftol_abs was reached.
#> Original fitted parameters: -0.3518392 0.7442816 -0.4586796 0.01731083
#> -0.4888207 0.7267998 -0.02381434 -0.009810166 -0.6883342 -1.418073 0.09609769
#> -0.2088286 -2.681139 0.2898059 -2.881283 -1.615149 -1.611839 -1.603597
#> 0.006324029 -0.04252988 0.1300433 0.1400118 0.3596041 0.1964667 0.07065217
#> 0.1435497 -0.1097269
#>
#> Transformed fitted parameters: -0.3518392 0.7442816 -0.4586796 0.01731083
#> -0.4888207 0.7267998 -0.02381434 -0.009810166 -0.6883342 0.2421803 0.02327296
#> -0.05057416 0.07072156 0.01498732 0.07237598 0.1988611 0.1995203 0.2011716
#> 0.006324029 -0.04252988 0.1300433 1.150287 0.413648 0.2259932 1.221957
#> 0.2353267 0.962594
#>
#> Doing end processing
#> Successful trial
#> Total Time: 42.69519
#> Backend Time: 42.68898
summary(fit)
#> Coefficients:
#> Estimate Std. Error t value ci.lower ci.upper Pr(>|t|)
#> phi_1_1 -0.351839 0.036416 -9.662 -0.423213 -0.280465 <2e-16 ***
#> phi_2_1 0.744282 0.021777 34.177 0.701599 0.786964 <2e-16 ***
#> phi_3_1 -0.458680 0.023534 -19.490 -0.504806 -0.412554 <2e-16 ***
#> phi_1_2 0.017311 0.031705 0.546 -0.044829 0.079451 0.2925
#> phi_2_2 -0.488821 0.019277 -25.358 -0.526602 -0.451039 <2e-16 ***
#> phi_3_2 0.726800 0.020871 34.824 0.685894 0.767706 <2e-16 ***
#> phi_1_3 -0.023814 0.024025 -0.991 -0.070903 0.023275 0.1608
#> phi_2_3 -0.009810 0.014718 -0.667 -0.038657 0.019036 0.2525
#> phi_3_3 -0.688334 0.016040 -42.913 -0.719773 -0.656896 <2e-16 ***
#> sigma_1_1 0.242180 0.006794 35.646 0.228864 0.255496 <2e-16 ***
#> sigma_2_1 0.023273 0.002545 9.146 0.018285 0.028261 <2e-16 ***
#> sigma_3_1 -0.050574 0.002749 -18.395 -0.055963 -0.045186 <2e-16 ***
#> sigma_2_2 0.070722 0.001907 37.093 0.066985 0.074458 <2e-16 ***
#> sigma_3_2 0.014987 0.001381 10.854 0.012281 0.017694 <2e-16 ***
#> sigma_3_3 0.072376 0.002099 34.475 0.068261 0.076491 <2e-16 ***
#> theta_1_1 0.198861 0.001170 169.909 0.196567 0.201155 <2e-16 ***
#> theta_2_2 0.199520 0.001000 199.500 0.197560 0.201480 <2e-16 ***
#> theta_3_3 0.201172 0.001016 198.052 0.199181 0.203162 <2e-16 ***
#> mu0_1_1 0.006324 0.111110 0.057 -0.211447 0.224095 0.4773
#> mu0_2_1 -0.042530 0.114320 -0.372 -0.266593 0.181533 0.3549
#> mu0_3_1 0.130043 0.102109 1.274 -0.070086 0.330172 0.1014
#> sigma0_1_1 1.150287 0.168811 6.814 0.819425 1.481150 <2e-16 ***
#> sigma0_2_1 0.413648 0.133495 3.099 0.152003 0.675293 0.0010 ***
#> sigma0_3_1 0.225993 0.123478 1.830 -0.016019 0.468006 0.0336 *
#> sigma0_2_2 1.221957 0.182233 6.705 0.864787 1.579128 <2e-16 ***
#> sigma0_3_2 0.235327 0.117629 2.001 0.004779 0.465875 0.0227 *
#> sigma0_3_3 0.962594 0.142152 6.772 0.683981 1.241207 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log-likelihood value at convergence = 429365.49
#> AIC = 429419.49
#> BIC = 429676.34
coefs <- coef(fit)
vcovs <- vcov(fit)Extract Matrices from the Fitted Model to use in cTMed
phi_names <- c(
"phi_1_1", "phi_2_1", "phi_3_1",
"phi_1_2", "phi_2_2", "phi_3_2",
"phi_1_3", "phi_2_3", "phi_3_3"
)
sigma_names <- c(
"sigma_1_1", "sigma_2_1", "sigma_3_1",
"sigma_2_1", "sigma_2_2", "sigma_3_2",
"sigma_3_1", "sigma_3_2", "sigma_3_3"
)
sigma_vech_names <- c(
"sigma_1_1", "sigma_2_1", "sigma_3_1",
"sigma_2_2", "sigma_3_2",
"sigma_3_3"
)
theta_names <- c(
"phi_1_1", "phi_2_1", "phi_3_1",
"phi_1_2", "phi_2_2", "phi_3_2",
"phi_1_3", "phi_2_3", "phi_3_3",
"sigma_1_1", "sigma_2_1", "sigma_3_1",
"sigma_2_2", "sigma_3_2",
"sigma_3_3"
)
phi <- matrix(
data = coefs[phi_names],
nrow = 3,
ncol = 3
)
sigma <- matrix(
data = coefs[sigma_names],
nrow = 3,
ncol = 3
)
theta <- coefs[theta_names]
vcov_phi_vec <- vcovs[phi_names, phi_names]
vcov_sigma_vech <- vcovs[sigma_vech_names, sigma_vech_names]
vcov_theta <- vcovs[theta_names, theta_names]Estimated Drift Matrix with Corresponding Sampling Covariance Matrix
phi
#> [,1] [,2] [,3]
#> [1,] -0.3518392 0.01731083 -0.023814339
#> [2,] 0.7442816 -0.48882067 -0.009810166
#> [3,] -0.4586796 0.72679980 -0.688334177
vcov_phi_vec
#> phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2
#> phi_1_1 1.326121e-03 9.158296e-05 -2.258193e-04 -1.108000e-03 -8.829144e-05
#> phi_2_1 9.158296e-05 4.742430e-04 -3.596064e-06 -6.566903e-05 -4.021299e-04
#> phi_3_1 -2.258193e-04 -3.596064e-06 5.538551e-04 1.845433e-04 1.077730e-05
#> phi_1_2 -1.108000e-03 -6.566903e-05 1.845433e-04 1.005190e-03 7.060966e-05
#> phi_2_2 -8.829144e-05 -4.021299e-04 1.077730e-05 7.060966e-05 3.715859e-04
#> phi_3_2 1.994853e-04 -3.470399e-06 -4.716662e-04 -1.780693e-04 -3.561584e-06
#> phi_1_3 7.414387e-04 3.913166e-05 -1.151977e-04 -6.965449e-04 -4.344659e-05
#> phi_2_3 6.974956e-05 2.704005e-04 -1.452237e-05 -5.931610e-05 -2.595864e-04
#> phi_3_3 -1.424724e-04 7.176532e-06 3.208390e-04 1.318148e-04 -2.963173e-06
#> phi_3_2 phi_1_3 phi_2_3 phi_3_3
#> phi_1_1 1.994853e-04 7.414387e-04 6.974956e-05 -1.424724e-04
#> phi_2_1 -3.470399e-06 3.913166e-05 2.704005e-04 7.176532e-06
#> phi_3_1 -4.716662e-04 -1.151977e-04 -1.452237e-05 3.208390e-04
#> phi_1_2 -1.780693e-04 -6.965449e-04 -5.931610e-05 1.318148e-04
#> phi_2_2 -3.561584e-06 -4.344659e-05 -2.595864e-04 -2.963173e-06
#> phi_3_2 4.355884e-04 1.163975e-04 1.028072e-05 -3.075120e-04
#> phi_1_3 1.163975e-04 5.772234e-04 4.373490e-05 -1.045143e-04
#> phi_2_3 1.028072e-05 4.373490e-05 2.166149e-04 -2.958830e-06
#> phi_3_3 -3.075120e-04 -1.045143e-04 -2.958830e-06 2.572924e-04Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix
sigma
#> [,1] [,2] [,3]
#> [1,] 0.24218026 0.02327296 -0.05057416
#> [2,] 0.02327296 0.07072156 0.01498732
#> [3,] -0.05057416 0.01498732 0.07237598
vcov_sigma_vech
#> sigma_1_1 sigma_2_1 sigma_3_1 sigma_2_2 sigma_3_2
#> sigma_1_1 4.616001e-05 -3.795706e-07 -5.585381e-06 -1.692509e-07 3.054065e-07
#> sigma_2_1 -3.795706e-07 6.475588e-06 1.002384e-07 -4.387889e-07 -9.708718e-07
#> sigma_3_1 -5.585381e-06 1.002384e-07 7.558692e-06 9.058849e-08 -3.163361e-07
#> sigma_2_2 -1.692509e-07 -4.387889e-07 9.058849e-08 3.635137e-06 3.523268e-08
#> sigma_3_2 3.054065e-07 -9.708718e-07 -3.163361e-07 3.523268e-08 1.906780e-06
#> sigma_3_3 6.193436e-07 3.360567e-08 -1.947016e-06 -4.509608e-08 4.773520e-08
#> sigma_3_3
#> sigma_1_1 6.193436e-07
#> sigma_2_1 3.360567e-08
#> sigma_3_1 -1.947016e-06
#> sigma_2_2 -4.509608e-08
#> sigma_3_2 4.773520e-08
#> sigma_3_3 4.407331e-06Estimated Drift Matrix and Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix
theta
#> phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2 phi_3_2
#> -0.351839184 0.744281649 -0.458679613 0.017310833 -0.488820674 0.726799800
#> phi_1_3 phi_2_3 phi_3_3 sigma_1_1 sigma_2_1 sigma_3_1
#> -0.023814339 -0.009810166 -0.688334177 0.242180259 0.023272963 -0.050574160
#> sigma_2_2 sigma_3_2 sigma_3_3
#> 0.070721559 0.014987323 0.072375983
vcov_theta
#> phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2
#> phi_1_1 1.326121e-03 9.158296e-05 -2.258193e-04 -1.108000e-03 -8.829144e-05
#> phi_2_1 9.158296e-05 4.742430e-04 -3.596064e-06 -6.566903e-05 -4.021299e-04
#> phi_3_1 -2.258193e-04 -3.596064e-06 5.538551e-04 1.845433e-04 1.077730e-05
#> phi_1_2 -1.108000e-03 -6.566903e-05 1.845433e-04 1.005190e-03 7.060966e-05
#> phi_2_2 -8.829144e-05 -4.021299e-04 1.077730e-05 7.060966e-05 3.715859e-04
#> phi_3_2 1.994853e-04 -3.470399e-06 -4.716662e-04 -1.780693e-04 -3.561584e-06
#> phi_1_3 7.414387e-04 3.913166e-05 -1.151977e-04 -6.965449e-04 -4.344659e-05
#> phi_2_3 6.974956e-05 2.704005e-04 -1.452237e-05 -5.931610e-05 -2.595864e-04
#> phi_3_3 -1.424724e-04 7.176532e-06 3.208390e-04 1.318148e-04 -2.963173e-06
#> sigma_1_1 -1.982932e-04 -2.184773e-05 3.456919e-05 1.543789e-04 1.849510e-05
#> sigma_2_1 1.502203e-05 -3.437997e-05 -6.620570e-06 -1.591184e-05 2.574007e-05
#> sigma_3_1 1.794495e-05 4.674418e-06 -4.279970e-05 -1.117672e-05 -3.814434e-06
#> sigma_2_2 4.361454e-07 1.424515e-05 1.010140e-06 -1.579752e-07 -1.451605e-05
#> sigma_3_2 -3.741407e-06 1.607279e-06 8.491226e-06 3.585007e-06 4.201294e-07
#> sigma_3_3 -1.043117e-06 -8.607286e-07 3.442663e-06 7.089937e-08 4.858637e-07
#> phi_3_2 phi_1_3 phi_2_3 phi_3_3 sigma_1_1
#> phi_1_1 1.994853e-04 7.414387e-04 6.974956e-05 -1.424724e-04 -1.982932e-04
#> phi_2_1 -3.470399e-06 3.913166e-05 2.704005e-04 7.176532e-06 -2.184773e-05
#> phi_3_1 -4.716662e-04 -1.151977e-04 -1.452237e-05 3.208390e-04 3.456919e-05
#> phi_1_2 -1.780693e-04 -6.965449e-04 -5.931610e-05 1.318148e-04 1.543789e-04
#> phi_2_2 -3.561584e-06 -4.344659e-05 -2.595864e-04 -2.963173e-06 1.849510e-05
#> phi_3_2 4.355884e-04 1.163975e-04 1.028072e-05 -3.075120e-04 -2.844219e-05
#> phi_1_3 1.163975e-04 5.772234e-04 4.373490e-05 -1.045143e-04 -9.445101e-05
#> phi_2_3 1.028072e-05 4.373490e-05 2.166149e-04 -2.958830e-06 -1.266847e-05
#> phi_3_3 -3.075120e-04 -1.045143e-04 -2.958830e-06 2.572924e-04 1.859927e-05
#> sigma_1_1 -2.844219e-05 -9.445101e-05 -1.266847e-05 1.859927e-05 4.616001e-05
#> sigma_2_1 6.267263e-06 1.038566e-05 -1.506230e-05 -4.238014e-06 -3.795706e-07
#> sigma_3_1 3.275583e-05 1.446600e-06 2.154429e-06 -1.897455e-05 -5.585381e-06
#> sigma_2_2 -1.392458e-06 1.948765e-08 9.194261e-06 1.027499e-06 -1.692509e-07
#> sigma_3_2 -8.374628e-06 -2.326402e-06 -2.643065e-06 5.086543e-06 3.054065e-07
#> sigma_3_3 6.621115e-07 1.340588e-06 1.970612e-07 -5.788873e-06 6.193436e-07
#> sigma_2_1 sigma_3_1 sigma_2_2 sigma_3_2 sigma_3_3
#> phi_1_1 1.502203e-05 1.794495e-05 4.361454e-07 -3.741407e-06 -1.043117e-06
#> phi_2_1 -3.437997e-05 4.674418e-06 1.424515e-05 1.607279e-06 -8.607286e-07
#> phi_3_1 -6.620570e-06 -4.279970e-05 1.010140e-06 8.491226e-06 3.442663e-06
#> phi_1_2 -1.591184e-05 -1.117672e-05 -1.579752e-07 3.585007e-06 7.089937e-08
#> phi_2_2 2.574007e-05 -3.814434e-06 -1.451605e-05 4.201294e-07 4.858637e-07
#> phi_3_2 6.267263e-06 3.275583e-05 -1.392458e-06 -8.374628e-06 6.621115e-07
#> phi_1_3 1.038566e-05 1.446600e-06 1.948765e-08 -2.326402e-06 1.340588e-06
#> phi_2_3 -1.506230e-05 2.154429e-06 9.194261e-06 -2.643065e-06 1.970612e-07
#> phi_3_3 -4.238014e-06 -1.897455e-05 1.027499e-06 5.086543e-06 -5.788873e-06
#> sigma_1_1 -3.795706e-07 -5.585381e-06 -1.692509e-07 3.054065e-07 6.193436e-07
#> sigma_2_1 6.475588e-06 1.002384e-07 -4.387889e-07 -9.708718e-07 3.360567e-08
#> sigma_3_1 1.002384e-07 7.558692e-06 9.058849e-08 -3.163361e-07 -1.947016e-06
#> sigma_2_2 -4.387889e-07 9.058849e-08 3.635137e-06 3.523268e-08 -4.509608e-08
#> sigma_3_2 -9.708718e-07 -3.163361e-07 3.523268e-08 1.906780e-06 4.773520e-08
#> sigma_3_3 3.360567e-08 -1.947016e-06 -4.509608e-08 4.773520e-08 4.407331e-06References
Ou, L., Hunter, M. D., & Chow, S.-M. (2019). What’s for dynr: A package for linear and nonlinear dynamic
modeling in R. The R Journal, 11(1), 91.
https://doi.org/10.32614/rj-2019-012
Pesigan, I. J. A., Russell, M. A., & Chow, S.-M. (2025). Inferences
and effect sizes for direct, indirect, and total effects in
continuous-time mediation models. Psychological Methods. https://doi.org/10.1037/met0000779
R Core Team. (2024). R: A language and environment for
statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/