import numpy as np
import matplotlib
matplotlib.use("Agg") # non-interactive backend: render to figures, never pop a window
import matplotlib.pyplot as plt
rng = np.random.default_rng(42) # seed every random draw for reproducibility
np.set_printoptions(precision=4, suppress=True)
print("Setup complete.")nb1.2 — Covariance Matrices & Correlation Heatmaps
In the chapter, Maya learned that a portfolio’s expected return is a boring straight line in the weight w, but its risk is not — risk carries a cross-term, 2w(1-w)\sigma_{12}, and that single term is where diversification comes from. Here we make that claim concrete with code.
We will:
- Reproduce Maya’s three covariance scenarios (assets moving together / uncorrelated / opposed) and watch the cross-term do all the work.
- See correlation as the cosine of an angle between two return vectors.
- Scale up to a small multi-asset covariance matrix \mathbf{\Sigma}, draw a correlation heatmap, confirm \mathbf{\Sigma} is positive semidefinite (\mathbf{w}'\mathbf{\Sigma}\mathbf{w}\ge 0), and find the minimum-variance weights.
The trick to keep in your head the whole way: a covariance matrix is just bookkeeping for one quadratic, \mathbf{w}'\mathbf{\Sigma}\mathbf{w}, and that quadratic is the variance of a portfolio. Everything below is that one fact, viewed from different angles.
1. Maya’s two-asset budget
From the chapter, Maya’s two assets are fixed population quantities:
| Quantity | Asset 1 | Asset 2 |
|---|---|---|
| Expected return \mathbb{E}[R] | 0.08 | 0.04 |
| Volatility \operatorname{sd}(R) | 0.20 | 0.10 |
| Variance \operatorname{Var}(R) | 0.04 | 0.01 |
Asset 1 is the high-octane choice — more expected return, but double the volatility. The one thing left unpinned is how the two assets move together, the covariance \sigma_{12}. She holds weight w in asset 1 and 1-w in asset 2, so
\operatorname{Var}(R_p) = w^2(0.04) + (1-w)^2(0.01) + 2w(1-w)\,\sigma_{12}.
Let’s encode the knowns.
mu = np.array([0.08, 0.04]) # expected returns
sd = np.array([0.20, 0.10]) # volatilities (standard deviations)
var = sd**2 # variances: [0.04, 0.01]
print("variances:", var)
def portfolio_stats(w, cov12):
'Expected return and variance of the two-asset portfolio (w in asset 1).'
w_vec = np.array([w, 1 - w])
er = w_vec @ mu
var_p = w**2 * var[0] + (1 - w)**2 * var[1] + 2 * w * (1 - w) * cov12
return er, var_p
er, vp = portfolio_stats(0.5, 0.0)
print(f"At w=0.5, cov=0: E[R_p]={er:.4f}, Var={vp:.4f}, sd={np.sqrt(vp):.4f}")2. Three scenarios at the even split — the cross-term reveals itself
Hold w = 0.5, so w^2 = (1-w)^2 = 0.25 and 2w(1-w) = 0.5. The only thing we change between scenarios is \sigma_{12}. Watch the expected return stay nailed at 0.06 while volatility slides around — that contrast is the whole point.
- Scenario A — assets move together perfectly: \sigma_{12} = 0.02 (the largest allowed, \rho = +1).
- Scenario B — uncorrelated: \sigma_{12} = 0.
- Scenario C — move oppositely: \sigma_{12} = -0.01 (\rho = -0.5).
print(f"{'scenario':<12}{'cov12':>8}{'E[R_p]':>10}{'Var':>10}{'sd':>9}")
for name, cov12 in [("A together", 0.02), ("B uncorr", 0.0), ("C opposed", -0.01)]:
er, vp = portfolio_stats(0.5, cov12)
print(f"{name:<12}{cov12:>+8.3f}{er:>10.4f}{vp:>10.4f}{np.sqrt(vp):>9.3f}")
# Chapter's worked numbers:
# cov=+0.020 Var=0.0225 sd=0.150 (no diversification: just averages the risk)
# cov=+0.000 Var=0.0125 sd=0.112 (below the average of the two vols, for free)
# cov=-0.010 Var=0.0075 sd=0.087 (below asset 2 alone, the SAFER asset!)Read the last column. Expected return never moved — linearity of expectation does not care how the assets are related. But volatility fell from 0.150 to 0.112 to 0.087 purely because the covariance dropped. In Scenario C the blend (0.087) is less volatile than asset 2 alone (0.10), the safer ingredient. Maya built something calmer than either piece by mixing in the riskier asset. That is diversification, and it lives entirely in the sign of one term.
3. Risk is a parabola in the weight
Maya is not stuck at w = 0.5. Because \operatorname{Var}(R_p) is a quadratic in w, it traces a parabola with a genuine bottom. Calculus gives the minimizing weight
w^\star = \frac{\sigma_2^2 - \sigma_{12}}{\sigma_1^2 + \sigma_2^2 - 2\sigma_{12}}.
Let’s plot the risk curve for all three scenarios and mark each minimum. We sweep w over a grid (allowing a little leverage past [0,1] so the parabola’s shape is visible).
def w_star(cov12):
return (var[1] - cov12) / (var[0] + var[1] - 2 * cov12)
w_grid = np.linspace(-0.2, 1.2, 200)
fig, ax = plt.subplots(figsize=(7, 4.5))
for name, cov12 in [("A: cov=+0.02 (rho=+1)", 0.02),
("B: cov=0 (rho=0)", 0.0),
("C: cov=-0.01 (rho=-0.5)", -0.01)]:
sd_curve = np.array([np.sqrt(portfolio_stats(w, cov12)[1]) for w in w_grid])
line, = ax.plot(w_grid, sd_curve, label=name)
ws = w_star(cov12)
ax.plot(ws, np.sqrt(portfolio_stats(ws, cov12)[1]), "o", color=line.get_color())
ax.axhline(sd[1], ls="--", color="gray", lw=1, label="asset 2 vol = 0.10")
ax.set_xlabel("weight w in asset 1")
ax.set_ylabel("portfolio volatility sd(R_p)")
ax.set_title("Risk is a parabola in w; the bottom is set by the covariance")
ax.legend(fontsize=8)
fig.tight_layout()
fig.savefig("nb12_risk_parabola.png", dpi=110)
plt.close(fig)
for cov12 in (0.02, 0.0, -0.01):
ws = w_star(cov12)
print(f"cov={cov12:+.3f} w*={ws:.3f} min Var={portfolio_stats(ws, cov12)[1]:.4f} "
f"min sd={np.sqrt(portfolio_stats(ws, cov12)[1]):.4f}")
# Scenario B reproduces the chapter: w*=0.200, min Var=0.0080, min sd≈0.0894. Correlation is the cosine of an angle
The chapter’s big reveal: centered random variables behave like vectors, with covariance as the inner product, standard deviation as length, and correlation as the cosine of the angle between them:
\rho_{XY} = \frac{\operatorname{Cov}(X,Y)}{\operatorname{sd}(X)\operatorname{sd}(Y)} = \cos\theta.
Maya’s three scenarios are three angles between an arrow of length 0.20 (asset 1) and an arrow of length 0.10 (asset 2): \rho=+1 is 0^\circ (parallel), \rho=0 is 90^\circ (perpendicular — “uncorrelated” is “orthogonal”), \rho=-0.5 is 120^\circ. Let’s draw them and confirm the angle equals \arccos\rho.
scenarios = [("rho=+1", 0.02), ("rho=0", 0.0), ("rho=-0.5", -0.01)]
fig, axes = plt.subplots(1, 3, figsize=(11, 3.8))
for ax, (label, cov12) in zip(axes, scenarios):
rho = cov12 / (sd[0] * sd[1]) # correlation = cosine
theta = np.arccos(rho) # angle between the two return vectors
v1 = np.array([sd[0], 0.0]) # asset 1: along x-axis, length 0.20
v2 = sd[1] * np.array([np.cos(theta), np.sin(theta)]) # asset 2: length 0.10 at angle theta
for v, c in [(v1, "C0"), (v2, "C1")]:
ax.annotate("", xy=v, xytext=(0, 0),
arrowprops=dict(arrowstyle="->", color=c, lw=2))
ax.set_title(f"{label} (theta={np.degrees(theta):.0f} deg)", fontsize=10)
ax.set_xlim(-0.12, 0.22); ax.set_ylim(-0.05, 0.22)
ax.set_aspect("equal"); ax.axhline(0, color="gray", lw=.5); ax.axvline(0, color="gray", lw=.5)
fig.suptitle("Correlation = cos(angle between centered return vectors)")
fig.tight_layout()
fig.savefig("nb12_correlation_angle.png", dpi=110)
plt.close(fig)
for label, cov12 in scenarios:
rho = cov12 / (sd[0] * sd[1])
print(f"{label:<9} cov={cov12:+.3f} rho={rho:+.2f} angle={np.degrees(np.arccos(rho)):.0f} deg")5. The covariance matrix \mathbf{\Sigma} for two assets
Two assets need three numbers — two variances and one covariance — and we pack them into a symmetric matrix with variances on the diagonal and the covariance off it:
\mathbf{\Sigma} = \begin{pmatrix} 0.04 & \sigma_{12} \\ \sigma_{12} & 0.01 \end{pmatrix}.
The payoff is that any portfolio’s variance becomes one tidy expression, \mathbf{w}'\mathbf{\Sigma}\mathbf{w}. Let’s verify it reproduces the scalar formula exactly for Scenario C.
def sigma_2asset(cov12):
'Two-asset covariance matrix from a given off-diagonal covariance.'
return np.array([[var[0], cov12],
[cov12, var[1]]])
cov12 = -0.01
Sigma2 = sigma_2asset(cov12)
print("Sigma =\n", Sigma2)
w = 0.5
w_vec = np.array([w, 1 - w])
quad_form = w_vec @ Sigma2 @ w_vec # w' Sigma w
scalar = portfolio_stats(w, cov12)[1] # the Section-2 formula
print(f"\nw'Sigma w = {quad_form:.6f}")
print(f"scalar formula = {scalar:.6f}")
assert np.isclose(quad_form, scalar), "matrix and scalar variance must agree"
print("Match: the matrix is just packaging for the same quadratic.")6. Scaling up: a five-asset covariance matrix
Maya’s lender tracks thousands of loans; a real portfolio holds dozens of assets. The matrix machinery does not change — only its size. Let’s simulate five assets with a known correlation structure, draw monthly returns, and estimate \mathbf{\Sigma} from the sample (this is exactly what you do with real return histories: the population probabilities become observed frequencies).
We build a correlation matrix \mathbf{C} by hand so we know the truth, convert it to a covariance matrix via \mathbf{\Sigma} = \mathbf{D}\,\mathbf{C}\,\mathbf{D} where \mathbf{D} = \operatorname{diag}(\text{vols}), then draw returns from a multivariate normal.
names = ["Stocks", "Bonds", "Gold", "Crypto", "REIT"]
vols = np.array([0.18, 0.06, 0.15, 0.60, 0.20]) # annualized-ish vols
# Hand-built TRUE correlation matrix (symmetric, ones on diagonal)
C_true = np.array([
[ 1.00, -0.20, 0.05, 0.45, 0.65], # Stocks
[-0.20, 1.00, 0.10, -0.15, -0.10], # Bonds (the diversifier)
[ 0.05, 0.10, 1.00, 0.20, 0.08], # Gold
[ 0.45, -0.15, 0.20, 1.00, 0.40], # Crypto
[ 0.65, -0.10, 0.08, 0.40, 1.00], # REIT
])
D = np.diag(vols)
Sigma_true = D @ C_true @ D # Sigma = D C D
# Draw a sample of monthly returns and ESTIMATE Sigma from it
n_obs = 600
returns = rng.multivariate_normal(mean=np.zeros(5), cov=Sigma_true, size=n_obs)
Sigma_hat = np.cov(returns, rowvar=False) # each column is an asset -> rowvar=False
print("Estimated covariance matrix (Sigma_hat):\n", Sigma_hat)
print("\nSampled", n_obs, "monthly returns for", len(names), "assets.")7. The correlation heatmap
A covariance matrix mixes units (Crypto’s variance dwarfs Bonds’), so to see relationships we rescale to the correlation matrix: divide each entry by the product of the relevant standard deviations. The diagonal becomes all ones; every off-diagonal entry is a cosine in [-1, 1]. A heatmap is then a grid of cosines — hot where assets move together, cool where they move apart. This is the single most honest one-glance picture of where diversification can and cannot help.
# Correlation from the ESTIMATED covariance: C_ij = Sigma_ij / (sd_i sd_j)
sd_hat = np.sqrt(np.diag(Sigma_hat))
C_hat = Sigma_hat / np.outer(sd_hat, sd_hat)
fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(C_hat, cmap="RdBu_r", vmin=-1, vmax=1)
ax.set_xticks(range(len(names)), names, rotation=45, ha="right")
ax.set_yticks(range(len(names)), names)
for i in range(len(names)): # annotate each cell with its correlation
for j in range(len(names)):
ax.text(j, i, f"{C_hat[i, j]:.2f}", ha="center", va="center",
color="white" if abs(C_hat[i, j]) > 0.5 else "black", fontsize=9)
fig.colorbar(im, ax=ax, label="correlation (cosine of angle)")
ax.set_title("Correlation heatmap (estimated from sampled returns)")
fig.tight_layout()
fig.savefig("nb12_correlation_heatmap.png", dpi=110)
plt.close(fig)
print("Heatmap saved. Note Bonds (negative correlations) is the natural diversifier.")
print("Recovered correlations are close to the truth, e.g. Stocks-REIT:",
f"{C_hat[0, 4]:.2f} vs true 0.65")8. Positive semidefiniteness: no portfolio can have negative variance
Because \mathbf{w}'\mathbf{\Sigma}\mathbf{w} is a variance, it can never be negative — for any weight vector \mathbf{w}. That property is called positive semidefiniteness, and it is equivalent to “all eigenvalues of \mathbf{\Sigma} are \ge 0.” Let’s check both: throw 10,000 random weight vectors at \widehat{\mathbf{\Sigma}} and confirm none produce a negative variance, then look at the eigenvalues directly.
# Test 1: random weights should never give negative variance
W = rng.normal(size=(10_000, 5)) # 10k arbitrary weight vectors
quad_forms = np.einsum("ij,jk,ik->i", W, Sigma_hat, W) # w' Sigma w for each row
print(f"min of w'Sigma w over 10,000 random w: {quad_forms.min():.6e} (>= 0 required)")
assert quad_forms.min() >= 0, "found a negative variance -- impossible for a valid Sigma"
# Test 2: eigenvalues all non-negative
eigvals = np.linalg.eigvalsh(Sigma_hat) # eigvalsh: symmetric matrices
print("eigenvalues:", eigvals)
print("all eigenvalues >= 0:", np.all(eigvals >= -1e-12))
print("smallest eigenvalue:", f"{eigvals.min():.6e}")9. The minimum-variance portfolio
With many assets, the weight that minimizes \mathbf{w}'\mathbf{\Sigma}\mathbf{w} subject to the budget constraint \mathbf{1}'\mathbf{w} = 1 (weights sum to 100%) has a clean closed form:
\mathbf{w}^\star = \frac{\mathbf{\Sigma}^{-1}\mathbf{1}}{\mathbf{1}'\mathbf{\Sigma}^{-1}\mathbf{1}}.
This needs \mathbf{\Sigma}^{-1} to exist — which it does precisely when \mathbf{\Sigma} is non-singular (no asset is a perfect linear combination of the others). Here is Maya’s experiment at scale: find the lowest-risk blend of all five assets.
ones = np.ones(5)
Sigma_inv = np.linalg.inv(Sigma_hat)
w_mv = Sigma_inv @ ones / (ones @ Sigma_inv @ ones) # min-variance weights
print("Minimum-variance weights:")
for nm, wi in zip(names, w_mv):
print(f" {nm:<8}{wi:>8.3f}")
print(f" {'sum':<8}{w_mv.sum():>8.3f}")
mv_var = w_mv @ Sigma_hat @ w_mv
print(f"\nMin-variance portfolio: variance={mv_var:.5f}, volatility={np.sqrt(mv_var):.4f}")
print("Individual asset volatilities:", dict(zip(names, np.round(sd_hat, 3))))
print("The blend is calmer than any single asset -- diversification, many-asset edition.")Your Turn
You now have all the machinery. Pick at least one:
Add a third asset to Maya’s two-asset world. Suppose a third asset has volatility 0.12 and is uncorrelated with both originals. Build the 3\times3 covariance matrix \mathbf{\Sigma} by hand, then use the closed form \mathbf{w}^\star = \mathbf{\Sigma}^{-1}\mathbf{1} / (\mathbf{1}'\mathbf{\Sigma}^{-1}\mathbf{1}) to find the minimum-variance blend. Is it more diversified than the two-asset answer?
Find the min-variance portfolio that allows no short-selling (every weight \ge 0). The closed form can hand back negative weights; clip them to zero, renormalize, and compare the volatility. (Real constrained optimization waits in a later week — this is a first taste.)
Break \mathbf{\Sigma} on purpose. Make a “clone” asset that is 0.999 correlated with Stocks and rebuild \mathbf{\Sigma}. Print its eigenvalues and its condition number (
np.linalg.cond). Watch the smallest eigenvalue collapse toward zero and the condition number explode — the matrix is becoming near-singular, andnp.linalg.invstarts producing wild, unstable weights. This is a preview of collinearity in Week 2: when regressors are nearly linearly dependent, the same matrix wall makes OLS estimates blow up.
A scaffold for option 3 to get you started:
# Your Turn -- option 3 scaffold: near-singular Sigma (preview of Week 2 collinearity)
clone_corr = 0.999
C_break = np.array([
[1.00, clone_corr], # Stocks
[clone_corr, 1.00], # near-clone of Stocks
])
D_break = np.diag([0.18, 0.18])
Sigma_break = D_break @ C_break @ D_break
eig = np.linalg.eigvalsh(Sigma_break)
print("eigenvalues:", eig)
print(f"condition number: {np.linalg.cond(Sigma_break):.1f} (huge => near-singular)")
print("smallest eigenvalue is nearly zero -> Sigma barely positive definite.")
# Try it: bump clone_corr to 1.0 exactly and np.linalg.inv will raise / give garbage.
# That same wall is what kills OLS under perfect collinearity in Week 2.