Skip to contents

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 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 <- 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/