Project 02 · Risk

Monte Carlo risk engine

A portfolio risk system built from scratch: correlated return simulation through a hand-written Cholesky decomposition, VaR and CVaR by Monte Carlo and historical simulation, risk attribution by holding, correlation stress tests, and statistical backtesting.

Fig. 1. Simulated 10-day returns of an equal-weight portfolio of four assets, generated through a Cholesky factor. The tail beyond the 99% VaR is shaded. Illustrative parameters.

(1)Correlated simulation

1.1

Returns to covariance

Prices become daily log returns, \(r_t = \ln(P_t / P_{t-1})\). They add across days, so variance scales linearly with the horizon. Taking the sample covariance of four assets' returns, with Bessel's correction, gives \(\Sigma\): variances on the diagonal, covariances off it.

\[ \Sigma = \frac{1}{T-1} \sum_{t=1}^{T} (r_t - \bar r)(r_t - \bar r)^\top \]

Each term is an outer product: one day's deviation vector times its own transpose, a rank-one matrix. Averaging enough of them recovers the full covariance, which the figure does over ten simulated years.

1.2

Split it apart

\(\Sigma\) mixes two things: how much each asset moves, and how the assets move together. Pulling the volatilities out onto a diagonal matrix \(D\) leaves the correlation matrix \(\rho\), with ones on its diagonal. Watch the entries of \(\Sigma\) turn into correlations.

\[ \Sigma = D\,\rho\,D, \qquad D = \operatorname{diag}(\sigma_1, \dots, \sigma_n), \qquad \Sigma_{ij} = \sigma_i\,\rho_{ij}\,\sigma_j \]

This split comes back in the stress test.

1.3

Factor it, column by column

To simulate, \(\Sigma\) needs a square root. The Cholesky-Banachiewicz algorithm finds the lower-triangular \(L\) with \(LL^\top = \Sigma\), one column at a time: the diagonal entry first, then everything below it. Each entry only uses entries already found: a dot product of two partial rows of \(L\), taken off \(\Sigma_{ij}\), then divided by \(L_{jj}\) (or square-rooted on the diagonal). The figure works through each one. The engine implements this by hand rather than calling a library.

\[ \begin{aligned} L_{jj} &= \sqrt{\Sigma_{jj} - \textstyle\sum_{k<j} L_{jk}^2} \\ L_{ij} &= \frac{1}{L_{jj}}\Big(\Sigma_{ij} - \textstyle\sum_{k<j} L_{ik} L_{jk}\Big), \quad i > j \end{aligned} \]
1.4

Check the factor

Multiplying \(L\) by its transpose rebuilds \(\Sigma\). The difference is at machine precision, which is the engine's first built-in check. The residual on the left is computed live in your browser, in double precision.

\[ \big\lVert LL^\top - \Sigma \big\rVert_\infty \approx 10^{-15} \]

The algorithm costs \(O(n^3/3)\), about twice as fast as LU, because symmetry means only the lower triangle is computed.

1.5

Correlate the noise

Independent standard normals \(Z\) have no structure: a round cloud. Multiplying by \(L\) stretches and shears them into returns with exactly the target covariance, because \(\operatorname{Cov}(LZ) = L\,I\,L^\top = \Sigma\). Over an \(h\)-day horizon:

\[ X = \mu h + L Z \sqrt{h}, \qquad Z \sim \mathcal{N}(0, I) \]

Geometrically, \(L\) maps the unit circle of independent draws onto the covariance ellipse \(\{x : x^\top \Sigma^{-1} x = 1\}\), and each draw follows: the figure sends one \(z\) through \(L\) row by row.

(2)VaR and CVaR

2.1

Build the distribution

Each simulated scenario gives one portfolio return, the weighted sum of its asset returns. Thousands of them build up the distribution of 10-day outcomes.

\[ R_p = w^\top X = \sum_i w_i X_i \]
2.2

Value at risk

VaR at confidence \(\alpha\) is the loss not exceeded with probability \(\alpha\). Slide down to the 1% quantile and read off the loss there.

\[ \operatorname{VaR}_\alpha = -Q_{1-\alpha}(R_p) = -F^{-1}(1-\alpha), \qquad F(x) = \mathbb{P}(R_p \le x) \]

A quantile is the cumulative distribution read backward: stack the histogram into \(F\), find where it reaches 1%, and drop down to the loss.

2.3

Expected shortfall

CVaR averages everything beyond that line, the tail's center of mass, so it measures how bad the tail gets, not just where it starts. It's a coherent risk measure. VaR isn't subadditive in general, so under VaR a diversified portfolio can look riskier than the sum of its parts.

\[ \operatorname{CVaR}_\alpha = -\,\mathbb{E}\big[\,R_p \mid R_p \le -\operatorname{VaR}_\alpha\big] \]

The engine computes both at 95% and 99%, by Monte Carlo and by historical simulation.

(3)Component VaR

3.1

Standalone risk overcounts

Measured one at a time, each holding's VaR ignores how it offsets the others. Their sum overshoots the portfolio's VaR, and the gap is the diversification benefit.

3.2

The gradient

