Simulation Study Tools

Trajectory Generators

dcvar provides several trajectory generators for simulation studies. All return a numeric vector of length T-1 (rho values for time steps 2 through n_time):

library(dcvar)

# Smooth trajectories
rho_constant(100, rho = 0.5)
#>  [1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
#> [20] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
#> [39] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
#> [58] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
#> [77] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
#> [96] 0.5 0.5 0.5 0.5
rho_decreasing(100, rho_start = 0.7, rho_end = 0.3)
#>  [1] 0.6682246 0.6667309 0.6651737 0.6635508 0.6618602 0.6600998 0.6582675
#>  [8] 0.6563613 0.6543790 0.6523188 0.6501787 0.6479566 0.6456508 0.6432596
#> [15] 0.6407811 0.6382139 0.6355564 0.6328074 0.6299655 0.6270298 0.6239994
#> [22] 0.6208736 0.6176519 0.6143340 0.6109199 0.6074099 0.6038044 0.6001040
#> [29] 0.5963100 0.5924234 0.5884461 0.5843798 0.5802269 0.5759898 0.5716715
#> [36] 0.5672751 0.5628042 0.5582625 0.5536542 0.5489837 0.5442557 0.5394751
#> [43] 0.5346470 0.5297770 0.5248706 0.5199336 0.5149719 0.5099917 0.5049990
#> [50] 0.5000000 0.4950010 0.4900083 0.4850281 0.4800664 0.4751294 0.4702230
#> [57] 0.4653530 0.4605249 0.4557443 0.4510163 0.4463458 0.4417375 0.4371958
#> [64] 0.4327249 0.4283285 0.4240102 0.4197731 0.4156202 0.4115539 0.4075766
#> [71] 0.4036900 0.3998960 0.3961956 0.3925901 0.3890801 0.3856660 0.3823481
#> [78] 0.3791264 0.3760006 0.3729702 0.3700345 0.3671926 0.3644436 0.3617861
#> [85] 0.3592189 0.3567404 0.3543492 0.3520434 0.3498213 0.3476812 0.3456210
#> [92] 0.3436387 0.3417325 0.3399002 0.3381398 0.3364492 0.3348263 0.3332691
#> [99] 0.3317754
rho_increasing(100, rho_start = 0.3, rho_end = 0.7)
#>  [1] 0.3317754 0.3332691 0.3348263 0.3364492 0.3381398 0.3399002 0.3417325
#>  [8] 0.3436387 0.3456210 0.3476812 0.3498213 0.3520434 0.3543492 0.3567404
#> [15] 0.3592189 0.3617861 0.3644436 0.3671926 0.3700345 0.3729702 0.3760006
#> [22] 0.3791264 0.3823481 0.3856660 0.3890801 0.3925901 0.3961956 0.3998960
#> [29] 0.4036900 0.4075766 0.4115539 0.4156202 0.4197731 0.4240102 0.4283285
#> [36] 0.4327249 0.4371958 0.4417375 0.4463458 0.4510163 0.4557443 0.4605249
#> [43] 0.4653530 0.4702230 0.4751294 0.4800664 0.4850281 0.4900083 0.4950010
#> [50] 0.5000000 0.5049990 0.5099917 0.5149719 0.5199336 0.5248706 0.5297770
#> [57] 0.5346470 0.5394751 0.5442557 0.5489837 0.5536542 0.5582625 0.5628042
#> [64] 0.5672751 0.5716715 0.5759898 0.5802269 0.5843798 0.5884461 0.5924234
#> [71] 0.5963100 0.6001040 0.6038044 0.6074099 0.6109199 0.6143340 0.6176519
#> [78] 0.6208736 0.6239994 0.6270298 0.6299655 0.6328074 0.6355564 0.6382139
#> [85] 0.6407811 0.6432596 0.6456508 0.6479566 0.6501787 0.6523188 0.6543790
#> [92] 0.6563613 0.6582675 0.6600998 0.6618602 0.6635508 0.6651737 0.6667309
#> [99] 0.6682246
rho_random_walk(100, sigma_omega = 0.05, seed = 42)
#>  [1] 0.5142921 0.4932249 0.5068409 0.5299766 0.5443565 0.5406119 0.5919066
#>  [8] 0.5888233 0.6508570 0.6490459 0.6852217 0.7412392 0.7083085 0.7012937
#> [15] 0.6978902 0.7138414 0.7068001 0.6339170 0.5552325 0.5992140 0.5892968
#> [22] 0.5281115 0.5218850 0.5646684 0.6257621 0.6124901 0.6043889 0.5454457
#> [29] 0.5614033 0.5390959 0.5550515 0.5789585 0.6123358 0.5929500 0.6090761
#> [36] 0.5522416 0.5243916 0.4928611 0.3961906 0.3977121 0.4063472 0.3911655
#> [43] 0.4227892 0.3924969 0.3331205 0.3522192 0.3161875 0.3796156 0.3610025
#> [50] 0.3891692 0.4027414 0.3693972 0.4353762 0.4610597 0.4645863 0.4753594
#> [57] 0.5012216 0.5045772 0.3848780 0.3969452 0.3813648 0.3892513 0.4136515
#> [64] 0.4699416 0.4411282 0.4920471 0.5046685 0.5423476 0.5740295 0.5976966
#> [71] 0.5631266 0.5600394 0.5810633 0.5486098 0.5293552 0.5499424 0.5761682
#> [78] 0.5914518 0.5619012 0.5231146 0.5758648 0.5844202 0.5873244 0.5833507
#> [85] 0.5425839 0.5638161 0.5563651 0.5500238 0.5817360 0.6082704 0.6502762
#> [92] 0.6363219 0.6552740 0.6931749 0.6631985 0.6383964 0.6036535 0.5552328
#> [99] 0.5579929

