Chapter 17 — Multi-Asset Models and Correlation
"Individual stock prices are noise. Their correlation structure is signal."
After this chapter you will be able to:
- Validate and repair correlation matrices using PSD checks and Higham's alternating projection algorithm
- Simulate correlated multi-asset paths using Cholesky decomposition and explain the geometric meaning of each factor
- Price rainbow options (best-of, worst-of) and basket options using Monte Carlo
- Explain why correlations rise during crises (correlation breakdown) and describe three approaches to stress-testing correlation
- Apply Gaussian and Student's $t$ copula models to multi-asset correlation structures
In a well-diversified portfolio under normal market conditions, the correlations between asset returns stabilise around moderate levels: 0.3 to 0.5 for individual equities within a sector, lower across sectors, lower still across asset classes. But during crises, correlations spike. During the market turmoil of October 2008, intraday correlations between S&P 500 constituent stocks approached 0.9 — portfolios that had been genuinely diversified under normal conditions suddenly behaved as single concentrated bets on the overall market. This correlation breakdown is one of the most dangerous phenomena in multi-asset risk management, and copula models (Chapter 16) were specifically developed to address it.
The mathematical machinery of multi-asset modelling is linear algebra. Correlated asset paths require a correlation matrix that is symmetric positive semi-definite (PSD) — a necessary condition for the covariance matrix to define a valid multivariate distribution. When correlation matrices are estimated from data, they often violate PSD numerically (due to missing data, different observation frequencies, or matrix approximation errors). Higham's nearest-correlation-matrix algorithm projects any symmetric matrix onto the cone of PSD matrices by iteratively alternating projections.
Option pricing over multiple assets confronts the curse of dimensionality: a 10-dimensional PDE has no tractable finite-difference solution. Monte Carlo is the natural tool, using correlated paths generated via Cholesky decomposition. Where speed is critical, analytical approximations based on moment-matching provide fast closed forms at the cost of some accuracy.
17.1 Correlation Fundamentals
For $n$ assets, the correlation matrix $\rho$ must be symmetric positive semi-definite (PSD). Cholesky decomposition $\rho = L L^T$ is used to simulate correlated asset paths.
Why Correlation Matrices Must Be PSD
A correlation matrix that is not PSD would imply a portfolio with negative variance — a mathematical impossibility for any real sum of squared returns. In practice, estimated correlation matrices are often not exactly PSD due to missing data (not all assets trade every day), asynchronous closing prices across time zones, or matrix approximation errors (common when constructing large correlation matrices from factor models). When a correlation matrix fails the PSD check, Higham's (2002) alternating projection algorithm finds the nearest valid matrix in the Frobenius-norm sense.
Cholesky Decomposition: Intuition
Cholesky computes a lower-triangular matrix $L$ such that $\rho = LL^T$. Think of $L$ as a square root of the correlation matrix: each column of $L$ captures one independent source of risk. The first column describes the common factor affecting all assets (systematic risk). The second column — which is orthogonal to the first — captures variation unexplained by the first factor. Each subsequent column adds an independent contribution, orthogonal to all previous ones.
This decomposition makes correlated simulation immediate: if $Z \sim N(0, I_n)$ is a vector of independent standard normals, then $W = LZ$ has covariance $E[WW^T] = L E[ZZ^T] L^T = L I L^T = LL^T = \rho$. The correlations are built into $W$ by construction. For an $n$-asset portfolio, Cholesky reduces an $O(n^2)$ correlated sampling problem to $n$ independent draws plus a matrix-vector multiplication.
Correlation Regimes: The Crisis Problem
The deepest practical challenge of multi-asset modelling is that correlations are not stationary. Under normal market conditions, a diversified equity portfolio might have average pairwise correlations of 0.30–0.45. During the 2008 financial crisis, intraday correlations between S&P 500 constituent stocks approached 0.90 — a portfolio that was genuinely diversified in normal times became essentially a single concentrated market bet at the worst possible moment.
This correlation breakdown phenomenon occurs because in crises, the dominant risk factor becomes a single global factor ("will the financial system survive?") that overwhelms idiosyncratic components. All correlations move towards 1 simultaneously. The Gaussian copula and standard multi-asset GBM models assume constant correlations and therefore dramatically understate tail risk in stress scenarios. Practitioners address this by:
- Stress testing: using historical crisis correlation matrices (2008, 2020, 1987) rather than average correlations
- Regime-switching models: fitting separate correlation matrices for normal and stressed regimes and modelling transitions between them
- Factor models: expressing correlations through a small number of common factors whose loadings can shift; the correlation matrix is then $\rho \approx BB^T + D$ where $B$ is the factor loading matrix and $D$ is idiosyncratic variance
- Copula approaches: using Student's $t$ copula (covered in §17.4) which has positive tail dependence and therefore produces higher correlations when both assets are in extreme tail events
For risk management purposes, always validate model outputs under elevated correlation assumptions. A portfolio that looks diversified at $\rho = 0.3$ may have VaR 40–60% larger than expected at $\rho = 0.8$.
module Correlation = struct
(** Validate that a matrix is symmetric PSD (all eigenvalues ≥ 0).
Uses Gershgorin circles as a quick check, then Cholesky. *)
let is_valid_correlation m =
let n = Array.length m in
(* Check diagonal = 1, off-diagonal in [-1, 1] *)
let basic = ref true in
for i = 0 to n - 1 do
if Float.abs (m.(i).(i) -. 1.0) > 1e-8 then basic := false;
for j = 0 to n - 1 do
if Float.abs m.(i).(j) > 1.0 then basic := false
done
done;
if not !basic then false
else begin
(* Attempt Cholesky — fails if not PSD *)
try
let _ = Numerics.cholesky m in true
with _ -> false
end
(** Nearest correlation matrix by Higham (2002) alternating projections *)
let nearest_correlation ?(tol = 1e-8) ?(max_iter = 1000) m =
let n = Array.length m in
let x = Array.map Array.copy m in (* current iterate *)
let ds = Array.make_matrix n n 0.0 in (* dual variable *)
let converged = ref false in
for _ = 0 to max_iter - 1 do
if not !converged then begin
let r = Array.init n (fun i -> Array.init n (fun j -> x.(i).(j) -. ds.(i).(j))) in
(* Project R onto PSD cone (eigendecomposition, zero out negative eigenvalues) *)
let r_psd = Owl.Mat.(let m = of_arrays r in
let v, d = eig m in
let d' = Mat.map (Float.max 0.0) (Mat.re d) in
to_arrays (Mat.re (v *@ (Mat.diagm d') *@ (Mat.transpose v)))) in
(* Update dual *)
for i = 0 to n - 1 do
for j = 0 to n - 1 do
ds.(i).(j) <- r_psd.(i).(j) -. r.(i).(j)
done
done;
(* Project onto correlation matrices (diagonal = 1) *)
for i = 0 to n - 1 do
for j = 0 to n - 1 do
x.(i).(j) <- if i = j then 1.0 else r_psd.(i).(j)
done
done;
(* Check convergence *)
let diff = ref 0.0 in
for i = 0 to n - 1 do
for j = 0 to n - 1 do
diff := !diff +. (x.(i).(j) -. m.(i).(j)) *. (x.(i).(j) -. m.(i).(j))
done
done;
if !diff < tol then converged := true
end
done;
x
(** Cholesky decomposition of correlation matrix *)
let cholesky_decomp = Numerics.cholesky
(** Generate correlated standard normals: z = L * w, w ~ N(0,I) *)
let correlated_normals chol =
let n = Array.length chol in
let w = Array.init n (fun _ -> Mc.std_normal ()) in
Array.init n (fun i ->
Array.fold_left (fun sum j -> sum +. chol.(i).(j) *. w.(j)) 0.0 (Array.init (i+1) Fun.id)
)
end
17.2 Multi-Asset GBM
For $n$ correlated assets:
$$dS_i = (r - q_i) S_i\cdot dt + \sigma_i S_i\cdot dW_i, \quad \text{Cov}(dW_i, dW_j) = \rho_{ij}\cdot dt$$
Using Cholesky: $dW = L\cdot dZ$ where $Z$ has independent components.
module Multi_gbm = struct
type asset = {
spot : float;
vol : float;
div_yield : float;
}
(** Simulate terminal prices of n correlated assets *)
let simulate_terminal ~assets ~rate ~tau ~corr_chol () =
let n = Array.length assets in
let wvec = Correlation.correlated_normals corr_chol in
Array.init n (fun i ->
let a = assets.(i) in
let drift = (rate -. a.div_yield -. 0.5 *. a.vol *. a.vol) *. tau in
a.spot *. exp (drift +. a.vol *. sqrt tau *. wvec.(i))
)
(** Full path simulation (all time steps) *)
let simulate_paths ~assets ~rate ~tau ~n_steps ~corr_chol () =
let n = Array.length assets in
let dt = tau /. float_of_int n_steps in
let paths = Array.init n (fun i ->
Array.make (n_steps + 1) assets.(i).spot
) in
for t = 0 to n_steps - 1 do
let wvec = Correlation.correlated_normals corr_chol in
for i = 0 to n - 1 do
let a = assets.(i) in
let drift = (rate -. a.div_yield -. 0.5 *. a.vol *. a.vol) *. dt in
paths.(i).(t + 1) <- paths.(i).(t) *. exp (drift +. a.vol *. sqrt dt *. wvec.(i))
done
done;
paths
end
17.3 Rainbow and Basket Options
module Rainbow = struct
(** Best-of call: max(max(S1,S2,...,Sn) - K, 0) *)
let best_of_call ~assets ~rate ~tau ~n_paths ~strike ~corr_chol () =
let df = exp (-. rate *. tau) in
let payoffs = Array.init n_paths (fun _ ->
let terminals = Multi_gbm.simulate_terminal ~assets ~rate ~tau ~corr_chol () in
let max_s = Array.fold_left Float.max (-. Float.max_float) terminals in
Float.max 0.0 (max_s -. strike)
) in
df *. Array.fold_left (+.) 0.0 payoffs /. float_of_int n_paths
(** Worst-of call: max(min(S1,...,Sn) - K, 0) *)
let worst_of_call ~assets ~rate ~tau ~n_paths ~strike ~corr_chol () =
let df = exp (-. rate *. tau) in
let payoffs = Array.init n_paths (fun _ ->
let terminals = Multi_gbm.simulate_terminal ~assets ~rate ~tau ~corr_chol () in
let min_s = Array.fold_left Float.min Float.max_float terminals in
Float.max 0.0 (min_s -. strike)
) in
df *. Array.fold_left (+.) 0.0 payoffs /. float_of_int n_paths
end
module Basket = struct
(** Basket call: max(Σ wᵢ Sᵢ - K, 0) *)
let call ~assets ~weights ~rate ~tau ~n_paths ~strike ~corr_chol () =
let df = exp (-. rate *. tau) in
let payoffs = Array.init n_paths (fun _ ->
let terminals = Multi_gbm.simulate_terminal ~assets ~rate ~tau ~corr_chol () in
let basket = Array.fold_left2 (fun acc w s -> acc +. w *. s) 0.0 weights terminals in
Float.max 0.0 (basket -. strike)
) in
df *. Array.fold_left (+.) 0.0 payoffs /. float_of_int n_paths
(** Moment matching approximation: match basket to single log-normal *)
let call_approx ~assets ~weights ~rate ~tau ~strike ~corr_matrix =
let n = Array.length assets in
(* First moment of basket: sum of forwards *)
let m1 = Array.fold_left2 (fun acc w asset ->
acc +. w *. asset.Multi_gbm.spot
*. exp ((rate -. asset.Multi_gbm.div_yield) *. tau)
) 0.0 weights assets in
(* Second moment: E[B²] = sum_ij wi wj Si Sj exp((ri+rj)T + ρij σi σj T) *)
let m2 = ref 0.0 in
for i = 0 to n - 1 do
for j = 0 to n - 1 do
let ai = assets.(i) and aj = assets.(j) in
m2 := !m2 +. weights.(i) *. weights.(j)
*. ai.Multi_gbm.spot *. aj.Multi_gbm.spot
*. exp ((rate -. ai.Multi_gbm.div_yield +. rate -. aj.Multi_gbm.div_yield) *. tau
+. corr_matrix.(i).(j) *. ai.Multi_gbm.vol *. aj.Multi_gbm.vol *. tau)
done
done;
let basket_vol = sqrt ((log (!m2 /. (m1 *. m1))) /. tau) in
(* Price as log-normal with these moments *)
Black_scholes.call ~spot:m1 ~strike
~rate:0.0 ~div_yield:0.0 ~vol:basket_vol ~tau
end
17.4 Copula Models
Beyond Gaussian, we can use Gumbel, Clayton, or Student-t copulas for non-linear dependence.
module Copula = struct
(** Student-t copula: fatter joint tails *)
let simulate_t_copula ~n ~nu ~corr_chol =
let chi2 = (* sum of nu independent N(0,1)² *)
let s = ref 0.0 in
for _ = 1 to nu do
let z = Mc.std_normal () in
s := !s +. z *. z
done;
sqrt (!s /. float_of_int nu)
in
let normals = Correlation.correlated_normals corr_chol in
(* Student-t marginals: t_i = z_i / sqrt(χ²_ν / ν) *)
let t_vars = Array.map (fun z -> z /. chi2) normals in
(* Transform to uniform via t_ν CDF *)
Array.map (fun t ->
(* Approximate t-CDF, or use Regularized incomplete beta *)
let x = float_of_int nu /. (float_of_int nu +. t *. t) in
let ib = (* regularized incomplete beta — simplified *) 0.5 *. (1.0 +. Numerics.norm_cdf t) in
ignore x;
ib (* Replace with proper Student-t CDF in production *)
) t_vars
(** Gumbel copula for upper-tail dependence *)
let gumbel_bivariate ~u ~v ~theta =
(* C(u,v) = exp(-((−ln u)^θ + (−ln v)^θ)^{1/θ}) *)
let a = (-. log u) ** theta and b = (-. log v) ** theta in
exp (-. (a +. b) ** (1.0 /. theta))
(** Clayton copula for lower-tail dependence *)
let clayton_bivariate ~u ~v ~theta =
(* C(u,v) = max(u^{-θ} + v^{-θ} - 1, 0)^{-1/θ} *)
let s = u ** (-. theta) +. v ** (-. theta) -. 1.0 in
(Float.max 1e-10 s) ** (-. 1.0 /. theta)
end
17.5 Spread and Exchange Options
Spread option: $\max(S_1 - S_2 - K, 0)$. When $K=0$: Margrabe's formula applies exactly.
$$C_{\text{exchange}} = S_1 N(d_1) - S_2 N(d_2)$$
$$d_{1,2} = \frac{\ln(S_1/S_2) \pm \frac{1}{2}\sigma^2 T}{\sigma\sqrt{T}}, \quad \sigma = \sqrt{\sigma_1^2 - 2\rho\sigma_1\sigma_2 + \sigma_2^2}$$
module Margrabe = struct
(** Exchange option: pays max(S1 - S2, 0) at T *)
let exchange_call ~s1 ~s2 ~vol1 ~vol2 ~rho ~tau =
let sigma = sqrt (vol1 *. vol1 -. 2.0 *. rho *. vol1 *. vol2 +. vol2 *. vol2) in
let d1 = (log (s1 /. s2) +. 0.5 *. sigma *. sigma *. tau) /. (sigma *. sqrt tau) in
let d2 = d1 -. sigma *. sqrt tau in
s1 *. Numerics.norm_cdf d1 -. s2 *. Numerics.norm_cdf d2
(** General spread with Kirk's approximation for K > 0 *)
let spread_call_kirk ~s1 ~s2 ~strike ~rate ~vol1 ~vol2 ~rho ~tau =
let df = exp (-. rate *. tau) in
let f1 = s1 *. exp (rate *. tau) in
let f2 = s2 *. exp (rate *. tau) in
let s2k = f2 +. strike /. df in
(* Treat (F2 + K/df) as modified second asset *)
let sigma2k = vol2 *. f2 /. s2k in (* adjusted vol *)
let sigma = sqrt (vol1 *. vol1 -. 2.0 *. rho *. vol1 *. sigma2k +. sigma2k *. sigma2k) in
let d1 = (log (f1 /. s2k) +. 0.5 *. sigma *. sigma *. tau) /. (sigma *. sqrt tau) in
let d2 = d1 -. sigma *. sqrt tau in
df *. (f1 *. Numerics.norm_cdf d1 -. s2k *. Numerics.norm_cdf d2)
end
17.6 Chapter Summary
Multi-asset derivatives and portfolio risk management both depend centrally on the joint distribution of asset returns — and specifically on the correlation structure, which determines how assets move together in both normal and stressed conditions.
The positive semi-definite requirement for correlation matrices is not merely a technical condition: it ensures that no linear combination of assets has negative variance, which is the fundamental consistency requirement of probability theory. In practice, matrices estimated from market data must be cleaned and regularized before use, and Higham's algorithm provides the canonical solution for projection onto the PSD cone.
Correlated Monte Carlo simulation via Cholesky decomposition is the workhorse for multi-asset option pricing. Given $n$ independent standard normal draws $\mathbf{z}$, the correlated vector $\mathbf{x} = L\mathbf{z}$ achieves the target covariance $\Sigma = LL^\top$. This approach scales cleanly to any dimension, limited only by the cost of simulating paths and the need for sufficient samples to converge. For basket options, moment-matching to the first two moments of the basket sum provides fast analytical approximations that are accurate near the money but can fail in the tails.
Copulas separate the marginal distributions of individual assets from their joint dependency structure. The Gaussian copula, ubiquitous in credit modelling, implies symmetric tail dependence. The Clayton copula has stronger lower tail dependence (joint crashes are more likely than joint rallies). Margrabe's formula for exchange options and Kirk's approximation for spread options are exact or near-exact results that avoid Monte Carlo entirely for specific two-asset structures.
Exercises
17.1 Simulate 5-asset basket option prices (equal weights, K=ATM) as correlation increases from 0 to 1. Confirm that at ρ=1 it reduces to a single-asset call.
17.2 Implement Cholesky decomposition from scratch (forward substitution algorithm) and verify against Owl for a 5×5 correlation matrix.
17.3 Validate Kirk's approximation against MC for a spread option $S_1 - S_2 - K$ over a range of strikes K ∈ [-20, +20].
17.4 Implement the Student-t copula properly using the regularized incomplete beta function for the CDF, and compare the joint tail dependence to the Gaussian copula.
Next: Chapter 18 — Market Risk