Under a normal model the decomposition has a closed form. Portfolio volatility is \(\sigma_p = \sqrt{w^\top \Sigma_h w}\) over the horizon, and its gradient is one matrix-vector product:

\[ \frac{\partial \sigma_p}{\partial w} = \frac{\Sigma_h w}{\sigma_p}, \qquad \operatorname{Component}_i = z_{\alpha}\, \frac{w_i\, (\Sigma_h w)_i}{\sigma_p} \]

The figure computes \(\Sigma_h w\) row by row, weights it by \(w\), and scales it into each holding's share of VaR.

3.3

Euler decomposition

The engine doesn't assume normality. It estimates marginal VaR from the simulation itself, at the VaR threshold, using a narrow band of \(\pm 0.5\%\,\sigma_p\) around the quantile:

\[ \begin{aligned} \operatorname{MVaR}_i &= \frac{\partial \operatorname{VaR}_p}{\partial w_i} = -\,\mathbb{E}\big[r_i \mid R_p = -\operatorname{VaR}_p\big] \\ \operatorname{CVaR}^{\text{comp}}_i &= w_i \operatorname{MVaR}_i, \qquad \textstyle\sum_i \operatorname{CVaR}^{\text{comp}}_i = \operatorname{VaR}_p \end{aligned} \]

Because VaR is positively homogeneous, Euler's theorem makes the components add up exactly. A position that hedges the portfolio shows a negative contribution.

(4)Stress testing

4.1

Two channels of risk

Splitting the covariance into volatilities \(D\) and correlations \(\rho\) isolates the correlation channel. A stress test can change how the assets move together without touching how much each one moves.

\[ \Sigma_{\text{stress}} = D\,\rho_{\text{stress}}\,D \]
4.2

Correlations converge

In a crisis, assets that normally diversify each other start falling together. The crisis scenario sets every off-diagonal correlation to 0.8, the convergence seen in 2008 and in March 2020, while \(D\) stays exactly as it was.

\[ (\rho_{\text{stress}})_{ij} = 0.8 \quad (i \ne j) \]
4.3

Rebuild the covariance

Multiplying by a diagonal matrix on the left scales rows; on the right, it scales columns. So \(D\rho D\) multiplies each correlation by both assets' volatilities:

\[ (D \rho D)_{ij} = \sigma_i\, \rho_{ij}\, \sigma_j \]

Against the calm matrix, every variance on the diagonal is unchanged and every covariance off it is larger.

4.4

Same volatilities, worse losses

Every asset is exactly as volatile as before, but the cloud collapses toward a line, losses stop offsetting, and VaR rises. Here the starting point is the matrix estimated in section 1; the engine's own calm scenario uses \(\rho = 0.2\) off the diagonal, which the buttons on the figure at the top of the page run.

(5)Backtesting

5.1

Forecast, then wait

Each day's VaR is forecast from the previous 252 days only, so nothing from the future leaks in. The window then rolls forward one day. In this simulated series volatility jumps partway through, and a rolling window is slow to catch up.

5.2

Count the violations

A violation is a day whose loss exceeds that morning's forecast. At 99%, about one day in a hundred should be one. Too many means the model understates risk; too few means it's too conservative.

\[ V = \sum_{t} \mathbf{1}\big\{ r_t < -\operatorname{VaR}_t \big\} \]
5.3

The Kupiec test

The proportion-of-failures test asks whether the violation rate \(\hat p = V/T\) is consistent with the target \(p = 1 - \alpha\). Under the null the statistic follows \(\chi^2\) with one degree of freedom:

\[ \operatorname{LR}_{\text{POF}} = -2 \ln\!\left[ \frac{p^{V} (1-p)^{T-V}}{\hat p^{V} (1-\hat p)^{T-V}} \right] \sim \chi^2_1 \]

The statistic is twice the gap in log-likelihood between the observed rate and the target. The figure measures that gap on the likelihood curve and lays it along the \(\chi^2_1\) axis:

\[ \ell(p) = V \ln p + (T - V)\ln(1-p), \qquad \operatorname{LR}_{\text{POF}} = 2\big[\ell(\hat p) - \ell(p)\big] \]

A p-value above 0.05 means the model isn't rejected. Basel III runs the same kind of test over 250 days.

Here it passes even though most violations bunch up after the jump. Kupiec counts violations but ignores when they happen; an independence test, such as Christoffersen's, is what catches clustering.

(6)Checks

Running python risk_engine.py executes five built-in checks:

  1. The Cholesky residual is about \(10^{-15}\), machine precision.
  2. The simulated covariance converges to its target as the number of paths grows.
  3. \(\operatorname{VaR}_{0.99} \ge \operatorname{VaR}_{0.95}\), and \(\operatorname{CVaR}_\alpha \ge \operatorname{VaR}_\alpha\) at every level.
  4. The crisis scenario (\(\rho = 0.8\)) produces a higher VaR than the calm one (\(\rho = 0.2\)).
  5. Component VaR sums to portfolio VaR to machine precision.

(7)Code

Python, NumPy, pandas, SciPy, yfinance, Streamlit, and Plotly. The dashboard has six tabs: overview, Monte Carlo simulation, VaR and CVaR, stress test, backtesting, and diagnostics (eigenspectrum, Cholesky factor, rolling volatility).