| Title: | Regime-Switching Dynamic Correlation Models |
|---|---|
| Description: | Estimation, forecasting, simulation, and portfolio construction for regime-switching models with exogenous variables as in Pelletier (2006) <doi:10.1016/j.jeconom.2005.01.013>. |
| Authors: | David Ardia [aut, cre] (ORCID: <https://orcid.org/0000-0003-2823-782X>), Benjamin Seguin [aut] |
| Maintainer: | David Ardia <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1-2 |
| Built: | 2026-06-03 07:42:54 UTC |
| Source: | https://github.com/ardiad/rsdc |
Daily returns for a green and a brown portfolios constructed following the equal-weighted 10-90 percentile approach.
data(greenbrown)data(greenbrown)
A data frame with 2266 rows and three columns:
Dates ranging from 2014-01-02 to 2024-12-30.
Numeric returns for the green portfolio.
Numeric returns for the brown portfolio.
Originally data in inst/extdata/green-brown-ptf.xlsx.
data("greenbrown") str(greenbrown) head(greenbrown)data("greenbrown") str(greenbrown) head(greenbrown)
The RSDC package provides a comprehensive framework for modeling, estimating, and forecasting correlation structures in multivariate time series under regime-switching dynamics. It supports both fixed transition probabilities and time-varying transition probabilities (TVTP) driven by exogenous variables.
The methodology is particularly suited to empirical asset pricing and portfolio management applications, enabling users to incorporate macroeconomic, financial, or climate-related predictors into the regime dynamics. The package integrates the full workflow — from model estimation to covariance matrix reconstruction and portfolio optimization — in a single, reproducible pipeline.
Model estimation and filtering:
rsdc_hamilton (Hamilton filter),
rsdc_likelihood (likelihood computation),
rsdc_estimate (parameter estimation).
Correlation and covariance forecasting:
rsdc_forecast.
Portfolio construction:
rsdc_minvar (minimum-variance portfolios),
rsdc_maxdiv (maximum-diversification portfolios).
Simulation:
rsdc_simulate (simulate TVTP regime-switching series).
David Ardia and Benjamin Seguin
Engle RF (2002).
“Dynamic conditional correlation: A simple class of multivariate generalized autoregressive conditional heteroskedasticity models.”
Journal of Business & Economic Statistics, 20(3), 339–350.
doi:10.1198/073500102288618487.
Hamilton JD (1989).
“A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle.”
Econometrica, 57(2), 357–384.
doi:10.2307/1912559.
Pelletier D (2006). “Regime switching for dynamic correlations.” Journal of Econometrics, 131(1-2), 445–473. doi:10.1016/j.jeconom.2005.01.004.
Unified front-end that dispatches to one of three estimators:
f_optim() — TVTP specification (method = "tvtp").
f_optim_noX() — fixed transition matrix (method = "noX").
f_optim_const() — constant correlation, single regime (method = "const").
rsdc_estimate( method = c("tvtp", "noX", "const"), residuals, N = 2, X = NULL, out_of_sample = FALSE, control = list() )rsdc_estimate( method = c("tvtp", "noX", "const"), residuals, N = 2, X = NULL, out_of_sample = FALSE, control = list() )
method |
Character. One of |
residuals |
Numeric matrix |
N |
Integer. Number of regimes. Ignored when |
X |
Numeric matrix |
out_of_sample |
Logical. If |
control |
Optional list. Currently forwards |
Method selection: match.arg() validates method.
Inputs: "tvtp" requires non-NULL X; N is ignored for "const".
Split: If out_of_sample = TRUE, the first 70\
transition_matrixEstimated transition matrix ( for "const").
correlationsRegime lower-triangular correlations.
covariancesArray of full correlation matrices.
log_likelihoodMaximized log-likelihood.
betaTVTP coefficients (only for "tvtp").
Mullen K, Ardia D, Gil D, Windover D, Ulrich J (2011).
“DEoptim: An R Package for Global Optimization by Differential Evolution.”
Journal of Statistical Software, 40(6), 1–26.
doi:10.18637/jss.v040.i06.
Hamilton JD (1989).
“A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle.”
Econometrica, 57(2), 357–384.
doi:10.2307/1912559.
Pelletier D (2006). “Regime switching for dynamic correlations.” Journal of Econometrics, 131(1-2), 445–473. doi:10.1016/j.jeconom.2005.01.004.
rsdc_hamilton and rsdc_likelihood.
y <- scale(matrix(rnorm(100 * 3), 100, 3)) rsdc_estimate("const", residuals = y) rsdc_estimate("noX", residuals = y, N = 2) X <- cbind(1, scale(seq_len(nrow(y)))) rsdc_estimate("tvtp", residuals = y, N = 2, X = X)y <- scale(matrix(rnorm(100 * 3), 100, 3)) rsdc_estimate("const", residuals = y) rsdc_estimate("noX", residuals = y, N = 2) X <- cbind(1, scale(seq_len(nrow(y)))) rsdc_estimate("tvtp", residuals = y, N = 2, X = X)
Generates per-period correlation and covariance matrices from a fitted model:
"const" (constant correlation), "noX" (fixed transition matrix), or
"tvtp" (time-varying transition probabilities).
rsdc_forecast( method = c("tvtp", "noX", "const"), N, residuals, X = NULL, final_params, sigma_matrix, value_cols, out_of_sample = FALSE, control = list() )rsdc_forecast( method = c("tvtp", "noX", "const"), N, residuals, X = NULL, final_params, sigma_matrix, value_cols, out_of_sample = FALSE, control = list() )
method |
Character. One of |
N |
Integer. Number of regimes (ignored for |
residuals |
Numeric matrix |
X |
Optional numeric matrix |
final_params |
List with fitted parameters (e.g., from |
sigma_matrix |
Numeric matrix |
value_cols |
Character/integer vector of columns in |
out_of_sample |
Logical. If |
control |
Optional list; supports |
Forecast horizon: If out_of_sample = TRUE, filter on the first threshold
fraction and forecast on the remainder.
Correlation paths:
"const" — empirical correlation of residuals, repeated across time.
"noX"/"tvtp" — smoothed-probability weighted average of regime correlations.
Covariance build: Reconstruct from the pairwise vector (columns ordered by
combn(K, 2)), set , and
.
BIC: Parameter count is
N * ncol(X) + N * K * (K - 1) / 2 for "tvtp",
N * (N - 1) + N * K * (K - 1) / 2 for "noX",
and K * (K - 1) / 2 for "const".
smoothed_probs smoothed probabilities ("noX"/"tvtp" only).
sigma_matrix slice aligned to the forecast horizon.
cov_matricesList of covariance matrices .
predicted_correlations pairwise correlations in combn(K, 2) order.
BICBayesian Information Criterion .
y residual matrix aligned to the forecast horizon.
rsdc_hamilton, rsdc_estimate,
rsdc_minvar, rsdc_maxdiv
set.seed(123) T <- 60; K <- 3; N <- 2 y <- scale(matrix(rnorm(T*K), T, K)) vols <- matrix(0.2 + 0.05*abs(sin(seq_len(T)/7)), T, K) rho <- rbind(c(0.10, 0.05, 0.00), c(0.60, 0.40, 0.30)) Pfix <- matrix(c(0.9, 0.1, 0.2, 0.8), 2, 2, byrow = TRUE) rsdc_forecast("noX", N, y, NULL, list(correlations = rho, transition_matrix = Pfix, log_likelihood = -200), vols, 1:K)set.seed(123) T <- 60; K <- 3; N <- 2 y <- scale(matrix(rnorm(T*K), T, K)) vols <- matrix(0.2 + 0.05*abs(sin(seq_len(T)/7)), T, K) rho <- rbind(c(0.10, 0.05, 0.00), c(0.60, 0.40, 0.30)) Pfix <- matrix(c(0.9, 0.1, 0.2, 0.8), 2, 2, byrow = TRUE) rsdc_forecast("noX", N, y, NULL, list(correlations = rho, transition_matrix = Pfix, log_likelihood = -200), vols, 1:K)
Runs the Hamilton (1989) filter for a multivariate regime-switching correlation model.
Supports either a fixed (time-invariant) transition matrix or time-varying transition
probabilities (TVTP) built from exogenous covariates X via a logistic link.
Returns filtered/smoothed regime probabilities and the log-likelihood.
rsdc_hamilton(y, X = NULL, beta = NULL, rho_matrix, K, N, P = NULL)rsdc_hamilton(y, X = NULL, beta = NULL, rho_matrix, K, N, P = NULL)
y |
Numeric matrix |
X |
Optional numeric matrix |
beta |
Optional numeric matrix |
rho_matrix |
Numeric matrix |
K |
Integer. Number of observed series (columns of |
N |
Integer. Number of regimes. |
P |
Optional |
Correlation rebuild: For regime , a correlation matrix is reconstructed
from rho_matrix[m, ] (lower-triangular fill + symmetrization). Non-PD proposals are penalized.
Transition dynamics:
Fixed P: If X or beta is missing, a constant is used
(user-provided via P; otherwise uniform rows).
TVTP: With X and beta, diagonal entries use
plogis(X[t, ] %*% beta[i, ]). Off-diagonals are equal and sum to .
For , .
Numerical safeguards: A small ridge is added before inversion; if filtering
degenerates at a time step, log_likelihood = -Inf is returned.
A list with:
matrix of filtered probabilities
.
matrix of smoothed probabilities
.
Scalar log-likelihood of the model given y.
TVTP uses a logistic link on the diagonal; off-diagonals are equal by construction.
Hamilton JD (1989). “A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle.” Econometrica, 57(2), 357–384. doi:10.2307/1912559.
rsdc_likelihood and rsdc_estimate.
set.seed(1) T <- 50; K <- 3; N <- 2 y <- scale(matrix(rnorm(T * K), T, K), center = TRUE, scale = TRUE) # Example rho: two regimes with different average correlations rho <- rbind(c(0.10, 0.05, 0.00), c(0.60, 0.40, 0.30)) # lower-tri order for K=3 # Fixed-P filtering Pfix <- matrix(c(0.9, 0.1, 0.2, 0.8), nrow = 2, byrow = TRUE) out_fix <- rsdc_hamilton(y = y, X = NULL, beta = NULL, rho_matrix = rho, K = K, N = N, P = Pfix) str(out_fix$filtered_probs) # TVTP filtering (include an intercept yourself) X <- cbind(1, scale(seq_len(T))) beta <- rbind(c(1.2, 0.0), c(0.8, -0.1)) out_tvtp <- rsdc_hamilton(y = y, X = X, beta = beta, rho_matrix = rho, K = K, N = N) out_tvtp$log_likelihoodset.seed(1) T <- 50; K <- 3; N <- 2 y <- scale(matrix(rnorm(T * K), T, K), center = TRUE, scale = TRUE) # Example rho: two regimes with different average correlations rho <- rbind(c(0.10, 0.05, 0.00), c(0.60, 0.40, 0.30)) # lower-tri order for K=3 # Fixed-P filtering Pfix <- matrix(c(0.9, 0.1, 0.2, 0.8), nrow = 2, byrow = TRUE) out_fix <- rsdc_hamilton(y = y, X = NULL, beta = NULL, rho_matrix = rho, K = K, N = N, P = Pfix) str(out_fix$filtered_probs) # TVTP filtering (include an intercept yourself) X <- cbind(1, scale(seq_len(T))) beta <- rbind(c(1.2, 0.0), c(0.8, -0.1)) out_tvtp <- rsdc_hamilton(y = y, X = X, beta = beta, rho_matrix = rho, K = K, N = N) out_tvtp$log_likelihood
Computes the negative log-likelihood for a multivariate correlation-only regime-switching model, with either a fixed (time-invariant) transition matrix or time-varying transition probabilities (TVTP) driven by exogenous covariates. Likelihood evaluation uses the Hamilton (1989) filter.
rsdc_likelihood(params, y, exog = NULL, K, N)rsdc_likelihood(params, y, exog = NULL, K, N)
params |
Numeric vector of model parameters packed as:
|
y |
Numeric matrix |
exog |
Optional numeric matrix |
K |
Integer. Number of observed series (columns of |
N |
Integer. Number of regimes. |
Transition dynamics:
Fixed P (no exog): params begins with transition
parameters. For , the implementation maps them to
.
TVTP: with exog, diagonal persistence is
; off-diagonals are equal
and sum to .
Correlation build: per regime, the lower-triangular vector is filled into
a symmetric correlation matrix. Non-positive-definite proposals or
are penalized via a large objective value.
Evaluation: delegates to rsdc_hamilton; if the filter returns
log_likelihood = -Inf, a large penalty is returned.
Numeric scalar: the negative log-likelihood to be minimized by an optimizer.
The function is written for use inside optimizers; it performs inexpensive validation and returns large penalties for invalid parameterizations instead of stopping with errors.
rsdc_hamilton (filter),
optim, and DEoptim
# Small toy example (N = 2, K = 3), fixed P (no exog) set.seed(1) T <- 40; K <- 3; N <- 2 y <- scale(matrix(rnorm(T * K), T, K), center = TRUE, scale = TRUE) # Pack parameters: trans (p11, p22), then rho by regime (lower-tri order) p11 <- 0.9; p22 <- 0.8 rho1 <- c(0.10, 0.05, 0.00) # (2,1), (3,1), (3,2) rho2 <- c(0.60, 0.40, 0.30) params <- c(p11, p22, rho1, rho2) nll <- rsdc_likelihood(params, y = y, exog = NULL, K = K, N = N) nll # TVTP example: add X and beta (pack beta row-wise, then rho) X <- cbind(1, scale(seq_len(T))) beta <- rbind(c(1.2, 0.0), c(0.8, -0.1)) params_tvtp <- c(as.vector(t(beta)), rho1, rho2) nll_tvtp <- rsdc_likelihood(params_tvtp, y = y, exog = X, K = K, N = N) nll_tvtp# Small toy example (N = 2, K = 3), fixed P (no exog) set.seed(1) T <- 40; K <- 3; N <- 2 y <- scale(matrix(rnorm(T * K), T, K), center = TRUE, scale = TRUE) # Pack parameters: trans (p11, p22), then rho by regime (lower-tri order) p11 <- 0.9; p22 <- 0.8 rho1 <- c(0.10, 0.05, 0.00) # (2,1), (3,1), (3,2) rho2 <- c(0.60, 0.40, 0.30) params <- c(p11, p22, rho1, rho2) nll <- rsdc_likelihood(params, y = y, exog = NULL, K = K, N = N) nll # TVTP example: add X and beta (pack beta row-wise, then rho) X <- cbind(1, scale(seq_len(T))) beta <- rbind(c(1.2, 0.0), c(0.8, -0.1)) params_tvtp <- c(as.vector(t(beta)), rho1, rho2) nll_tvtp <- rsdc_likelihood(params_tvtp, y = y, exog = X, K = K, N = N) nll_tvtp
Computes rolling maximum-diversification (MaxDiv) portfolio weights from a sequence of per-period covariance matrices implied by forecasted volatilities and correlations. Falls back to equal weights if the nonlinear solver fails.
rsdc_maxdiv(sigma_matrix, value_cols, predicted_corr, y, long_only = TRUE)rsdc_maxdiv(sigma_matrix, value_cols, predicted_corr, y, long_only = TRUE)
sigma_matrix |
Numeric matrix |
value_cols |
Character/integer vector naming columns in |
predicted_corr |
Numeric matrix |
y |
Numeric matrix |
long_only |
Logical. If |
Covariance build: For each , reconstruct from the
pairwise vector; set and
.
Objective (MaxDiv): maximize
subject to and bounds on . Implemented by minimizing
the negative ratio.
Solver: Rsolnp::solnp with equality and
bounds by long_only; on error, weights default to .
weights matrix of weights.
returnsVector of realized portfolio returns sum(y[t,] * weights[t,]).
diversification_ratiosVector of realized diversification ratios.
mean_diversificationAverage diversification ratio.
KNumber of assets.
assetsAsset names.
volatilityStandard deviation of realized portfolio returns.
# Toy example with K = 3 if (requireNamespace("Rsolnp", quietly = TRUE)) { T <- 50; K <- 3 set.seed(42) vols <- matrix(0.2 + 0.05*abs(sin(seq_len(T)/7)), T, K) colnames(vols) <- paste0("A", 1:K) # simple, stationary correlations (order: (2,1), (3,1), (3,2)) pred_corr <- cbind(rep(0.20, T), rep(0.10, T), rep(0.05, T)) rets <- matrix(rnorm(T*K, sd = 0.01), T, K); colnames(rets) <- colnames(vols) mx <- rsdc_maxdiv(sigma_matrix = vols, value_cols = colnames(vols), predicted_corr = pred_corr, y = rets, long_only = TRUE) head(mx$weights) mx$mean_diversification }# Toy example with K = 3 if (requireNamespace("Rsolnp", quietly = TRUE)) { T <- 50; K <- 3 set.seed(42) vols <- matrix(0.2 + 0.05*abs(sin(seq_len(T)/7)), T, K) colnames(vols) <- paste0("A", 1:K) # simple, stationary correlations (order: (2,1), (3,1), (3,2)) pred_corr <- cbind(rep(0.20, T), rep(0.10, T), rep(0.05, T)) rets <- matrix(rnorm(T*K, sd = 0.01), T, K); colnames(rets) <- colnames(vols) mx <- rsdc_maxdiv(sigma_matrix = vols, value_cols = colnames(vols), predicted_corr = pred_corr, y = rets, long_only = TRUE) head(mx$weights) mx$mean_diversification }
Computes rolling minimum-variance (MV) portfolio weights from a sequence of per-period covariance matrices implied by forecasted volatilities and pairwise correlations. Supports long-only or unconstrained MV. If the QP solver fails at a time step, the routine falls back to equal weights.
rsdc_minvar(sigma_matrix, value_cols, predicted_corr, y, long_only = TRUE)rsdc_minvar(sigma_matrix, value_cols, predicted_corr, y, long_only = TRUE)
sigma_matrix |
Numeric matrix |
value_cols |
Character or integer vector giving the columns in |
predicted_corr |
Numeric matrix |
y |
Numeric matrix |
long_only |
Logical. If |
Covariance build: For each , a correlation matrix
is reconstructed ... Let
and .
Optimization: Minimize subject to
and, if long_only,
(solved with quadprog::solve.QP).
Failure handling: If the QP fails at time \(t\), weights default to equal allocation \(w_i = 1/K\).
An object of class "minvar_portfolio":
matrix of MV weights (one row per time).
List of length with the per-period covariance matrices.
Realized standard deviation of portfolio returns (see Note on units).
The input y matrix (coerced to ).
Number of assets.
The realized portfolio return at time \(t\) is computed as sum(y[t, ] * weights[t, ]) / 100.
This assumes y is expressed in \
remove the / 100 in the implementation or convert inputs accordingly.
rsdc_maxdiv (maximum diversification),
solve.QP
# Toy example with K = 3 T <- 50; K <- 3 set.seed(42) vols <- matrix(0.2 + 0.05*abs(sin(seq_len(T)/7)), T, K) colnames(vols) <- paste0("A", 1:K) # simple, stationary correlations pred_corr <- cbind(rep(0.20, T), rep(0.10, T), rep(0.05, T)) # order: (2,1), (3,1), (3,2) rets <- matrix(rnorm(T*K, sd = 0.01), T, K); colnames(rets) <- colnames(vols) mv <- rsdc_minvar(sigma_matrix = vols, value_cols = colnames(vols), predicted_corr= pred_corr, y = rets, long_only = TRUE) head(mv$weights) mv$volatility# Toy example with K = 3 T <- 50; K <- 3 set.seed(42) vols <- matrix(0.2 + 0.05*abs(sin(seq_len(T)/7)), T, K) colnames(vols) <- paste0("A", 1:K) # simple, stationary correlations pred_corr <- cbind(rep(0.20, T), rep(0.10, T), rep(0.05, T)) # order: (2,1), (3,1), (3,2) rets <- matrix(rnorm(T*K, sd = 0.01), T, K); colnames(rets) <- colnames(vols) mv <- rsdc_minvar(sigma_matrix = vols, value_cols = colnames(vols), predicted_corr= pred_corr, y = rets, long_only = TRUE) head(mv$weights) mv$volatility
Simulates a multivariate time series from a regime-switching model with
time-varying transition probabilities (TVTP) driven by covariates X.
Transition probabilities are generated via a multinomial logistic (softmax) link;
observations are drawn from regime-specific Gaussian distributions.
rsdc_simulate(n, X, beta, mu, sigma, N, seed = NULL)rsdc_simulate(n, X, beta, mu, sigma, N, seed = NULL)
n |
Integer. Number of time steps to simulate. |
X |
Numeric matrix |
beta |
Numeric array |
mu |
Numeric matrix |
sigma |
Numeric array |
N |
Integer. Number of regimes. |
seed |
Optional integer. If supplied, sets the RNG seed for reproducibility. |
Initial state and first draw: The initial regime is sampled
uniformly; the first observation is drawn from
.
TVTP via softmax: For , the row of is
computed with log-sum-exp stabilization.
Sampling: Given , draw from the categorical
distribution with probabilities and
.
A list with:
Integer vector of length n; the simulated regime index at each time.
Numeric matrix ; the simulated observations.
Array ; the transition matrix
used at each time step (with undefined by construction; see Details).
Requires mvtnorm for multivariate normal sampling (called as mvtnorm::rmvnorm).
rsdc_hamilton (filter/evaluation),
rsdc_estimate (estimators),
rsdc_forecast (forecasting)
set.seed(123) n <- 200; K <- 3; N <- 2; p <- 2 X <- cbind(1, scale(seq_len(n))) beta <- array(0, dim = c(N, N, p)) beta[1, 1, ] <- c(1.2, 0.0) beta[2, 2, ] <- c(1.0, -0.1) mu <- rbind(c(0, 0, 0), c(0, 0, 0)) rho <- rbind(c(0.10, 0.05, 0.00), c(0.60, 0.40, 0.30)) Sig <- array(0, dim = c(K, K, N)) for (m in 1:N) { R <- diag(K); R[lower.tri(R)] <- rho[m, ]; R[upper.tri(R)] <- t(R)[upper.tri(R)] Sig[, , m] <- R } sim <- rsdc_simulate(n = n, X = X, beta = beta, mu = mu, sigma = Sig, N = N, seed = 99)set.seed(123) n <- 200; K <- 3; N <- 2; p <- 2 X <- cbind(1, scale(seq_len(n))) beta <- array(0, dim = c(N, N, p)) beta[1, 1, ] <- c(1.2, 0.0) beta[2, 2, ] <- c(1.0, -0.1) mu <- rbind(c(0, 0, 0), c(0, 0, 0)) rho <- rbind(c(0.10, 0.05, 0.00), c(0.60, 0.40, 0.30)) Sig <- array(0, dim = c(K, K, N)) for (m in 1:N) { R <- diag(K); R[lower.tri(R)] <- rho[m, ]; R[upper.tri(R)] <- t(R)[upper.tri(R)] Sig[, , m] <- R } sim <- rsdc_simulate(n = n, X = X, beta = beta, mu = mu, sigma = Sig, N = N, seed = 99)