## 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.

```{r, include = FALSE, message = FALSE, warning = FALSE, results = 'hide'}
pkg <- c("dynr", "expm", "simStateSpace", "cTMed")
for (i in seq_along(pkg)) {
  if (!require(pkg[i], character.only = TRUE, quietly = TRUE)) {
    install.packages(pkg[i], dependencies = TRUE)
  }
}
```

```{r}
library(dynr)
library(expm)
library(simStateSpace)
```

## 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)

```{r}
#| echo = FALSE
set.seed(42)
n <- 20
time <- 100
delta_t <- 0.10
k <- p <- 3
iden <- diag(k)
null_vec <- rep(x = 0, times = k)
null_mat <- matrix(
  data = 0,
  nrow = p,
  ncol = p
)
mu0 <- null_vec
sigma0 <- matrix(
  data = c(
    1.0, 0.2, 0.2,
    0.2, 1.0, 0.2,
    0.2, 0.2, 1.0
  ),
  nrow = 3
)
sigma0_l <- t(chol(sigma0))
mu <- null_vec
phi <- matrix(
  data = c(
    -0.357, 0.771, -0.450,
    0.0, -0.511, 0.729,
    0, 0, -0.693
  ),
  nrow = p
)
beta_var1 <- expm::expm(
  phi
)
sigma <- matrix(
  data = c(
    0.24455556, 0.02201587, -0.05004762,
    0.02201587, 0.07067800, 0.01539456,
    -0.05004762, 0.01539456, 0.07553061
  ),
  nrow = p
)
sigma_l <- t(chol(sigma))
nu <- null_vec
lambda <- iden
theta <- 0.2 * iden
theta_l <- t(chol(theta))
```

## 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.

```{r}
n
time
delta_t
mu0
sigma0
sigma0_l # sigma0_l <- t(chol(sigma0))
mu
phi
sigma
sigma_l # sigma_l <- t(chol(sigma))
nu
lambda
theta
theta_l # theta_l <- t(chol(theta))
```

```{r}
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
)
plot(sim)
data <- as.data.frame(sim)
colnames(data) <- c("id", "time", "x", "m", "y")
head(data)
```

## Fitting 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

```{r}
dynr_data <- dynr.data(
  dataframe = data,
  id = "id",
  time = "time",
  observed = c("x", "m", "y")
)
```

### Prepare the Initial Condition

```{r}
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

```{r}
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

```{r}
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

```{r}
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.

```{r}
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
```

### Model Estimation

We estimate the parameters using `dynr.cook()`.

```{r}
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)
}
summary(fit)
```

### Extracting 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.

```{r}
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

```{r}
phi
vcov_phi_vec
```

#### Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix

```{r}
sigma
vcov_sigma_vech
```

#### Estimated Drift Matrix and Process Noise Covariance Matrix with Corresponding Sampling Covariance Matrix

```{r}
theta
vcov_theta
```

```{r}
#| include = FALSE
unlink(c("fit-ct-var-dynr.c", "fit-ct-var-dynr.s", "fit-ct-var-dynr.so", "fit-ct-var-dynr.o"))
```