Many correlated outcomes: full covariance, factor models and copulas

A risk manager holds all twelve US industries. Twelve good models of one industry each do not answer her questions: what is the 1% worst day for the whole book? and how often do nine or more industries crash on the same day? Both depend on the joint distribution: on how the industries move together, and on whether they move together more in a crash than on an ordinary day.

This notebook builds that joint distribution in three steps, each with an honest failure:

1. a full covariance (LKJCholeskyCov): flexible, but the parameter count grows with the square of the dimension;
2. a factor model: a few common shocks plus industry-specific noise. It has few parameters, but it is not identified, and the symmetries break the sampler's diagnostics in instructive ways;
3. a copula: keep the dependence, but give each industry heavy-tailed (Student-t) margins.

Finally the three are tested out of sample on 2020-2024, a period with a pandemic crash and a bear market in it.

#1 · Data: twelve industries, and a daily volatility filter

The French Data Library assigns every CRSP stock to an industry and publishes daily returns of the resulting portfolios. The 12-industry file has two blocks, value-weighted and then equal-weighted, each introduced by a line of text. We keep the value-weighted block: the rows before the first index entry that is not an 8-digit date.

Daily returns are not independent draws from a fixed distribution. Their volatility clusters: calm months, then weeks in which every day moves 3%. A static joint model fitted to raw returns would average the calm and the storm and be wrong in both. The standard fix, and the one we use, is to divide each return by a volatility forecast made the day before. Here that is RiskMetrics' exponentially weighted moving average,

$$ s_{i,t}^2 = 0.94\, s_{i,t-1}^2 + 0.06\, r_{i,t-1}^2, \qquad z_{i,t} = r_{i,t} / s_{i,t}. $$

The $z_{i,t}$ are "returns in units of yesterday's expected volatility". Everything below models the vector $z_t$ of the 12 industries as independent draws from one joint distribution, and a forecast for day $t$ multiplies back by the (known) $s_{i,t}$. This is the constant-correlation idea of Bollerslev (1990) with the simplest possible volatility model. It is a real limitation, and section 6 will show where it bites.

The filter does its job on the scale (every sd z is close to 1) and removes a good part of the fat tails, but not all of them: the standardised returns still have excess kurtosis of roughly 1-5 (a Normal has 0). Remember that for section 5.

The dependence is strong. Below is the correlation matrix of the training $z$: every pair is positively correlated, half of the pairs above 0.6, and utilities are the odd one out.

#2 · A full covariance matrix with LKJCholeskyCov

The first model is a multivariate Normal with its own mean for each industry and an unrestricted covariance matrix:

$$ z_t \sim \text{MvNormal}(\mu, \Sigma), \qquad \Sigma = \text{diag}(\sigma)\, \Omega\, \text{diag}(\sigma), \qquad \Omega \sim \text{LKJ}(\eta = 2),\ \sigma_i \sim \text{HalfNormal}(1.5). $$

