Hypothetical Stress Scenario Engine

Overview

The Hypothetical Stress Scenario Engine is a production-grade, statistically coherent, multi-factor stress-propagation engine built for RBC Capital Markets, summer 2026. The user pins a handful of factors to chosen shocks (e.g. "equities −8%, oil −15%, 10Y rates +25bp"); the engine returns statistically coherent shocks for every remaining factor, calibrated on history and conditioned through the joint covariance and copula.

Two universes

• Universe A — 41 weakly correlated cross-asset factors (equity, crypto, energy, metals, ags, FX, rates). The Gaussian / copula core handles this. This MVP delivers Universe A end-to-end.
• Universe B — highly collinear single-market factors (a rates curve, a gas-hub complex). Direct conditioning is numerically unstable, so we condition in PCA space (demonstrated on the rates complex).

Key Features

• Bloomberg native: a PX_LAST panel spanning 1996–2026, with a Bloomberg-driven segment / tier map (0 unknown segments).
• Modern-regime universe: 41 Tier A+B factors over 2018-06 → 2026-06 (2,082 rows) — deliberately spans the COVID-2020 crash and includes crypto and the SOFR transition, which a long-history Tier-A-only window would exclude.
• Eight covariance estimators: sample, EWMA, EWMA+corr shrinkage, Ledoit–Wolf, cluster-structured block shrinkage, PCA routing, and PCA / Barra factor models — selected by a robust backtest.
• Three engine cores: conditional Gaussian, Gaussian copula, and Student-t copula (tail dependence) with empirical marginals.
• Contract-roll stitching: M1→M2→M3 continuous series with a grid-searched weighted roll, removing the calendar-basis jump.
• Out-of-sample validation: reconstruction of the 20 worst historical stress days, scored against a naive "leave-the-rest-flat" baseline.
• Interactive dashboard: search / pin factors, run Gaussian / copula, and view per-factor quantile bars via ipywidgets.

Motivation

Traditional stress scenarios at financial institutions often rely on:

• Ad-hoc scaling: "if equities drop 20%, then bonds move +500bp" (rule-of-thumb, no covariance).
• Isometric shocks: "all risk factors move the same percentage" (ignores correlations).
• Expert override: "we know what should happen" (no statistical grounding, hard to audit).

These approaches miss spillover effects and produce scenarios that violate historical covariance relationships.

Conditional Gaussian solution

By conditioning on a subset of "pinned" factors, we estimate the joint distribution of all remaining factors using the multivariate Gaussian framework:

μ_B|A = μ_B + Σ_BA Σ_AA^-1 (shock_A − μ_A)

This respects covariance, is repeatable, auditable, and supports validation against historical stress episodes. Because the mapping is linear in the pin, it also extrapolates: pinning a never-seen shock still yields a coherent scenario.

System Architecture

The architecture separates concerns so each piece can be productionised independently: a universe-specific Configuration layer, a universe-agnostic statistical Engine, and an out-of-sample Validation layer.

Layer	Responsibility	Status
1. Configuration	Segment / tier mapping, transforms, pinned-factor choices	Complete
2. Data pipeline	PX_LAST panel, holiday-aware bounded fill, contract-roll stitching	Complete
3. Statistics	Returns, clustering, eight covariance estimators + robust selection	Complete
4. Conditional engine	Schur-complement conditioning + antithetic Monte-Carlo	Complete
5. Copula engine	Gaussian & Student-t copula with empirical marginals	Complete
6. Universe B	PCA-space conditioning for collinear curves (rates / gas hub)	Demonstrated (rates complex)
7. Validation	20-worst-day out-of-sample reconstruction vs naive baseline	Complete
8. Application	ipywidgets dashboard; Streamlit deployment	Dashboard complete · Streamlit Phase 2

The methodology progresses as a strict generalisation at each step: Gaussian MVP → Gaussian copula → Student-t copula (tail dependence) → PCA-space conditioning (collinear blocks), with DCC-GARCH / vine copulas / entropy pooling on the roadmap.

Risk Factor Universe

The engine ingests a cross-asset Bloomberg PX_LAST panel and classifies every instrument into one of seven segments via a Bloomberg-driven segment map (0 unknowns after the remap). Segments are ordered risk-on → diversifier → safe-haven.

• Equity — index futures (S&P 500, Nasdaq, etc.)
• Crypto — BTC / ETH factors on their native 7-day calendar
• Energy — WTI crude (CL1), natural gas (NG1), TTF gas, products
• Metals — gold, silver, copper, and the LME complex
• Agriculture — wheat, corn
• FX — EUR/USD, JPY/USD, CAD/USD futures
• Rates — Treasury / Bund / Schatz curve and short-rate futures

