## ----setup, include = FALSE---------------------------------------------------
has_stan <- requireNamespace("rstan", quietly = TRUE)
eval_fits <- identical(Sys.getenv("DCVAR_EVAL_VIGNETTES"), "true") && has_stan
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = eval_fits
)

## ----trajectories, eval = TRUE------------------------------------------------
library(dcvar)

# Smooth trajectories
rho_constant(100, rho = 0.5)
rho_decreasing(100, rho_start = 0.7, rho_end = 0.3)
rho_increasing(100, rho_start = 0.3, rho_end = 0.7)
rho_random_walk(100, sigma_omega = 0.05, seed = 42)

# Step-function trajectories
rho_step(100, rho_before = 0.7, rho_after = 0.3, breakpoint = 0.5)
rho_double_step(100, rho_levels = c(0.7, 0.3, 0.7))

# Named scenarios
rho_scenario("double_relapse", n_time = 150)

## ----plot-trajectories, eval = TRUE-------------------------------------------
# Compare all built-in scenarios
plot_trajectories(100)

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

## ----simulate, eval = TRUE----------------------------------------------------
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, eval = TRUE--------------------------------------------------
# Single breakpoint
sim_bp <- simulate_breakpoint_data(n_time = 100, type = "single", seed = 42)
head(sim_bp$Y_df)
sim_bp$true_params$breakpoint

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

## ----prepare, eval = TRUE-----------------------------------------------------
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")])

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

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

## ----metrics-fit, include = FALSE---------------------------------------------
# fit <- dcvar(
#   sim$Y_df,
#   vars = c("y1", "y2"),
#   chains = 2,
#   iter_warmup = 250,
#   iter_sampling = 250,
#   seed = 99,
#   refresh = 0
# )

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

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