(With predictors in the mean - one regression per industry, correlated errors - this is Zellner's seemingly unrelated regressions. Nothing below changes.)

pm.LKJCholeskyCov samples the Cholesky factor of $\Sigma$ directly, which is what MvNormal wants, and hands back the correlations and standard deviations as by-products.

What does LKJ(2) say in 12 dimensions? The marginal prior of each correlation under LKJ($\eta$) in dimension $d$ is a Beta($\eta - 1 + d/2$, same) stretched to $(-1, 1)$. At $d = 12$ that is Beta(7, 7): centred on zero with a standard deviation of about 0.27. The same $\eta$ means a tighter prior in higher dimension, because a large matrix of big correlations is hard to keep positive definite. Our data (a median correlation of 0.6, from 1258 days) will overrule it, but check the prior against the data before trusting that.

The prior puts little mass where the data are, but it does not rule them out, and 1258 observations of a correlation are a lot of information (the sampling sd of a correlation of 0.7 is about 0.015). We fit.

Clean, fast, and the posterior correlations sit within a few hundredths of the sample correlations: the LKJ prior pulls them slightly towards zero, and the data win.

So why do anything else? Because of how the parameter count grows. A covariance matrix has $d(d+1)/2$ free entries: 78 here, but tens of thousands for a stock universe, and the data per parameter shrink accordingly.

With 500 stocks and five years of daily data (about 1260 days) the full covariance has 125,250 correlations against 630,000 numbers of data, and the sample covariance of 500 series from 1260 days is famously noisy. A factor model asks for $K$ loadings per series instead.

#3 · The factor model and its symmetries

Suppose each day's returns are driven by $K$ common shocks $f_t$ (a "market" shock, a "defensive versus cyclical" shock...) plus an industry-specific one:

$$ z_t = \mu + \Lambda f_t + \varepsilon_t, \qquad f_t \sim \text{Normal}(0, I_K),\quad \varepsilon_t \sim \text{Normal}(0, \text{diag}(\psi^2)). $$

$\Lambda$ is the $d \times K$ matrix of loadings. We never need the 1258 x $K$ latent factors: integrating them out gives

$$ z_t \sim \text{MvNormal}\big(\mu,\ \Lambda\Lambda^\top + \text{diag}(\psi^2)\big), $$

a covariance that is low-rank plus diagonal. NUTS then sees only $dK + 2d$ parameters.

The catch: $\Lambda\Lambda^\top = (\Lambda Q)(\Lambda Q)^\top$ for any orthogonal $K\times K$ matrix $Q$. Flipping the sign of a column, swapping two columns and rotating the columns all leave the likelihood unchanged. With iid Normal priors on the loadings the prior is unchanged too, so the posterior has exactly the same symmetry. The loadings are not identified; only $\Lambda\Lambda^\top$ is.

#3.1 · One factor: the sign flip

With $K = 1$ the only symmetry is $\Lambda \to -\Lambda$. Two mirror-image modes, separated by a wall of low probability (every loading would have to pass through zero at once). Each chain picks one.

The loadings' r_hat is far above 1.01, yet nothing is wrong with the model: the covariance entries have r_hat close to 1, because both modes imply the same $\Sigma$. Every prediction made from this fit is fine. Only the loadings - the thing you would want to interpret - are a mixture of two mirror images.

#3.2 · Three factors: the rotation

With $K = 3$ the symmetry is continuous: the posterior is constant along a whole 3-dimensional family of rotations. There is no wall to cross, so every chain wanders around the orbit.

This failure is quieter than the sign flip. Each chain explores the same rotation orbit, so r_hat is only a little above 1 and nothing is flagged loudly. But every posterior mean of a loading is near zero, and a table of them would say "no industry loads on any factor", which is absurd: the length of Manuf's loading vector (right) is pinned near 1. The sampler also pays for the walk: the loadings' ESS is low, and the fit takes several times longer than the full covariance model of section 2 (compare the two timings printed above).

(A third symmetry, column switching - factors 1 and 2 swapping roles between chains - is the $K > 1$ analogue of label switching in mixtures. It is one of the rotations, so the two fixes below handle it too.)

#3.3 · Fix 1: a lower-triangular constraint, done naively

The textbook fix (Geweke & Zhou 1996; Lopes & West 2004) removes the rotation by forcing $\Lambda$ to be lower triangular with a positive diagonal: industry 1 loads only on factor 1, industry 2 only on factors 1-2, and so on. That leaves $dK - K(K-1)/2$ loadings and, in principle, one mode. Here it is with the industries in the file's order.

Two different failures:

• $K = 1$. Some chains have the NoDur loading at about 0.7 and others have it stuck at zero. The constraint did not delete the mirror mode; it made it a local mode on the boundary. A chain that starts with most loadings negative cannot flip all of them together, so it keeps the other eleven negative and pushes NoDur's loading to the nearest legal value, 0. NoDur's variance is then explained by its own noise $\psi$. That fits worse, but it is a trap all the same.
• $K = 3$. Manufacturing, the third anchor, gets a diagonal loading of about 0.01 in every chain. Its loadings lie almost in the plane spanned by the first two anchors (durables and manufacturing are both cyclical), so the constraint pins factor 3 to a direction that the data barely identify. r_hat is 1.00 here and notices nothing (with another seed, in our tests, one chain found 0.05 instead and r_hat rose to 1.36). The model is "identified" only in name.

The constraint works when the anchors (the first $K$ rows) are well chosen, meaning that their loading vectors are far from collinear, and when the chains start in the right basin.

#3.4 · Fix 1, done properly: anchors from pivoted QR, chains started at the ML solution

1. Find the maximum-likelihood factor solution (15 lines of SciPy below).
2. Choose the anchors with a pivoted QR decomposition of $\hat\Lambda^\top$. It picks the industry with the longest loading vector, then the one most orthogonal to it, and so on. This is the standard way to pick well-conditioned rows.
3. Instead of reordering the data, impose the triangular pattern on the anchor rows wherever they are: anchor $j$ loads on factors $1..j$ only, positively on factor $j$.
4. Rotate $\hat\Lambda$ into that pattern and start every chain there, without nutpie's random jitter (compile_kwargs={"jitter_rvs": set()}).

Now the fit is clean and fast, and all chains agree. Be clear about what was needed, though: a data-driven choice of anchors and chains started in the right basin. With nutpie's default jitter the same model sent some chains into boundary modes in our tests. r_hat cannot flag a wrong basin that every chain shares, and with no jitter all four chains start at one point, so agreement between them is weaker evidence than usual. That is the price of this fix.

#3.5 · Fix 2: sample freely, align afterwards

The other fix leaves the sampler alone and deals with the symmetry after sampling. Since the posterior of the unconstrained model is exactly invariant under $\Lambda \to \Lambda Q$, we may rotate each draw by whatever $Q$ we like without changing the distribution of anything identified. Choose a reference loading matrix $\Lambda_{\text{ref}}$ and rotate each draw to it (an orthogonal Procrustes problem, one SVD per draw):

$$ Q_s = \arg\min_{Q^\top Q = I} \lVert \Lambda_s Q - \Lambda_{\text{ref}} \rVert_F = U V^\top, \quad \text{where } \Lambda_s^\top \Lambda_{\text{ref}} = U D V^\top. $$

For the reference we take the ML loadings after a varimax rotation: the rotation that makes each column have a few large and many small loadings, the easiest to read. The reference only chooses the coordinate system. The uncertainty comes from the draws.

The same draws that looked like noise now have r_hat about 1.00 and an ESS in the thousands. The sampler had been exploring the posterior of the identified quantities perfectly well all along; the loadings only needed a common coordinate system. The unconstrained fit and the anchored fit are also the same model up to rotation (their loadings live in different coordinates). Rotating the anchored draws to the same reference gives the same picture:

Read the columns (the numbers are loadings on standardised returns, so 0.9 means the factor carries most of that industry's daily variance):

• F1, the cyclical market. Manufacturing, finance, durables, business equipment and "other" load at 0.75-0.95. Utilities barely load at all.
• F2, defensives. Utilities, consumer non-durables and chemicals (a sector that includes household products) load strongly, cyclicals weakly. On days when F2 is up, safe havens move together.
• F3, growth and consumer-facing. Business equipment (tech), wholesale/retail and healthcare load at about 0.4-0.5, energy at zero.

Varimax is a choice: any rotation of these three columns fits equally well. "Market / defensive / growth" is a helpful description, not a discovery, and the fact that the posterior supports it is only a statement about this coordinate system.

#4 · How many factors?

$K$ is a modelling choice with a predictive answer. We compare, on the held-out 2020-2024 data, the log predictive density of each day's 12-vector:

$$ \text{lpd}_t = \log \frac{1}{S} \sum_{s=1}^S \text{MvNormal}(z_t \mid \mu_s, \Sigma_s), $$

summed over the 1258 test days. Predictive quantities only need $\Sigma$, which is identified, so we can use the fast anchored fits. We also compute PSIS-LOO on the training days for comparison. Rows of a daily time series are not exchangeable, so LOO here measures "predict a random missing day", not "predict next year".

Four findings, all of them worth stating plainly:

1. Dependence is almost everything. Independent margins (K = 0) are thousands of nats worse than any model with a factor. A single market factor recovers most of the gap.
2. More factors keep helping, and the full covariance is best on held-out data. With 1258 days and only 12 series there is plenty of data for all 66 correlations, and whatever structure a 3-factor model leaves out persists into the test period. K = 3 and 4 are about equally good. At K = 5-6 the factor model is within about two standard errors of the full model.
3. The anchored fits get harder as K grows. K = 4 and 5 printed dozens of divergences and ESS near 100 (the fourth and fifth factors are weak, and their anchors' loadings press against zero). Their held-out numbers are still sensible, since they are averages over draws of an identified $\Sigma$, but they are not fits to publish loadings from.
4. LOO on the training days agrees on the big picture (read the LOO column top to bottom), but at the top it ranks K = 6 slightly above the full model, which the held-out data do not. The difference is about 7 nats with a standard error of about 5: noise, as the comparison below shows. Every model also has a few training days with Pareto $k > 0.7$ (the last column), so the LOO estimates are themselves shaky on the most extreme days. For a time series, the held-out years are the comparison to believe. The difference at the top, with its standard error:

So at $d = 12$ the factor model is a tool for interpretation, not better prediction. When does it earn its keep as a predictor? When there are few days per dimension.

#4.1 · Where factor models win: 30 industries, one year of data

The same library splits the market into 30 industries. Fit only on 2019 (252 days, so 8 days per dimension instead of 105) and predict 2020-2024 again. The full covariance now has 435 correlations to learn from 252 days.

The ranking reverses. With 252 days for 30 series, a factor model with a handful of factors predicts the next five years better than the full covariance, by several standard errors. A likely reason: the full model spends its data on hundreds of correlations it cannot pin down, and its LKJ prior shrinks every one of them towards zero, the wrong direction for industries that all co-move. The factor model shrinks towards a structure instead.

The rule of thumb this illustrates: full covariance when days per dimension are plentiful; factor structure (or a full covariance shrunk towards a factor structure) when they are not. For a universe of hundreds of stocks the second case is the only case.

#5 · Heavy tails: a Gaussian copula with Student-t margins

Every model so far is multivariate Normal, and section 1 showed that the standardised returns are not. A copula separates the two questions: what does each industry's distribution look like (its margin), and how are the industries tied together (the copula)? Sklar's theorem says any joint distribution can be written this way.

The Gaussian copula with Student-t margins:

1. margins: $z_{i,t} \sim \text{Student-t}(\nu_i, \mu_i, \sigma_i)$ with CDF $F_i$;
2. probability integral transform: $u_{i,t} = F_i(z_{i,t})$ is Uniform(0, 1) if the margin is right;
3. normal scores: $x_{i,t} = \Phi^{-1}(u_{i,t})$ is standard Normal;
4. dependence: $x_t \sim \text{MvNormal}(0, R)$ with $R$ a correlation matrix.

The joint density of $z_t$ is the product of the margins times the copula density, and the copula density is the MvNormal of the scores divided by their independent Normal densities:

$$ \log p(z_t) = \sum_i \log f_i(z_{i,t})\; +\; \log \text{MvNormal}(x_t \mid 0, R) - \sum_i \log \text{Normal}(x_{i,t} \mid 0, 1). $$

In PyMC the margins are an ordinary observed StudentT, and the copula term is a pm.Potential. Two numerical details matter:

• pt.ndtri_exp(log_u) computes $\Phi^{-1}(e^{\log u})$ stably from the log-CDF, so a probability of $10^{-12}$ does not underflow.
• The upper tail loses precision ($u = 1 - 10^{-12}$ rounds). Use the symmetry of the t: $x = -\operatorname{sign}(w)\,\Phi^{-1}\big(F(-|w|)\big)$ for the standardised $w = (z - \mu)/\sigma$.

A PyMC 6.3 trap. pm.LKJCorr's documentation says to use it as a correlation matrix. Inside a model, however, its value is the Cholesky factor of the correlation matrix (its default transform, CholeskyCorrTransform, maps to the factor, and its logp expects one). Used as a correlation matrix, it makes nutpie fail with "All initialization points failed". Treat it as the Cholesky factor and form $R = L L^\top$ yourself.

Estimating $\nu$. With $\nu_i$ free inside the copula, each gradient needs the derivative of the incomplete beta function with respect to its parameter. That makes each gradient several times slower. We use the two-stage inference functions for margins approach (Joe 1997): fit the 12 margins on their own (seconds), fix $\nu_i$ at their posterior means, and then fit locations, scales and the copula jointly. Uncertainty about $\nu_i$ is not propagated. Section 7 returns to this.

Tail indices of about 4-9: tails that are clearly heavier than Normal and vary by industry (finance and tech are the heaviest). Now the copula:

A clean fit, but a slow one (see the timing): each gradient evaluates the Student-t CDF, through the incomplete beta function, at 15,096 points.

#5.1 · A contender: the multivariate Student-t

The copula makes each margin heavy-tailed but keeps Gaussian dependence. A Gaussian copula has no tail dependence: however strong the correlation, the probability that industry B crashes, given that A has, tends to zero as the crash gets more extreme. The multivariate Student-t is the opposite compromise. It has one $\nu$ for everybody, but its heavy tails are shared. It is a Normal whose whole covariance is scaled by a common random factor each day (a common volatility shock), so extremes arrive together.

#5.2 · Which joint distribution predicts 2020-2024 best?

• Heavy-tailed margins alone, without dependence, are hopeless. Dependence comes first.
• The Gaussian copula with t margins beats the full multivariate Normal by a few hundred nats. The only change is the margins, so this is what fat tails are worth one industry at a time.
• The multivariate t beats both by a much larger margin, even though its margins are cruder (a single $\nu$). What it adds is shared tails: the days when everything moves a lot at once. The EWMA filter does not remove those days, because it adapts only after the fact. On a day when the whole market is hit by news, every $z_{i,t}$ is large at once.

The tail-dependence picture makes that concrete. For a probability level $q$, take the conditional exceedance $\chi(q) = P(U_j < q \mid U_i < q)$, averaged over all 66 pairs. Under independence it equals $q$. For a Gaussian copula it falls to 0 as $q \to 0$. For a t it levels off at a positive value.