Comparing Copula VAR Models

Overview

The dcvar package provides three copula VAR models:

DC-VAR (dcvar()): Time-varying correlation via a random walk on the Fisher-z scale. Best for smooth, gradual changes.
HMM Copula (dcvar_hmm()): Discrete regime-switching with K hidden states. Best for abrupt changes between distinct coupling regimes.
Constant Copula (dcvar_constant()): Single time-invariant correlation. Serves as a null/baseline model.

Mathematical Background

For normal-margin fits, all three models share a common decomposition of the bivariate likelihood into marginal and copula components. The package also supports exponential, gamma, and skew-normal margins for these three models, but the copula structure below is the common part. All three single-level models additionally support per-variable margins, where each variable can use a different marginal family (see the mixed-margins section below).

VAR(1) Marginals

Each variable follows a VAR(1) process with innovations:

\[y_t = \mu + \Phi \, (y_{t-1} - \mu) + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \text{diag}(\sigma_{\varepsilon}^2))\]

The normal marginal densities provide the per-variable likelihoods in this special case. For non-normal margins, dcvar uses the same VAR recursion and Gaussian copula but swaps in the corresponding marginal likelihood and transform.

Gaussian Copula

The dependence between variables is captured by a Gaussian copula with time-varying correlation \(\rho_t\). The log-likelihood contribution of the copula at time \(t\) is:

\[\log c(u_{1t}, u_{2t}; \rho_t) = -\frac{1}{2}\log(1 - \rho_t^2) - \frac{\rho_t^2 \, (z_{1t}^2 + z_{2t}^2) - 2\rho_t \, z_{1t} z_{2t}}{2(1 - \rho_t^2)}\]

where \(z_{it} = \Phi^{-1}(u_{it})\) and \(u_{it} = \Phi(\varepsilon_{it} / \sigma_{\varepsilon,i})\).

How the Models Differ

DC-VAR: The correlation evolves as a random walk on the Fisher-z scale: \(z_t = z_{t-1} + \omega_t\) with \(\omega_t \sim N(0, \sigma_\omega^2)\), and \(\rho_t = \tanh(z_t)\).
HMM Copula: The correlation is state-dependent: \(\rho_t = \rho_{s_t}\), where \(s_t \in \{1, \ldots, K\}\) follows a Markov chain with transition matrix \(A\). State-specific correlations use an ordered constraint (\(\rho_1 < \rho_2 < \ldots\)) to prevent label switching.
Constant Copula: \(\rho_t = \rho\) for all \(t\). A special case of both the DC-VAR (\(\sigma_\omega = 0\)) and HMM (\(K = 1\)).

Per-variable (mixed) margins

A copula separates the marginal distributions from the dependence structure, so each variable is free to follow its own marginal family. dcvar_constant(), dcvar(), and dcvar_hmm() expose this by accepting a length-2 margins vector. For example, a roughly-symmetric first series paired with a positively-skewed second series can be modelled as normal and exponential margins under a single Gaussian copula:

library(dcvar)

sim_mix <- simulate_dcvar(
  n_time = 200,
  rho_trajectory = rho_constant(200, rho = 0.6),
  margins = c("normal", "exponential"),
  skew_direction = c(1, 1),
  seed = 11
)

fit_mix <- dcvar_constant(
  sim_mix$Y_df,
  vars = c("y1", "y2"),
  margins = c("normal", "exponential"),
  skew_direction = c(1, 1),
  seed = 11
)

# Each dimension is reported under its own family: sigma_eps for the normal
# variable, sigma_exp for the exponential variable.
coef(fit_mix)[c("sigma_eps", "sigma_exp")]

skew_direction is required whenever any dimension uses an exponential or gamma margin, and only those dimensions consult it. A single margins string still applies the same family to both variables, and an all-identical vector such as c("normal", "normal") is treated exactly like the scalar form (it reuses the specialised single-family model), so existing code is unaffected. Mixed margins currently require the Gaussian copula and are available for the three single-level models (dcvar_constant(), dcvar(), dcvar_hmm()).

Time-varying coefficients and scales

In dcvar() the copula correlation always evolves over time. The autoregressive dynamics and the residual scales can be allowed to evolve as well, each as a random walk around a constant baseline. Two flags turn these on:

