Case study: joint densities and Monte Carlo integration

The model

Take the annual returns $R_A$ and $R_B$ of the two assets as normal, with the numbers used since case study 3.2.5. Their expected returns are 6% and 8%, their volatilities 12% and 18%, and their correlation $\rho = 0.3$.

The numbers are illustrative, and normal annual returns are a simplification. They are used here because the answer can then be checked exactly.

Each asset alone is a one-variable question (lesson 2.6.3). A loses more than 10% when $R_A < -0.10$. Since $\frac{-0.10 - 0.06}{0.12} = -1.3333$, that is 1.3333 standard deviations below its mean, so the chance is $N(-1.3333) = 0.0912$. For B it is $N(-1) = 0.1587$. Their product is 0.01447.

The product is right only if the two returns are independent: if one tells you nothing about the other. They are correlated, so the question needs both returns at once.

A joint density is a surface with volume 1

Start with two independent standard normal draws $Z_1$ and $Z_2$, as case study 3.7.4 drew them. Their joint density is the surface

with $\varphi$ the standard normal density. Lesson 4.7.3 found that its volume over the whole plane is 1. The probability that the pair lands in a region is the volume over that region.

Two regions show the idea. Both draws are negative on the quarter plane $z_1 < 0$, $z_2 < 0$, a rectangle, so the product rule of lesson 4.7.1 gives $\frac{1}{2} \times \frac{1}{2} = \frac{1}{4}$. The pair lands within distance 2 of the origin with probability

in polar coordinates (lesson 4.7.3). That is 0.8647, the share of draws that the dashed circle of radius 2 in case study 3.7.4 should hold.

Correlated returns: the density by a change of variables

Case study 3.7.4 made correlated returns from independent draws with the Cholesky factor $L$ of the covariance matrix (lesson 3.7.2): $\mathbf{x} = \boldsymbol{\mu} + L\mathbf{z}$. For two assets,

This map multiplies every area by $\det L = \sigma_A\sigma_B\sqrt{1 - \rho^2} = 0.020605$ (lesson 4.7.2). A small cell keeps its probability as it moves, so its density, probability per unit area, is divided by $\det L$.

The peak, at the expected returns, is $\frac{1}{2\pi \det L} = 7.7241$ per unit of return squared. The units make sense in numbers. Both returns land within half a percentage point of their means, a square of side 0.01, with probability about $7.7241 \times 0.01^2 = 0.000772$; dblquad gives 0.000772 too.

The probability of a joint loss

Both lose more than 10% on the corner $x < -0.10$, $y < -0.10$. Its probability is the volume over the corner:

Both limits are infinite, as in lesson 2.6.3, and dblquad accepts -np.inf. Cutting the corner off 8 standard deviations below the means gives the same 0.028810: the tails beyond hold almost nothing.

Correlation does not change either asset alone: A still loses more than 10% in 9.12% of years. It moves volume towards the corners where both are low or both are high. The corner's probability rises with $\rho$:

At $\rho = 0$ the double integral is the product, 0.0145, as lesson 4.7.1's product rule says.

A region between curves

The regions need not be rectangles. A portfolio of half A and half B loses more than 10% when $0.5x + 0.5y < -0.10$, that is when $y < -0.2 - x$. The inner limit depends on $x$, as in lesson 4.7.1, and dblquad gives 0.082139.

The answer can be checked another way. The portfolio's return is itself normal (Level 5), with mean 7% and, by case study 1.2.4's formula, volatility 12.2229%. So the chance is $N\left(\frac{-0.10 - 0.07}{0.122229}\right) = 0.082139$.

The cut-off $-1.5$ is more than 8 standard deviations below both means, since $0.06 - 8 \times 0.12 = -0.9$ and $0.08 - 8 \times 0.18 = -1.36$.

Monte Carlo: the same probability by simulation

There is a second way to find a volume: simulate. Draw many years of returns with $\mathbf{x} = \boldsymbol{\mu} + L\mathbf{z}$, as in case study 3.7.4, and count the share that land in the corner. That share estimates the probability, and the estimate is called Monte Carlo integration.

With the seed 47, here is the share among the first $N$ of a million draws, against the double integral's 0.028810:

The standard error $\sqrt{\hat{p}(1 - \hat{p})/N}$, with $\hat{p}$ the estimate, is the typical size of the error; module 5.7 explains where it comes from. Each error here is within two standard errors. A Monte Carlo number should always be reported with its standard error.

The box version works for any integral. Lesson 4.7.1's volume under $x^2 + 2y^2$ over $[0, 2] \times [0, 1]$ is 4.

With 1,000 uniform points and the seed 49, the area 2 times the mean height gives 3.8932, with standard error 0.0837. With 100,000 points it gives 4.0011, with standard error 0.0084.

Error like one over the square root of $N$

One run can be lucky or unlucky, so run the simulation 200 times at each size and measure the spread of the estimates. The spread falls from 0.01896 at 100 draws to 0.00052 at 100,000, close to the formula $\sqrt{p(1 - p)/N}$ each time.

