Save the Python 3 code below as risk_contribution.py, install NumPy with python3 -m pip install numpy, then run python3 risk_contribution.py. It accepts an already aligned covariance matrix and weights; it neither retrieves market data nor infers a covariance matrix from trade fills. Set periods_per_year=1 for the annual matrix used here. Outputs are decimal volatility, component amounts, percentage shares and the annual covariance used.
The script checks analytic marginal contributions against central differences and prints the reconciliation. Tiny negative eigenvalues within the declared relative tolerance are clipped and reported; materially non-PSD matrices fail.python
import numpy as np
def risk_contributions(weights, covariance, periods_per_year=1.0):
w = np.asarray(weights, dtype=float)
s = np.asarray(covariance, dtype=float)
p = float(periods_per_year)
if w.ndim != 1 or w.size == 0 or s.shape != (w.size, w.size):
raise ValueError("Expected n weights and an n by n covariance")
if not (np.isfinite(w).all() and np.isfinite(s).all()
and np.isfinite(p) and p > 0):
raise ValueError("Missing, infinite or invalid input")
scale = float(np.max(np.abs(s)))
if scale == 0:
raise ValueError("Zero portfolio volatility: contributions undefined")
tol = 1e-10 * scale
if np.max(np.abs(s - s.T)) > tol:
raise ValueError("Covariance must be symmetric")
# Accept only roundoff-sized asymmetry and negative eigenvalues.
s = (s + s.T) / 2
eigenvalues, vectors = np.linalg.eigh(s)
if eigenvalues[0] < -tol:
raise ValueError("Covariance is not positive semidefinite")
repaired = bool(np.any(eigenvalues < 0))
if repaired:
s = (vectors * np.maximum(eigenvalues, 0)) @ vectors.T
with np.errstate(over="raise", invalid="raise", divide="raise"):
s = s * p
q = s @ w
variance = float(w @ q)
numerical_floor = 1e-12 * float(np.linalg.norm(s, 2)) * float(w @ w)
if variance <= numerical_floor:
raise ValueError("Zero/near-zero volatility: attribution unstable")
sigma = np.sqrt(variance)
mrc = q / sigma
crc = w * mrc
share = 100 * crc / sigma
if not all(np.isfinite(x).all() for x in (s, mrc, crc, share)):
raise ValueError("Non-finite calculation")
if not np.isclose(crc.sum(), sigma, rtol=1e-9, atol=0):
raise ArithmeticError("Euler reconciliation failed")
return sigma, mrc, crc, share, s, repaired
S = np.array([[0.04, -0.006], [-0.006, 0.01]])
for weights in ([0.8, 0.2], [0.6, 0.4], [1.2, -0.2]):
sigma, mrc, crc, share, used, repaired = risk_contributions(weights, S)
w = np.array(weights)
# Independent central differences of volatility, without renormalizing.
h = 1e-6
directions = np.eye(len(w))
def vol(x):
return np.sqrt(x @ used @ x)
fd = np.array([(vol(w + h*d) - vol(w - h*d)) / (2*h)
for d in directions])
assert np.allclose(fd, mrc, rtol=1e-7, atol=1e-9)
assert np.isclose(crc.sum(), sigma)
assert np.isclose(share.sum(), 100)
print(weights, "sigma", round(sigma, 10))
print("CRC", np.round(crc, 10), "share %", np.round(share, 6))
print("max derivative error", np.max(np.abs(fd - mrc)),
"roundoff PSD repair", repaired)
assert risk_contributions([0.8, 0.2], S)[2][1] < 0
assert risk_contributions([0.6, 0.4], S)[2][1] > 0
# Singular PSD is valid when this allocation has positive variance.
assert np.isclose(risk_contributions([0.5, 0.5], [[.04, .04], [.04, .04]])[0], .2)
# Cash contributes zero in this model, even though the matrix is singular.
assert risk_contributions([0.8, 0.2], [[.04, 0], [0, 0]])[2][1] == 0
