Simulation Study Tools

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.