On log scales the points lie on a line of slope $-\frac{1}{2}$, read as in lesson 2.8.2, where Simpson's rule had slope $-4$. Case study 3.7.4 met the same pattern for estimated correlations: ten times the draws narrowed the band only about $\sqrt{10} \approx 3.16$ times.

Rare events cost more. The relative error is $\sqrt{(1 - p)/(pN)}$: with 10,000 draws it is 5.8% for this corner. For a 1-in-1,000 event it is 31.6%, and a 10% relative error needs about 99,900 draws.

dblquad against simulation, and where Monte Carlo wins

In two dimensions the double integral wins easily. With the corner cut off 8 standard deviations out, dblquad gives 0.028810 to about eight digits from about 700 values of the density. A million simulated years give only about three digits, $0.02892 \pm 0.00017$.

More dimensions turn this round. A grid with $n$ points along each of $d$ axes needs $n^d$ values. With $n = 100$, that is $10^4$ in two dimensions, $10^6$ in three and $10^{20}$ in ten.

With a fixed budget of $N$ values, a midpoint grid's error shrinks like $N^{-2/d}$, while Monte Carlo's shrinks like $N^{-1/2}$ in any number of dimensions.

The two meet near four dimensions; beyond that, simulation wins. The rows give orders of size only, since each method's error also has a constant in front.

Three assets already show the cost. The chance that A, B and C (case study 3.2.5) all lose more than 10% is a triple integral.

SciPy's tplquad finds 0.013030, but needs about 29,000 values of the density instead of 700. A million simulated years give $0.013038 \pm 0.00011$ at the same cost as before. The independence product, 0.0029, is more than four times too small.

import math

from scipy.stats import norm

MU = (0.06, 0.08)
SIGMA = (0.12, 0.18)


def test_case_study():
    """Assets A and B with correlation 0.3 both lose more than 10% with probability about 0.028810"""
    expect_close(joint_loss_exact(MU, SIGMA, 0.3, -0.10), 0.028809628207340565, what="joint_loss_exact(A, B, 0.3, -0.10)", rel=1e-6)


def test_independent():
    """With correlation 0 the answer is the product of the single chances, about 0.014471"""
    product = norm.cdf(-0.16 / 0.12) * norm.cdf(-0.18 / 0.18)
    expect_close(joint_loss_exact(MU, SIGMA, 0.0, -0.10), product, what="joint_loss_exact with rho = 0", rel=1e-6)


def test_negative_correlation():
    """With correlation −0.3 a joint loss is rarer: about 0.004982"""
    expect_close(joint_loss_exact(MU, SIGMA, -0.3, -0.10), 0.004982408097745086, what="joint_loss_exact with rho = -0.3", rel=1e-6)


def test_standard_normals():
    """Standard normals with correlation 0.5 are both negative a third of the time"""
    expect_close(joint_loss_exact((0, 0), (1, 1), 0.5, 0.0), 1 / 3, what="joint_loss_exact((0, 0), (1, 1), 0.5, 0)", rel=1e-6)


def test_daily():
    """Daily returns (means 0.04% and 0.03%, volatilities 1% and 1.5%, correlation 0.6) both fall below −2% with probability about 0.011227"""
    expect_close(joint_loss_exact((0.0004, 0.0003), (0.01, 0.015), 0.6, -0.02), 0.011226545532604399, what="joint_loss_exact for daily returns", rel=1e-5)


def test_mc_close():
    """200,000 simulated pairs land within 4 standard errors of the exact 0.028810"""
    p, se = joint_loss_mc(MU, SIGMA, 0.3, -0.10, 200_000, seed=1)
    assert abs(p - 0.028809628207340565) < 4 * se, f"estimate {p:.5f} is more than 4 standard errors ({se:.5f}) from 0.028810"


def test_mc_standard_error():
    """The standard error is √(p(1 − p)/n) for the estimate p"""
    p, se = joint_loss_mc(MU, SIGMA, 0.3, -0.10, 50_000, seed=2)
    expect_close(se, math.sqrt(p * (1 - p) / 50_000), what="standard error", rel=1e-9)


def test_mc_seed():
    """The same seed gives the same estimate; another seed gives a different one"""
    a = joint_loss_mc(MU, SIGMA, 0.3, -0.10, 10_000, seed=5)
    b = joint_loss_mc(MU, SIGMA, 0.3, -0.10, 10_000, seed=5)
    c = joint_loss_mc(MU, SIGMA, 0.3, -0.10, 10_000, seed=6)
    expect_equal(a, b, what="two runs with seed 5")
    assert a != c, "seeds 5 and 6 gave the same estimate: is the seed used?"


def test_bad_inputs():
    """A volatility of 0, a correlation of 1 or n = 0 raises ValueError"""
    calls = [
        lambda: joint_loss_exact(MU, (0.0, 0.18), 0.3, -0.10),
        lambda: joint_loss_exact(MU, SIGMA, 1.0, -0.10),
        lambda: joint_loss_mc(MU, SIGMA, -1.0, -0.10, 1000, seed=1),
        lambda: joint_loss_mc(MU, SIGMA, 0.3, -0.10, 0, seed=1),
    ]
    for k, call in enumerate(calls, 1):
        try:
            call()
        except ValueError:
            continue
        raise AssertionError(f"case {k}: expected ValueError")