Universe Selection & the N-vs-T Trade-off

Stress lives in the modern regime — crypto, the SOFR transition, the 2022 hiking cycle, the 2025 Strait-of-Hormuz shock — so the engine is built on Tier A+B factors rather than only the long-history mature ones. Adding the young factors (crypto ~2017, SOFR ~2018) forces a shorter complete-case window: more factors N ⇒ fewer usable rows T.

We resolve this explicitly by scanning candidate start dates and choosing the one that maximises factor count subject to (i) ≥ 6 years (1,500 rows) of complete data and (ii) a start on/before 2020-01-01 so the window spans the COVID-2020 crash. This yields 41 factors over 2018-06 → 2026-06 (2,082 rows). The smaller T is then absorbed by shrinkage.

Prices → Returns

Models are built on returns, not price levels, because price series are non-stationary while returns are approximately stationary and comparable across instruments. We use log returns:

r_t = ln( P_t / P_t-1 )

• Time-additivity — multi-day returns are sums, so horizon scaling and aggregation are clean.
• Symmetry — a +10% then −10% round-trip nets to 0 in log space.
• Near-normality — log returns are closer to Gaussian, which the covariance / copula engines exploit.

Safe masking. A return is computed only when both P_t > 0 and P_t-1 > 0; otherwise it is left as NaN rather than producing ±∞.

Calendar & missing-value policy

• Weekend rows removed; U.S. federal holidays marked with USFederalHolidayCalendar.
• Prices reindexed onto a full Monday–Friday business-day grid.
• Bounded forward-fill (limit=3) inside each ticker's own observed history only — no artificial prices before a ticker starts or after it ends.
• A fill diagnostic separates holiday fills from other short-gap fills.

Collinearity & Hierarchical Clustering

We map the dependence structure of daily log-returns to find which factors carry redundant information (collinearity) and which form natural blocks (clusters). Pearson correlation is computed pairwise-complete with a minimum overlap of 120 observations, anchored to a common dense window (the first date on which ≥ 80% of factors are observed) so unequal-history pairs stay comparable.

We cluster on the proper correlation distance

d_ij = √( ½ (1 − ρ̂_ij) ) ∈ [0, 1]

using average-linkage agglomeration, cutting the tree at k = min(8, max(3, ⌊N/5⌋)) clusters. The dendrogram leaf order defines the block ordering of the annotated correlation heatmap. A within-segment pass then quantifies redundancy — near-1 blocks (e.g. the Treasury complex TY1/FV1/TU1) are substitutes and can be pruned; loose blocks carry genuine intra-segment diversification.

Cross-Asset Sector & Crisis Regimes

Diversification is a property of the current regime, not the long run — correlations tend to spike toward 1 in crises ("correlation breakdown"). We re-estimate the sector correlation matrix inside three stress windows and report the delta versus baseline, Δρ̄_ab = ρ̄ ^event − ρ̄^base:

Regime	Window	Shock
GFC 2008	2008-09-01 → 2009-03-31	Global financial crisis / Lehman
COVID 2020	2020-02-15 → 2020-05-31	Pandemic liquidation & policy response
Strait of Hormuz 2025	2025-06-13 → 2025-07-15	Israel–Iran conflict; energy-supply risk

Contract-Roll Stitching (M1 / M2 / M3 → continuous)

A Bloomberg generic front-month (M1) silently switches its underlying contract at each expiry, injecting the calendar basis (M2 − M1) as an artificial jump into returns. Those jumps inflate volatility, distort correlations and poison every downstream covariance estimate. We ingest the 1st/2nd/3rd generic contracts per root, reconcile M1 against the validated clean panel, then build smoothed continuous series with a staggered weighted roll:

r_cont(t) = (1 − w(t)) r_M1(t) + w(t) r_M2(t)

where w ramps 0 → 1 over K days ending s days before expiry. Blending in return space keeps levels positive and returns economically interpretable. The window K and offset s are chosen by grid-search (K ∈ {3,5,7,10} × s ∈ {0,…,10}), minimising a min–max-normalised weighted sum of four costs: residual jump (excess kurtosis), vol distortion, turnover (~1/√K) and tracking carry (~w̄). The parse is cache-aware (Parquet).

Covariance Estimation

The engine needs a covariance matrix Σ. With N assets there are N(N+1)/2 parameters, so the sample covariance is noisy and often ill-conditioned. We compare four base estimators and three structure-aware challengers, then let the data pick.

Base estimators