# Step-function trajectories
rho_step(100, rho_before = 0.7, rho_after = 0.3, breakpoint = 0.5)
#>  [1] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [20] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [39] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#> [58] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#> [77] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#> [96] 0.3 0.3 0.3 0.3
rho_double_step(100, rho_levels = c(0.7, 0.3, 0.7))
#>  [1] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [20] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.3 0.3 0.3 0.3 0.3
#> [39] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#> [58] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [77] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [96] 0.7 0.7 0.7 0.7

# Named scenarios
rho_scenario("double_relapse", n_time = 150)
#>   [1] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#>  [19] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#>  [37] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.3 0.3 0.3 0.3
#>  [55] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#>  [73] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#>  [91] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [109] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [127] 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
#> [145] 0.7 0.7 0.7 0.7 0.7

Visualising Trajectories

plot_trajectories() shows all built-in scenarios side by side:

# Compare all built-in scenarios
plot_trajectories(100)


# Subset of scenarios
plot_trajectories(100, scenarios = c("decreasing", "single_middle", "double_relapse"))

Simulating Data

simulate_dcvar() generates bivariate VAR(1) data with a given rho trajectory and known ground-truth parameters:

sim <- simulate_dcvar(
  n_time = 150,
  rho_trajectory = rho_decreasing(150),
  mu = c(0, 0),
  Phi = matrix(c(0.3, 0.1, 0.1, 0.3), 2, 2),
  sigma_eps = c(1, 1),
  seed = 42
)

Breakpoint Data Simulator

simulate_breakpoint_data() provides a simplified interface for generating data with single or double breakpoints:

# Single breakpoint
sim_bp <- simulate_breakpoint_data(n_time = 100, type = "single", seed = 42)
head(sim_bp$Y_df)
#>   time        y1        y2
#> 1    1 0.0000000 0.0000000
#> 2    2 1.3709584 0.5563958
#> 3    3 0.8300555 1.0101588
#> 4    4 0.7543009 0.5932529
#> 5    5 1.7971375 1.2438713
#> 6    6 2.6819521 1.9209849
sim_bp$true_params$breakpoint
#> NULL

# Double breakpoint
sim_bp2 <- simulate_breakpoint_data(n_time = 100, type = "double", seed = 42)

Data Preparation

Before fitting, data must be prepared with prepare_dcvar_data(), prepare_hmm_data(), or prepare_constant_data(). These functions handle standardisation, prior hyperparameters, and validation:

df <- data.frame(time = 1:100, y1 = rnorm(100), y2 = rnorm(100))

# DC-VAR data (includes sigma_omega prior)
stan_data <- prepare_dcvar_data(df, vars = c("y1", "y2"))
str(stan_data[c("n_time", "D", "sigma_mu_prior", "sigma_omega_prior")])
#> List of 4
#>  $ n_time           : int 100
#>  $ D                : int 2
#>  $ sigma_mu_prior   : num 2
#>  $ sigma_omega_prior: num 0.1

# HMM data (includes K, kappa, alpha_off)
hmm_data <- prepare_hmm_data(df, vars = c("y1", "y2"), K = 3)
str(hmm_data[c("K", "kappa", "alpha_off")])
#> List of 3
#>  $ K        : int 3
#>  $ kappa    : num 10
#>  $ alpha_off: num 1

# Constant data (simplest)
const_data <- prepare_constant_data(df, vars = c("y1", "y2"))

Parameter Recovery Metrics

After fitting, evaluate recovery with:

# Extract estimated rho
rho_df <- rho_trajectory(fit)

# Compute recovery metrics
metrics <- compute_rho_metrics(
  rho_true = sim$true_params$rho,
  rho_est = rho_df$mean,
  rho_lower = rho_df$q2.5,
  rho_upper = rho_df$q97.5
)

metrics$bias
metrics$coverage
metrics$correlation

Running a Mini Simulation Study

n_reps <- 5
results <- lapply(1:n_reps, function(i) {
  sim <- simulate_dcvar(
    n_time = 100,
    rho_trajectory = rho_decreasing(100),
    seed = i * 1000
  )

  fit <- dcvar(sim$Y_df, vars = c("y1", "y2"),
               chains = 2, iter_warmup = 500, iter_sampling = 1000,
               seed = i, refresh = 0)

  rho_df <- rho_trajectory(fit)

  list(rho = compute_rho_metrics(
    rho_true = sim$true_params$rho,
    rho_est = rho_df$mean,
    rho_lower = rho_df$q2.5,
    rho_upper = rho_df$q97.5
  ))
})

# Aggregate across replications
agg <- aggregate_metrics(results)
print(agg$rho)