Monte Carlo without a covariance matrix

Simulated one-day returns for three stocks, generated in your browser from a year of daily closing prices. The returns are drawn from a multivariate normal distribution fitted to the historical data — but the covariance matrixCovariance matrixA table holding, for every pair of assets, how much they tend to move together. With a thousand assets it has a million entries, which is why avoiding it is worth the trouble. is never built.

Why avoid the covariance matrix

The textbook route to correlated normal draws is to estimate the N × N covariance matrix S, factor it, and multiply the factor by a standard normal vector. Both standard factorizations have problems, and they are not the same problems.

CholeskyCholesky decompositionA way of splitting a covariance matrix into a factor times its own transpose, so that multiplying independent random numbers by that factor produces correlated ones. is brittle. The textbook factorization — L lower triangular with a strictly positive diagonal — exists and is unique exactly when S is positive definite. Zero eigenvalues are not fatal in principle: a positive semi-definite matrix still admits S = LL’ with zeros on the diagonal, and a pivoted, rank-revealing variant will find it. The unpivoted routine you actually call will not. Each pivot is computed as a difference, and at rank deficiency that difference is a rounding error of arbitrary sign; take its square root and the factorization fails outright.

That is not a corner case in this application: the sample covariance of T days of N factors has rank at most T, so whenever N exceeds T the matrix is singular by construction. Even short of that it is fragile, because near-linear dependence among factors is the norm rather than the exception — a stock and its ADR, two adjacent tenors on the same yield curve, an index and its largest constituent. Those pairs drive S toward singularity, and a Cholesky routine that survives them at all returns a factor contaminated by rounding.

SVD is robust but not cheap. A singular value decomposition handles rank deficiency gracefully, which is the honest fix for the brittleness above. It does not, however, address either cost that actually binds:

• You must still form S. That is $N^2$ entries of storage and $O(TN^2)$ work to fill them.
• You must still factor it, at $O(N^3)$ in the number of factors.

In a real risk system N is a factor count, not a handful of stocks — every tenor on every curve, every equity index, every exchange rate. At N in the thousands, the cube dominates everything else in the calculation.

The decay factor bounds the other dimension. Exponential weighting discounts day t by λ raised to its age, and RiskMetrics uses λ between 0.94 and 0.97 for daily data. Those correspond to half-lives of about 11 and 23 trading days, and the weight surviving a full year is $\lambda^{250}$ — roughly $2 \times 10^{-7}$ at λ = 0.94, and $5 \times 10^{-4}$ at λ = 0.97. Past a year of history you are adding rows that contribute nothing. So T is not a free parameter that grows with your ambitions; it is pinned near 250.

That is the whole argument for drawing per day instead of per factor. The method below never forms S and never factors it, so the $O(N^2)$ storage and the $O(N^3)$ factorization both disappear. What remains is $O(TN)$ per draw against a T that λ has already capped — linear in the number of factors, with a constant of about 250. Per draw that beats multiplying by a Cholesky factor, which costs $O(N^2)$, once N exceeds T; the setup cost and the storage it beats unconditionally.

Methodology

Start with a matrix R of historical log returnsLog returnThe natural logarithm of today’s price divided by yesterday’s. It behaves better than a percentage change when returns are added up over time.. Each row is a trading day, each column an asset, and each column has had its mean removed. With T days and N assets, R is T × N.

The usual way to simulate correlated normal returns is to estimate the N × N covariance matrix S, factor it as S = LL’ by Cholesky decomposition, and set r = Lz for a standard normal vector z of length N. That works, but it makes you form and factor S first.

You don’t have to. Draw one standard normal per day rather than per asset:

The covariance of $r$ is $R'R$, which is the sample covariance of the historical returns. So $r$ has exactly the distribution you wanted, and $S$ was never constructed — with the consequences for cost and conditioning described above.

This is the approach described in Peter Benson and Peter Zangari, “A general approach to calculating VaR without volatilities and correlations”, RiskMetrics Monitor, Second Quarter 1997.

The same method, implemented six times over — in Python, Go, TypeScript, Rust, C++, and Java — is on GitHub at pbenson/var-without-covariance. All six share one specified random number generator, run against the same price data used on this page, and are tested to produce the same numbers.

Weighting the observations

Equal weighting treats a return from eleven months ago as being as informative as yesterday’s. Exponential weighting instead discounts day t by λ raised to its age in days, so recent market conditions dominate.

Weighting drops into the same expression. Scale row t of R by the square root of its weight w(t), normalized so the weights sum to one:

which gives $r$ covariance $\sum_t w(t)\, r(t) r(t)'$ — the weighted sample covariance. Equal weighting is the special case $w(t) = 1/T$, where the scale factor collapses to the familiar $1/\sqrt{T}$.

What is simulated

Each histogram below shows one-day continuously compounded returns, centered on zero. Dropping the historical drift is the usual convention for one-day risk: over a single day the estimated mean return is far smaller than the noise around it, so including it adds bias without adding information.

The three assets are simulated jointly — every point in the run shares one draw of z — so the correlations among them are preserved even though each panel shows only a marginal distribution.

What to look for

Switch from equal to exponential weighting and watch the histograms breathe: if the last few weeks have been calmer than the year as a whole, the distributions tighten and the 95% value at riskValue at riskA loss level a portfolio should exceed only rarely — a 95% one-day VaR is the loss you expect to be worse than on about one trading day in twenty. pulls in toward zero. That responsiveness is the point of exponential weighting, and it is also its cost — the estimate reacts to recent quiet just as eagerly as to recent turbulence.

Then watch the correlogram under the same switch. Volatility and correlation move independently: a pair can grow more volatile while becoming less correlated, which is why a risk estimate built from volatilities alone can be right about each asset and wrong about the portfolio.

Raise the simulation count and the bars converge on their ticks. That convergence is the claim being tested — sampling error shrinks, and what remains is the historical covariance structure, reproduced exactly, from a matrix that was never built.