library(dynr)
library(expm)
library(simStateSpace)Fit CT-VAR Using dynr
Installing Required Packages
Before proceeding, ensure that the necessary packages are installed and loaded. We use dynr for model fitting, expm for matrix exponentiation, and simStateSpace for data generation.
Data Generation
We simulate continuous-time data from a three-variable system (\(X\), \(M\), \(Y\)) following an Ornstein–Uhlenbeck (OU) process, which underlies the CT-VAR model.
Each latent variable evolves continuously over time according to a system of stochastic differential equations characterized by:
- Drift matrix \(\boldsymbol{\Phi}\) (system dynamics),
- Process noise covariance \(\boldsymbol{\Sigma}\) (process variability), and
- Measurement noise covariance \(\boldsymbol{\Theta}\) (observation error)
Data Simulation Using simStateSpace
We generate data using the SimSSMOUFixed() function from the simStateSpace package, which simulates trajectories from a continuous-time state-space model with fixed parameters.
n[1] 20time[1] 100delta_t[1] 0.1mu0[1] 0 0 0sigma0 [,1] [,2] [,3]
[1,] 1.0 0.2 0.2
[2,] 0.2 1.0 0.2
[3,] 0.2 0.2 1.0sigma0_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.9660918mu[1] 0 0 0phi [,1] [,2] [,3]
[1,] -0.357 0.000 0.000
[2,] 0.771 -0.511 0.000
[3,] -0.450 0.729 -0.693sigma [,1] [,2] [,3]
[1,] 0.24455556 0.02201587 -0.05004762
[2,] 0.02201587 0.07067800 0.01539456
[3,] -0.05004762 0.01539456 0.07553061sigma_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.243975nu[1] 0 0 0lambda [,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1theta [,1] [,2] [,3]
[1,] 0.2 0.0 0.0
[2,] 0.0 0.2 0.0
[3,] 0.0 0.0 0.2theta_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.4472136sim <- 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
)
plot(sim)
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.036033Fitting the CT-VAR Model Using dynr
We now use the dynr package to estimate the parameters of the CT-VAR model from the simulated data.
Prepare the Data
dynr_data <- dynr.data(
dataframe = data,
id = "id",
time = "time",
observed = c("x", "m", "y")
)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 <- ubModel Estimation
We estimate the parameters using dynr.cook().
fn <- "fit-ct-var-dynr.Rds"
if (file.exists(fn)) {
fit <- readRDS(fn)
} else {
fit <- dynr.cook(
dynr_model,
verbose = FALSE
)
saveRDS(fit, file = fn)
}[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.2150421 0.7992531 -0.5372629 -0.1590345
-0.527781 0.7714035 0.09440367 0.0700308 -0.7375808 -1.531834 0.06888895
-0.2670268 -2.826667 0.1710893 -2.686909 -1.630224 -1.604964 -1.643162
-0.1764606 -0.2433819 0.331227 -0.3120474 -0.1145515 0.140546 0.1462391
0.01082431 0.268019
Transformed fitted parameters: -0.2150421 0.7992531 -0.5372629 -0.1590345
-0.527781 0.7714035 0.09440367 0.0700308 -0.7375808 0.2161389 0.01488958
-0.05771489 0.06023562 0.006154261 0.08523567 0.1958856 0.2008968 0.1933677
-0.1764606 -0.2433819 0.331227 0.7319468 -0.08384563 0.1028722 1.167078
0.0007446744 1.321966
Doing end processing
Successful trial
Total Time: 2.488782
Backend Time: 2.470807 summary(fit)Coefficients:
Estimate Std. Error t value ci.lower ci.upper Pr(>|t|)
phi_1_1 -0.2150421 0.1195297 -1.799 -0.4493160 0.0192318 0.0361 *
phi_2_1 0.7992531 0.0799146 10.001 0.6426233 0.9558830 <2e-16 ***
phi_3_1 -0.5372629 0.0900443 -5.967 -0.7137466 -0.3607793 <2e-16 ***
phi_1_2 -0.1590345 0.0980611 -1.622 -0.3512306 0.0331617 0.0525 .
phi_2_2 -0.5277810 0.0669206 -7.887 -0.6589430 -0.3966190 <2e-16 ***
phi_3_2 0.7714035 0.0748155 10.311 0.6247679 0.9180392 <2e-16 ***
phi_1_3 0.0944037 0.0879781 1.073 -0.0780302 0.2668375 0.1417
phi_2_3 0.0700308 0.0592130 1.183 -0.0460246 0.1860862 0.1185
phi_3_3 -0.7375808 0.0671581 -10.983 -0.8692084 -0.6059533 <2e-16 ***
sigma_1_1 0.2161389 0.0333941 6.472 0.1506876 0.2815902 <2e-16 ***
sigma_2_1 0.0148896 0.0142932 1.042 -0.0131246 0.0429038 0.1488
sigma_3_1 -0.0577149 0.0166399 -3.468 -0.0903285 -0.0251013 0.0003 ***
sigma_2_2 0.0602356 0.0123083 4.894 0.0361118 0.0843594 <2e-16 ***
sigma_3_2 0.0061543 0.0098537 0.625 -0.0131586 0.0254671 0.2662
sigma_3_3 0.0852357 0.0167118 5.100 0.0524810 0.1179903 <2e-16 ***
theta_1_1 0.1958856 0.0078438 24.973 0.1805121 0.2112592 <2e-16 ***
theta_2_2 0.2008968 0.0070651 28.435 0.1870494 0.2147442 <2e-16 ***
theta_3_3 0.1933677 0.0070266 27.519 0.1795958 0.2071396 <2e-16 ***
mu0_1_1 -0.1764606 0.1990029 -0.887 -0.5664990 0.2135778 0.1877
mu0_2_1 -0.2433819 0.2454784 -0.991 -0.7245107 0.2377469 0.1608
mu0_3_1 0.3312270 0.2619670 1.264 -0.1822190 0.8446729 0.1031
sigma0_1_1 0.7319468 0.2506898 2.920 0.2406039 1.2232897 0.0018 **
sigma0_2_1 -0.0838456 0.2213391 -0.379 -0.5176622 0.3499710 0.3524
sigma0_3_1 0.1028722 0.2350566 0.438 -0.3578302 0.5635746 0.3308
sigma0_2_2 1.1670776 0.3852839 3.029 0.4119350 1.9222202 0.0012 **
sigma0_3_2 0.0007447 0.2904693 0.003 -0.5685648 0.5700541 0.4990
sigma0_3_3 1.3219659 0.4366714 3.027 0.4661057 2.1778261 0.0012 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-2 log-likelihood value at convergence = 8608.96
AIC = 8662.96
BIC = 8814.18Extracting Matrices for Use in cTMed
After fitting, we extract the estimated drift matrix (\(\boldsymbol{\Phi}\)) and process noise covariance (\(\boldsymbol{\Sigma}\)), along with their sampling covariance matrices.
coefs <- coef(fit)
vcovs <- vcov(fit)
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
)
colnames(phi) <- rownames(phi) <- c("x", "m", "y")
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 x m y
x -0.2150421 -0.1590345 0.09440367
m 0.7992531 -0.5277810 0.07003080
y -0.5372629 0.7714035 -0.73758084vcov_phi_vec phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2
phi_1_1 0.0142873522 0.0006079799 -0.0030038849 -0.0097797943 -4.790491e-04
phi_2_1 0.0006079799 0.0063863490 -0.0004996187 -0.0001104026 -4.577211e-03
phi_3_1 -0.0030038849 -0.0004996187 0.0081079793 0.0018766011 4.121789e-04
phi_1_2 -0.0097797943 -0.0001104026 0.0018766011 0.0096159753 2.870124e-04
phi_2_2 -0.0004790491 -0.0045772109 0.0004121789 0.0002870124 4.478369e-03
phi_3_2 0.0020607367 0.0002166850 -0.0056919613 -0.0020116829 -3.373929e-04
phi_1_3 0.0066272153 -0.0000384865 -0.0010966857 -0.0065404706 -5.471564e-05
phi_2_3 0.0003941462 0.0031845277 -0.0003412820 -0.0003010944 -3.150796e-03
phi_3_3 -0.0014572097 -0.0001184132 0.0039009020 0.0014237474 1.646252e-04
phi_3_2 phi_1_3 phi_2_3 phi_3_3
phi_1_1 0.0020607367 6.627215e-03 0.0003941462 -0.0014572097
phi_2_1 0.0002166850 -3.848650e-05 0.0031845277 -0.0001184132
phi_3_1 -0.0056919613 -1.096686e-03 -0.0003412820 0.0039009020
phi_1_2 -0.0020116829 -6.540471e-03 -0.0003010944 0.0014237474
phi_2_2 -0.0003373929 -5.471564e-05 -0.0031507962 0.0001646252
phi_3_2 0.0055973556 1.189036e-03 0.0003301834 -0.0038753551
phi_1_3 0.0011890355 7.740141e-03 0.0002868850 -0.0016934695
phi_2_3 0.0003301834 2.868850e-04 0.0035061835 -0.0003437688
phi_3_3 -0.0038753551 -1.693470e-03 -0.0003437688 0.0045102156Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix
sigma [,1] [,2] [,3]
[1,] 0.21613889 0.014889581 -0.057714887
[2,] 0.01488958 0.060235625 0.006154261
[3,] -0.05771489 0.006154261 0.085235668vcov_sigma_vech sigma_1_1 sigma_2_1 sigma_3_1 sigma_2_2 sigma_3_2
sigma_1_1 1.115167e-03 4.004670e-05 -1.981839e-04 1.254314e-05 -1.364776e-05
sigma_2_1 4.004670e-05 2.042965e-04 -1.184950e-05 1.534016e-05 -4.509985e-05
sigma_3_1 -1.981839e-04 -1.184950e-05 2.768867e-04 -4.248350e-06 1.226255e-05
sigma_2_2 1.254314e-05 1.534016e-05 -4.248350e-06 1.514936e-04 -8.509441e-06
sigma_3_2 -1.364776e-05 -4.509985e-05 1.226255e-05 -8.509441e-06 9.709484e-05
sigma_3_3 3.967862e-05 7.242787e-06 -1.124499e-04 5.840479e-07 -1.087946e-05
sigma_3_3
sigma_1_1 3.967862e-05
sigma_2_1 7.242787e-06
sigma_3_1 -1.124499e-04
sigma_2_2 5.840479e-07
sigma_3_2 -1.087946e-05
sigma_3_3 2.792858e-04Estimated 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.215042116 0.799253150 -0.537262944 -0.159034458 -0.527781028 0.771403519
phi_1_3 phi_2_3 phi_3_3 sigma_1_1 sigma_2_1 sigma_3_1
0.094403674 0.070030797 -0.737580843 0.216138890 0.014889581 -0.057714887
sigma_2_2 sigma_3_2 sigma_3_3
0.060235625 0.006154261 0.085235668 vcov_theta phi_1_1 phi_2_1 phi_3_1 phi_1_2 phi_2_2
phi_1_1 1.428735e-02 6.079799e-04 -3.003885e-03 -9.779794e-03 -4.790491e-04
phi_2_1 6.079799e-04 6.386349e-03 -4.996187e-04 -1.104026e-04 -4.577211e-03
phi_3_1 -3.003885e-03 -4.996187e-04 8.107979e-03 1.876601e-03 4.121789e-04
phi_1_2 -9.779794e-03 -1.104026e-04 1.876601e-03 9.615975e-03 2.870124e-04
phi_2_2 -4.790491e-04 -4.577211e-03 4.121789e-04 2.870124e-04 4.478369e-03
phi_3_2 2.060737e-03 2.166850e-04 -5.691961e-03 -2.011683e-03 -3.373929e-04
phi_1_3 6.627215e-03 -3.848650e-05 -1.096686e-03 -6.540471e-03 -5.471564e-05
phi_2_3 3.941462e-04 3.184528e-03 -3.412820e-04 -3.010944e-04 -3.150796e-03
phi_3_3 -1.457210e-03 -1.184132e-04 3.900902e-03 1.423747e-03 1.646252e-04
sigma_1_1 -2.333538e-03 -3.147503e-04 5.590589e-04 1.229505e-03 1.558647e-04
sigma_2_1 1.860527e-04 -4.867247e-04 -7.184730e-05 -2.533846e-04 2.364259e-04
sigma_3_1 2.316315e-04 7.797280e-05 -6.592973e-04 -3.085935e-05 -2.232206e-05
sigma_2_2 -3.241493e-05 1.690887e-04 3.090526e-05 2.853755e-05 -2.312286e-04
sigma_3_2 -5.186949e-06 3.832628e-05 9.180382e-05 2.027968e-05 3.420616e-05
sigma_3_3 -3.757908e-05 -4.591419e-05 8.377694e-05 -7.921819e-07 1.697090e-05
phi_3_2 phi_1_3 phi_2_3 phi_3_3 sigma_1_1
phi_1_1 2.060737e-03 6.627215e-03 3.941462e-04 -1.457210e-03 -2.333538e-03
phi_2_1 2.166850e-04 -3.848650e-05 3.184528e-03 -1.184132e-04 -3.147503e-04
phi_3_1 -5.691961e-03 -1.096686e-03 -3.412820e-04 3.900902e-03 5.590589e-04
phi_1_2 -2.011683e-03 -6.540471e-03 -3.010944e-04 1.423747e-03 1.229505e-03
phi_2_2 -3.373929e-04 -5.471564e-05 -3.150796e-03 1.646252e-04 1.558647e-04
phi_3_2 5.597356e-03 1.189036e-03 3.301834e-04 -3.875355e-03 -2.910012e-04
phi_1_3 1.189036e-03 7.740141e-03 2.868850e-04 -1.693470e-03 -5.863538e-04
phi_2_3 3.301834e-04 2.868850e-04 3.506184e-03 -3.437688e-04 -7.760723e-05
phi_3_3 -3.875355e-03 -1.693470e-03 -3.437688e-04 4.510216e-03 1.448105e-04
sigma_1_1 -2.910012e-04 -5.863538e-04 -7.760723e-05 1.448105e-04 1.115167e-03
sigma_2_1 7.670478e-05 1.414938e-04 -1.119311e-04 -3.768986e-05 4.004670e-05
sigma_3_1 3.251735e-04 -1.831107e-04 -2.155920e-05 -9.609252e-05 -1.981839e-04
sigma_2_2 -3.096602e-05 -1.275493e-05 1.359967e-04 1.458799e-05 1.254314e-05
sigma_3_2 -1.288004e-04 -1.380574e-06 -1.015506e-04 6.378117e-05 -1.364776e-05
sigma_3_3 7.878627e-05 5.620888e-05 1.958059e-05 -2.839006e-04 3.967862e-05
sigma_2_1 sigma_3_1 sigma_2_2 sigma_3_2 sigma_3_3
phi_1_1 1.860527e-04 2.316315e-04 -3.241493e-05 -5.186949e-06 -3.757908e-05
phi_2_1 -4.867247e-04 7.797280e-05 1.690887e-04 3.832628e-05 -4.591419e-05
phi_3_1 -7.184730e-05 -6.592973e-04 3.090526e-05 9.180382e-05 8.377694e-05
phi_1_2 -2.533846e-04 -3.085935e-05 2.853755e-05 2.027968e-05 -7.921819e-07
phi_2_2 2.364259e-04 -2.232206e-05 -2.312286e-04 3.420616e-05 1.697090e-05
phi_3_2 7.670478e-05 3.251735e-04 -3.096602e-05 -1.288004e-04 7.878627e-05
phi_1_3 1.414938e-04 -1.831107e-04 -1.275493e-05 -1.380574e-06 5.620888e-05
phi_2_3 -1.119311e-04 -2.155920e-05 1.359967e-04 -1.015506e-04 1.958059e-05
phi_3_3 -3.768986e-05 -9.609252e-05 1.458799e-05 6.378117e-05 -2.839006e-04
sigma_1_1 4.004670e-05 -1.981839e-04 1.254314e-05 -1.364776e-05 3.967862e-05
sigma_2_1 2.042965e-04 -1.184950e-05 1.534016e-05 -4.509985e-05 7.242787e-06
sigma_3_1 -1.184950e-05 2.768867e-04 -4.248350e-06 1.226255e-05 -1.124499e-04
sigma_2_2 1.534016e-05 -4.248350e-06 1.514936e-04 -8.509441e-06 5.840479e-07
sigma_3_2 -4.509985e-05 1.226255e-05 -8.509441e-06 9.709484e-05 -1.087946e-05
sigma_3_3 7.242787e-06 -1.124499e-04 5.840479e-07 -1.087946e-05 2.792858e-04