• Sample — the MLE. Unbiased but high-variance and frequently non-invertible when N ≈ T.
• EWMA (RiskMetrics) — weights recent data by λ^age(t), capturing volatility clustering; we sweep λ ∈ {0.90,…,0.99} (half-lives ≈ 7–69 days).
• EWMA + correlation shrinkage — shrink the EWMA correlation toward the identity, then rescale by EWMA vols: responsive yet regularised.
• Ledoit–Wolf shrinkage — Σ̂ = (1−δ)S + δμI with δ chosen analytically to minimise expected MSE; guarantees a well-conditioned, invertible matrix.

Structure-aware challengers

• Cluster-structured (block) shrinkage — shrink toward the block structure from the clustering step.
• PCA routing — regularise via the leading principal components.
• PCA / Barra factor model — a low-rank factor + idiosyncratic decomposition.

Calibration = rolling out-of-sample backtest

Each candidate is scored by how well it predicts the next day it has not seen, via the one-step Gaussian negative log-likelihood:

NLL_t = ½ [ N ln 2π + ln det Σ̂_t + (x_t+1 − μ_t)^T Σ̂_t^-1 (x_t+1 − μ_t) ]

NLL is a proper scoring rule — it rewards both the right variances (the ln det term) and the right shape (the Mahalanobis term). The lowest-mean-NLL model becomes the baseline: here ewma λ=0.99 (~69-day half-life).

Robust Head-to-Head Selection

All eight estimators are scored on the duplicate-pruned 40-factor universe with a rolling out-of-sample harness that goes well beyond one likelihood:

• NLL — proper scoring of the full predictive density.
• QLIKE — a volatility loss that penalises under-forecasting risk more than over-forecasting — the asymmetry a risk desk wants.
• GMV realised volatility — the economic test: build the global-minimum-variance portfolio from each Σ̂ and measure its realised out-of-sample vol.
• PSD / condition-number gate — every Σ̂ must be positive semi-definite.
• Diebold–Mariano test — is #1 significantly better than #2, or just noise?

Result. The cluster-structured block_shrink target wins all three loss functions — NLL ≈ −134.9 (Diebold–Mariano p ≈ 0.002 vs the runner-up), QLIKE ≈ 1.94, and the lowest GMV realised vol ≈ 0.47% annualised — and all eight pass the PSD gate.

Conditional Gaussian Engine

The core engine conditions the joint Gaussian on the pinned factors and returns the conditional mean of the free factors, with Monte-Carlo draws for the full distribution.

1. Input: user specifies pinned factors A with shocks a; the rest are target factors B.
2. Split: partition μ and Σ into blocks μ_A, μ_B, Σ_AA, Σ_BA.
3. Multicollinearity check: warn if the condition number is high or pinned factors are highly correlated.
4. Condition (Schur complement): compute the conditional mean

E[B | A = a] = μ_B + Σ_BA Σ_AA^-1 (a − μ_A)

5. Sample: draw antithetic Monte-Carlo paths from the conditional law for the full distribution.
6. Uncertainty bands: report 5% / 95% conditional quantiles per factor.
7. Output: a table of implied shocks by factor, with warnings and diagnostics.

Conditional Copula Engines