tv_phi lets the VAR(1) coefficients drift. tv_phi = TRUE varies all four (the autoregressive coefficients phi11, phi22 and the cross-lagged coefficients phi12, phi21); a selector varies only a subset. Use "ar" to let only the autoregressive coefficients drift (for example to capture changing inertia), "cross" for only the cross-lagged coefficients, or pass specific names such as c("phi11", "phi22").
tv_sigma lets the residual scales drift. Normal and skew-normal dimensions scale multiplicatively on the log scale. Exponential and gamma dimensions use a soft-barrier transform so their scale can vary even though the shifted margin has a hard support boundary; tv_sigma_k controls how closely it approximates the exact margin.

With both flags off, dcvar() is the usual dynamic copula model with constant coefficients and scales. Smaller models are nested in larger ones, which suggests fitting a short ladder and comparing the fits.

library(dcvar)

# A series with declining coupling to illustrate the time-varying fits
sim_tv <- simulate_dcvar(
  n_time = 150,
  rho_trajectory = rho_decreasing(150, rho_start = 0.7, rho_end = 0.3),
  seed = 3
)

# Only the autoregressive coefficients evolve over time
fit_ar <- dcvar(
  sim_tv$Y_df,
  vars = c("y1", "y2"),
  tv_phi = "ar",
  seed = 4
)

# Time-varying coefficients and time-varying scales together, with a
# right-skewed second margin (the exponential dimension uses the soft barrier)
fit_full <- dcvar(
  sim_tv$Y_df,
  vars = c("y1", "y2"),
  margins = c("normal", "exponential"),
  skew_direction = c(1, 1),
  tv_phi = TRUE,
  tv_sigma = TRUE,
  seed = 5
)

The coefficient and scale paths come back as posterior summaries, in the same shape as the correlation trajectory, and have their own plotting helpers.

# Posterior summaries of phi11(t), phi12(t), phi21(t), phi22(t)
head(phi_trajectory(fit_ar))

# Per-variable residual-scale paths
head(sigma_trajectory(fit_full))

# Facetted trajectory plots
plot_phi_trajectory(fit_ar)
plot_sigma_trajectory(fit_full)

Fitting All Three Models

library(dcvar)

# Simulate data with a step-function trajectory
sim <- simulate_dcvar(
  n_time = 200,
  rho_trajectory = rho_step(200, rho_before = 0.7, rho_after = 0.3),
  seed = 42
)

# Fit all three models on the same data
fit_dc  <- dcvar(sim$Y_df, vars = c("y1", "y2"), seed = 1)
fit_hmm <- dcvar_hmm(sim$Y_df, vars = c("y1", "y2"), K = 2, seed = 2)
fit_con <- dcvar_constant(sim$Y_df, vars = c("y1", "y2"), seed = 3)

LOO-CV Comparison

dcvar_compare(dcvar = fit_dc, hmm = fit_hmm, constant = fit_con)

Extracting Tidy Summaries

All fit objects support as.data.frame() for exporting tidy parameter summaries. All five fit families also support fitted()/predict() for one-step-ahead values (for SEM and multilevel fits as well). Prediction intervals are restricted to normal margins for single-level, multilevel, and naive-SEM fits; indicator-SEM predictions are available for all supported SEM margins:

# Full parameter summary as a data frame
param_df <- as.data.frame(fit_dc)
head(param_df)

# Predictions with marginal intervals
pred_df <- predict(fit_hmm)
head(pred_df)

HMM-Specific Outputs

The HMM model provides additional outputs related to state inference:

# State posteriors, Viterbi path, transition matrix
states <- hmm_states(fit_hmm)

# State-specific rho values
states$rho_state

# Transition matrix
states$A

# Visualise state posteriors
plot(fit_hmm, type = "states")

# Transition matrix heatmap
plot(fit_hmm, type = "transition")

Clinical Interpretation

interpret_rho_trajectory() provides a model-aware narrative:

interpret_rho_trajectory(fit_dc)
interpret_rho_trajectory(fit_hmm)
interpret_rho_trajectory(fit_con)

When to Use Which Model

Smooth dynamics (therapy effects, developmental trends): Use DC-VAR
Abrupt regime switches (treatment onset, relapse): Use HMM Copula
No time variation expected: Use Constant Copula as baseline
Uncertain: Fit all three and compare with LOO-CV