Fit the Continuous-Time Vector Autoregressive Model Using the dynr Package
Ivan Jacob Agaloos Pesigan
2025-08-23
Source:vignettes/fit-ct-var-dynr.Rmd
fit-ct-var-dynr.Rmd
Data 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.036033
Model 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 <- ub
Fit 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.3518738 0.7442754 -0.4586681 0.01733165
#> -0.4888138 0.7267855 -0.02382127 -0.009818147 -0.6883228 -1.418055 0.09609475
#> -0.2088293 -2.681166 0.2898423 -2.881308 -1.615149 -1.611839 -1.603597
#> 0.006214893 -0.04250065 0.1300648 0.1401327 0.3595736 0.1964612 0.07051581
#> 0.1435534 -0.1096431
#>
#> Transformed fitted parameters: -0.3518738 0.7442754 -0.4586681 0.01733165
#> -0.4888138 0.7267855 -0.02382127 -0.009818147 -0.6883228 0.2421847 0.02327267
#> -0.05057526 0.07071965 0.01498933 0.07237617 0.198861 0.1995204 0.2011716
#> 0.006214893 -0.04250065 0.1300648 1.150426 0.4136629 0.2260142 1.221804
#> 0.2353104 0.9626701
#>
#> Doing end processing
#> Successful trial
#> Total Time: 2.306963
#> Backend Time: 2.306795
summary(fit)
#> Coefficients:
#> Estimate Std. Error t value ci.lower ci.upper Pr(>|t|)
#> phi_1_1 -0.3518738 0.0345599 -10.182 -0.4196100 -0.2841377 <2e-16 ***
#> phi_2_1 0.7442754 0.0213025 34.938 0.7025233 0.7860275 <2e-16 ***
#> phi_3_1 -0.4586681 0.0224194 -20.459 -0.5026093 -0.4147269 <2e-16 ***
#> phi_1_2 0.0173316 0.0304068 0.570 -0.0422645 0.0769278 0.2843
#> phi_2_2 -0.4888138 0.0188377 -25.949 -0.5257350 -0.4518927 <2e-16 ***
#> phi_3_2 0.7267855 0.0198254 36.659 0.6879284 0.7656426 <2e-16 ***
#> phi_1_3 -0.0238213 0.0233767 -1.019 -0.0696387 0.0219962 0.1541
#> phi_2_3 -0.0098181 0.0144220 -0.681 -0.0380848 0.0184485 0.2480
#> phi_3_3 -0.6883228 0.0152967 -44.998 -0.7183038 -0.6583418 <2e-16 ***
#> sigma_1_1 0.2421847 0.0064484 37.557 0.2295461 0.2548232 <2e-16 ***
#> sigma_2_1 0.0232727 0.0025183 9.241 0.0183369 0.0282084 <2e-16 ***
#> sigma_3_1 -0.0505753 0.0026989 -18.739 -0.0558649 -0.0452856 <2e-16 ***
#> sigma_2_2 0.0707197 0.0018899 37.420 0.0670155 0.0744238 <2e-16 ***
#> sigma_3_2 0.0149893 0.0013541 11.070 0.0123354 0.0176433 <2e-16 ***
#> sigma_3_3 0.0723762 0.0020868 34.683 0.0682861 0.0764663 <2e-16 ***
#> theta_1_1 0.1988610 0.0011588 171.608 0.1965898 0.2011323 <2e-16 ***
#> theta_2_2 0.1995204 0.0009996 199.609 0.1975613 0.2014795 <2e-16 ***
#> theta_3_3 0.2011716 0.0010145 198.288 0.1991832 0.2031601 <2e-16 ***
#> mu0_1_1 0.0062149 0.1219710 0.051 -0.2328439 0.2452737 0.4797
#> mu0_2_1 -0.0425006 0.1208284 -0.352 -0.2793200 0.1943187 0.3625
#> mu0_3_1 0.1300648 0.1041761 1.249 -0.0741167 0.3342462 0.1059
#> sigma0_1_1 1.1504264 0.1810442 6.354 0.7955863 1.5052665 <2e-16 ***
#> sigma0_2_1 0.4136629 0.1418008 2.917 0.1357384 0.6915875 0.0018 **
#> sigma0_3_1 0.2260142 0.1187792 1.903 -0.0067889 0.4588172 0.0285 *
#> sigma0_2_2 1.2218038 0.1985814 6.153 0.8325914 1.6110161 <2e-16 ***
#> sigma0_3_2 0.2353104 0.1267133 1.857 -0.0130432 0.4836640 0.0317 *
#> sigma0_3_3 0.9626701 0.1594441 6.038 0.6501655 1.2751747 <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.3518738 0.01733165 -0.023821271
#> [2,] 0.7442754 -0.48881381 -0.009818147
#> [3,] -0.4586681 0.72678550 -0.688322785
vcov_phi_vec
#> phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2
#> phi_1_1 1.194386e-03 9.434520e-05 -1.832722e-04 -1.004483e-03 -8.964831e-05
#> phi_2_1 9.434520e-05 4.537958e-04 -4.543459e-06 -6.913945e-05 -3.836601e-04
#> phi_3_1 -1.832722e-04 -4.543459e-06 5.026282e-04 1.483225e-04 1.026409e-05
#> phi_1_2 -1.004483e-03 -6.913945e-05 1.483225e-04 9.245720e-04 7.282409e-05
#> phi_2_2 -8.964831e-05 -3.836601e-04 1.026409e-05 7.282409e-05 3.548578e-04
#> phi_3_2 1.607647e-04 -2.804753e-06 -4.250983e-04 -1.450161e-04 -2.919274e-06
#> phi_1_3 6.774418e-04 4.293486e-05 -9.128814e-05 -6.468942e-04 -4.622427e-05
#> phi_2_3 7.077153e-05 2.572020e-04 -1.302596e-05 -6.094418e-05 -2.476033e-04
#> phi_3_3 -1.143637e-04 6.496454e-06 2.868236e-04 1.076101e-04 -3.253344e-06
#> phi_3_2 phi_1_3 phi_2_3 phi_3_3
#> phi_1_1 1.607647e-04 6.774418e-04 7.077153e-05 -1.143637e-04
#> phi_2_1 -2.804753e-06 4.293486e-05 2.572020e-04 6.496454e-06
#> phi_3_1 -4.250983e-04 -9.128814e-05 -1.302596e-05 2.868236e-04
#> phi_1_2 -1.450161e-04 -6.468942e-04 -6.094418e-05 1.076101e-04
#> phi_2_2 -2.919274e-06 -4.622427e-05 -2.476033e-04 -3.253344e-06
#> phi_3_2 3.930470e-04 9.435689e-05 8.776401e-06 -2.762247e-04
#> phi_1_3 9.435689e-05 5.464686e-04 4.564710e-05 -8.808095e-05
#> phi_2_3 8.776401e-06 4.564710e-05 2.079953e-04 -1.947278e-06
#> phi_3_3 -2.762247e-04 -8.808095e-05 -1.947278e-06 2.339894e-04
Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix
sigma
#> [,1] [,2] [,3]
#> [1,] 0.24218466 0.02327267 -0.05057526
#> [2,] 0.02327267 0.07071965 0.01498933
#> [3,] -0.05057526 0.01498933 0.07237617
vcov_sigma_vech
#> sigma_1_1 sigma_2_1 sigma_3_1 sigma_2_2 sigma_3_2
#> sigma_1_1 4.158154e-05 6.272750e-09 -5.206212e-06 -3.547992e-07 4.766463e-08
#> sigma_2_1 6.272750e-09 6.341788e-06 4.599255e-08 -3.676482e-07 -9.306099e-07
#> sigma_3_1 -5.206212e-06 4.599255e-08 7.283856e-06 8.207137e-08 -2.042095e-07
#> sigma_2_2 -3.547992e-07 -3.676482e-07 8.207137e-08 3.571711e-06 7.140861e-08
#> sigma_3_2 4.766463e-08 -9.306099e-07 -2.042095e-07 7.140861e-08 1.833491e-06
#> sigma_3_3 4.762538e-07 6.122437e-08 -1.970287e-06 -3.991608e-08 7.177255e-08
#> sigma_3_3
#> sigma_1_1 4.762538e-07
#> sigma_2_1 6.122437e-08
#> sigma_3_1 -1.970287e-06
#> sigma_2_2 -3.991608e-08
#> sigma_3_2 7.177255e-08
#> sigma_3_3 4.354813e-06
Estimated 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.351873815 0.744275429 -0.458668101 0.017331646 -0.488813811 0.726785502
#> phi_1_3 phi_2_3 phi_3_3 sigma_1_1 sigma_2_1 sigma_3_1
#> -0.023821271 -0.009818147 -0.688322785 0.242184655 0.023272674 -0.050575256
#> sigma_2_2 sigma_3_2 sigma_3_3
#> 0.070719655 0.014989335 0.072376167
vcov_theta
#> phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2
#> phi_1_1 1.194386e-03 9.434520e-05 -1.832722e-04 -1.004483e-03 -8.964831e-05
#> phi_2_1 9.434520e-05 4.537958e-04 -4.543459e-06 -6.913945e-05 -3.836601e-04
#> phi_3_1 -1.832722e-04 -4.543459e-06 5.026282e-04 1.483225e-04 1.026409e-05
#> phi_1_2 -1.004483e-03 -6.913945e-05 1.483225e-04 9.245720e-04 7.282409e-05
#> phi_2_2 -8.964831e-05 -3.836601e-04 1.026409e-05 7.282409e-05 3.548578e-04
#> phi_3_2 1.607647e-04 -2.804753e-06 -4.250983e-04 -1.450161e-04 -2.919274e-06
#> phi_1_3 6.774418e-04 4.293486e-05 -9.128814e-05 -6.468942e-04 -4.622427e-05
#> phi_2_3 7.077153e-05 2.572020e-04 -1.302596e-05 -6.094418e-05 -2.476033e-04
#> phi_3_3 -1.143637e-04 6.496454e-06 2.868236e-04 1.076101e-04 -3.253344e-06
#> sigma_1_1 -1.729920e-04 -2.237575e-05 2.852008e-05 1.340106e-04 1.876584e-05
#> sigma_2_1 1.322993e-05 -3.283686e-05 -5.521064e-06 -1.446626e-05 2.438253e-05
#> sigma_3_1 1.541527e-05 4.680862e-06 -3.927416e-05 -9.108262e-06 -3.735018e-06
#> sigma_2_2 1.607424e-06 1.345058e-05 1.017099e-06 -1.221604e-06 -1.378053e-05
#> sigma_3_2 -1.924164e-06 1.728793e-06 7.082798e-06 2.066926e-06 2.530022e-07
#> sigma_3_3 6.711228e-09 -9.215972e-07 4.134147e-06 -8.068151e-07 5.584090e-07
#> phi_3_2 phi_1_3 phi_2_3 phi_3_3 sigma_1_1
#> phi_1_1 1.607647e-04 6.774418e-04 7.077153e-05 -1.143637e-04 -1.729920e-04
#> phi_2_1 -2.804753e-06 4.293486e-05 2.572020e-04 6.496454e-06 -2.237575e-05
#> phi_3_1 -4.250983e-04 -9.128814e-05 -1.302596e-05 2.868236e-04 2.852008e-05
#> phi_1_2 -1.450161e-04 -6.468942e-04 -6.094418e-05 1.076101e-04 1.340106e-04
#> phi_2_2 -2.919274e-06 -4.622427e-05 -2.476033e-04 -3.253344e-06 1.876584e-05
#> phi_3_2 3.930470e-04 9.435689e-05 8.776401e-06 -2.762247e-04 -2.291365e-05
#> phi_1_3 9.435689e-05 5.464686e-04 4.564710e-05 -8.808095e-05 -8.158519e-05
#> phi_2_3 8.776401e-06 4.564710e-05 2.079953e-04 -1.947278e-06 -1.283264e-05
#> phi_3_3 -2.762247e-04 -8.808095e-05 -1.947278e-06 2.339894e-04 1.465593e-05
#> sigma_1_1 -2.291365e-05 -8.158519e-05 -1.283264e-05 1.465593e-05 4.158154e-05
#> sigma_2_1 5.290129e-06 9.403917e-06 -1.409233e-05 -3.513111e-06 6.272750e-09
#> sigma_3_1 2.962004e-05 2.206458e-07 2.046372e-06 -1.675644e-05 -5.206212e-06
#> sigma_2_2 -1.398553e-06 8.112480e-07 8.663912e-06 9.962591e-07 -3.547992e-07
#> sigma_3_2 -7.103547e-06 -1.323634e-06 -2.461897e-06 4.199434e-06 4.766463e-08
#> sigma_3_3 -2.118570e-08 1.882929e-06 1.338675e-07 -5.198197e-06 4.762538e-07
#> sigma_2_1 sigma_3_1 sigma_2_2 sigma_3_2 sigma_3_3
#> phi_1_1 1.322993e-05 1.541527e-05 1.607424e-06 -1.924164e-06 6.711228e-09
#> phi_2_1 -3.283686e-05 4.680862e-06 1.345058e-05 1.728793e-06 -9.215972e-07
#> phi_3_1 -5.521064e-06 -3.927416e-05 1.017099e-06 7.082798e-06 4.134147e-06
#> phi_1_2 -1.446626e-05 -9.108262e-06 -1.221604e-06 2.066926e-06 -8.068151e-07
#> phi_2_2 2.438253e-05 -3.735018e-06 -1.378053e-05 2.530022e-07 5.584090e-07
#> phi_3_2 5.290129e-06 2.962004e-05 -1.398553e-06 -7.103547e-06 -2.118570e-08
#> phi_1_3 9.403917e-06 2.206458e-07 8.112480e-07 -1.323634e-06 1.882929e-06
#> phi_2_3 -1.409233e-05 2.046372e-06 8.663912e-06 -2.461897e-06 1.338675e-07
#> phi_3_3 -3.513111e-06 -1.675644e-05 9.962591e-07 4.199434e-06 -5.198197e-06
#> sigma_1_1 6.272750e-09 -5.206212e-06 -3.547992e-07 4.766463e-08 4.762538e-07
#> sigma_2_1 6.341788e-06 4.599255e-08 -3.676482e-07 -9.306099e-07 6.122437e-08
#> sigma_3_1 4.599255e-08 7.283856e-06 8.207137e-08 -2.042095e-07 -1.970287e-06
#> sigma_2_2 -3.676482e-07 8.207137e-08 3.571711e-06 7.140861e-08 -3.991608e-08
#> sigma_3_2 -9.306099e-07 -2.042095e-07 7.140861e-08 1.833491e-06 7.177255e-08
#> sigma_3_3 6.122437e-08 -1.970287e-06 -3.991608e-08 7.177255e-08 4.354813e-06
References
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
R Core Team. (2024). R: A language and environment for
statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/