The Gaussian engine assumes elliptical dependence and thin tails. The copula engines separate the dependence structure from the marginals, using empirical marginals (so each factor's own fat tails and skew are preserved) with either a Gaussian or a Student-t copula for the joint. The Student-t copula adds tail dependence: factors crash together more often than a Gaussian would predict, which is exactly the behaviour a stress engine must capture. A Gaussian-vs-t tail table quantifies the difference in the extremes.

PCA-Space Conditioning (Universe B)

In a single market (the rates curve, a gas-hub complex) factors are highly collinear — the correlation matrix is near-singular, so Σ_AA^-1 is numerically unstable and tiny estimation errors explode into wild peripheral shocks. The fix is to condition in an orthogonal basis.

Diagonalise the standardised covariance Σ = WΛW^T and work in principal-component scores s = W^Tx, which are uncorrelated. For curves the PCs have a clean economic reading: PC1 = level, PC2 = slope, PC3 = curvature. Pinning an original factor to a shock c is a single linear constraint a^Ts = c, and conditioning the diagonal Gaussian on it is closed-form and inversion-free in the collinear directions:

E[s | a^Ts = c] = (Λa / a^TΛa) c ,    x̂ = W · E[s | ·]

Because we never invert the near-singular Σ_AA directly, the result is stable even when the condition number is enormous — demonstrated on the rates complex.

Multicollinearity Checks

Diagnostic layers run before outputting results:

Layer 1: Direction conflicts

Flag when highly correlated pinned factors (|r| > 0.75) are shocked in opposite directions. Note: cross-asset relationships such as equity vs treasury futures are not hard-coded as always opposite — stocks and bonds can move together in inflation, liquidity, or policy-shock regimes — so equity-rate sign conflicts are a warning / regime-review, not an automatic hard inconsistency.

Layer 2: Magnitude extremes

Check if any shock falls in the historical tail (> 4σ soft warning, > 6σ hard warning).

Layer 3: System fragility

If the condition number is very high, the system is fragile — the pseudo-inverse / ridge fallback engages and the user is prompted to re-check pins or switch covariance method.

Out-of-Sample Validation Harness

The engine produces hypothetical scenarios that never happened — so there is no ground truth for them. The only honest test is reconstruction on history: pick real crisis days, pin the factors that actually moved, and ask the engine to predict the rest.

Protocol (no look-ahead)

1. Fit the covariance on the window strictly before day d (last lookback rows).
2. Pin the realised shocks of the anchor factors A = {ES1, CO1, TY1}.
3. Predict the conditional mean of all peripheral factors B with the closed-form Gaussian formula.
4. Score against the realised y_B.

Stress days are the 20 largest joint anchor shocks (max_i |z_i| over A) inside the engine window — the days the model most needs to get right.

Metrics

• RMSE — magnitude error.
• Sign-match % — did we get the direction right?
• Cross-sectional correlation & Kendall-τ — did we rank winners/losers correctly?

We benchmark against a naive baseline that leaves all peripheral factors flat. The engine earns its keep by beating naive on direction and ranking: it cuts RMSE ≈ 30% below the naive baseline, with sign-match ≈ 0.73 and Kendall-τ ≈ 0.41.

Output Contract

The engine implements all three output formats so the desk can choose:

Output	Definition
Conditional mean	μB + ΣBAΣAA-1(a − μA) — the central propagated shock
Mean + adverse add-on	conditional mean plus a tail band (5% / 95% quantile of the conditional law)
Sampled distribution	full Monte-Carlo draws (Gaussian and Student-t copula with empirical margins)

The dashboard exposes the mean and the 5/95 add-on per factor; the CSVs carry the full quantile grid.

API & Interface

Python API

from rbc_scenario.engine import ConditionalGaussianEngine

engine = ConditionalGaussianEngine(
    returns_df=returns_clean,     # safe log returns
    cov_method="ewma",            # baseline: ewma lambda=0.99
    lookback=750,
)

scenario = {
    "ES1 Index":  {"shock": -0.08, "unit": "return"},
    "CL1 Comdty": {"shock": -0.15, "unit": "return"},
    "TY1 Comdty": {"shock":  0.25, "unit": "bp"},
}

result = engine.run(scenario, validate=True, n_paths=20000)
print(result.summary_table())        # conditional mean + 5/95 bands
result.to_csv("MVP/outputs/scenario.csv")

Interactive Jupyter dashboard

An ipywidgets-based UI: search and pin factor shocks, run the Gaussian or Gaussian-copula engine, compare implied target-factor quantiles, and inspect per-factor quantile bars for the selected method.

Reproduction

Kernel: Python 3.11. Dependencies: pandas, numpy, scipy, scikit-learn, seaborn, matplotlib, ipywidgets, openpyxl.

cd RBC
pip install pandas numpy scipy scikit-learn seaborn matplotlib ipywidgets openpyxl
jupyter notebook "RBC_MVP_Aiden_modified (1).ipynb"

Input: data/Products_BB_Reference.xlsx (consolidated), data/Products_BB_Raw.xlsx (authoritative Bloomberg reference), and the multi-contract workbooks data/Products_BB (3) (1).xlsx + data/Second.xlsx (M1/M2/M3 blocks for the roll). Run top-to-bottom; cells share kernel state, and all artifacts are written to MVP/outputs/.

Roadmap

1. Time-varying covariance (DCC-GARCH) — correlations that react to vol regimes; directly fixes the "69-day half-life is too slow for stress" tension.
2. Nonlinear Ledoit–Wolf shrinkage (2020) — drop-in upgrade over linear shrinkage, better-conditioned eigenvalues.
3. Vine / Student-t copulas at the engine core for richer tail dependence; entropy pooling (Meucci) for view-blending.
4. Universe B at scale — extend PCA conditioning to the full gas-hub / rates curve using the roll-stitched continuous contracts.
5. Deep-learning covariance (GNN / Transformer / VAE) — documented stretch goal only; at N ≈ 40, T ≈ 2.7k it overfits and fights PSD constraints.
6. Productionisation — extract versioned engine / data_pipeline / validation modules behind clean interfaces, wrap the dashboard as a Streamlit app, add a hashed scenario registry for audit, unit-test the math, and schedule the data refresh + calibration.