An interactive workbench for portfolio construction: seven allocators on identical inputs, six covariance estimators scored against a known truth, exact constrained mean-variance optimization with per-point optimality certificates, and a Monte Carlo study that separates estimation error from luck.
python3 run.pyNumPy is the only dependency. No build step, no SciPy, no CVXPY — the quadratic programs, the hierarchical clustering, the Marchenko-Pastur fit and the HTTP server are all written directly against NumPy and the standard library.
Mean-variance optimization has a famous failure mode: it is an error maximizer. With N assets there are N(N+1)/2 covariance parameters to estimate; the optimizer inverts that matrix, so it loads up precisely on the directions that are worst estimated. Meanwhile the expected-return vector — the input it is most sensitive to — has a standard error comparable to the estimate itself.
So the interesting questions are not "what is the tangency portfolio". They are:
- How much of an efficient frontier is real? On five years of the built-in dataset, James-Stein shrinkage sets its intensity to 1.00 — it judges the entire cross-sectional spread in average returns to be sampling error. Shrink the means accordingly and the frontier collapses to a single point. The app shows that collapse and explains it rather than hiding it behind a raw sample mean.
- Which allocator actually wins, and can you tell from a backtest? No: a backtest is one sample path and conflates estimation error with luck. So the app ships a Monte Carlo that draws fresh samples from a known distribution, allocates on each, and scores the weights under the true parameters. What comes out is unambiguous — maximum Sharpe has the best oracle and the worst shortfall, with 4× the trial-to-trial spread of risk parity.
- Is your optimizer's answer actually optimal? Every solution here ships with its KKT residual, so optimality is verified rather than asserted. That check caught a hand-rolled Critical Line Algorithm that looked plausible and violated its own weight caps.
Six tabs, all scoping off one dataset and estimator panel in the left rail:
| Tab | What it does |
|---|---|
| Universe | Asset statistics with the t-statistic on each mean, growth curves, correlation heatmap. On synthetic data, the true expected return sits beside the sample estimate. |
| Frontier | The box-constrained efficient frontier with per-point KKT residuals, the capital market line, named reference portfolios, and weight composition along the frontier. |
| Allocators | All seven methods on identical inputs — weights and risk contributions as heat tables, concentration by weight versus by risk, and the HRP clustering tree. |
| Covariance Lab | Eigenvalue spectrum against the Marchenko-Pastur law, raw versus denoised correlations, estimator comparison, and a Monte Carlo accuracy study. |
| Black-Litterman | A view editor (absolute and relative, with confidence sliders) showing prior, posterior, per-asset view impact and the resulting weight change. |
| Backtest | Strict walk-forward with turnover costs, drawdowns, in-sample-versus-realized decay, and the estimation-error Monte Carlo. |
Every optimizer here reduces to repeated Euclidean projection onto
{ w : sum(w) = 1, lo <= w <= hi }
The KKT conditions give w = clip(v − θ, lo, hi) for one scalar θ fixed by the
budget. f(θ) = sum(clip(v − θ, lo, hi)) is piecewise linear and non-increasing,
with kinks only at v_i − hi_i and v_i − lo_i. So rather than bisecting θ to a
tolerance, the code locates the segment between consecutive kinks where f crosses
1 and solves the linear equation there exactly — O(n log n), accurate to 3e-15
against a reference bisection across 6,000 randomized cases including infinite
bounds and short-selling.
This matters: FISTA calls it twice per iteration, thousands of times per solve. The iterative version made the whole frontier roughly 50× slower for no accuracy gain.
min ½w'Σw − λμ'w s.t. 1'w = 1, lo ≤ w ≤ hi, solved by FISTA (Nesterov-accelerated
projected gradient) with adaptive restart on top of that exact projection. The
problem is convex and the projection is exact, so it converges to the global
optimum; sweeping λ ≥ 0 traces the whole frontier.
Every solution returns its KKT residual: the largest violation of
g_i + ν = 0 where lo_i < w_i < hi_i
g_i + ν >= 0 where w_i == lo_i
g_i + ν <= 0 where w_i == hi_i
for g = Σw − λμ, with ν recovered from the interior set. Observed residuals are
~1e-10 to 3e-10 across universes from 6 to 100 assets. Iteration stops on that
certificate rather than on step size, because the certificate is what callers care
about — and the last thousands of FISTA steps move the weights without improving it.
Verification, all in the test suite:
- Random search over 200k long-only portfolios never beats the frontier at matched return.
- Unconstrained minimum variance matches the closed form
Σ⁻¹1/1'Σ⁻¹1to 1e-7. - With a 20% cap: the solver's Sharpe is 0.529 against 0.458 for the best of 400k random capped portfolios, and the cap binds at exactly 0.20000000.
- The residual rejects a merely-feasible point (equal weight scores >1e-3 where the optimum scores <1e-7).
The first version of this used Markowitz's Critical Line Algorithm, which computes the frontier's exact turning points. It was subtly wrong: brute force beat it by 0.34 in volatility, and with a 20% cap it produced weights of 0.2075 and capped Sharpe of 0.334 against 0.458 for random sampling. Rather than ship a plausible-looking "exact" algorithm, it was replaced with the certified QP above. The KKT residual is what made the bug undeniable, and is why every point now carries one.
Seven methods, organized by how much they trust their inputs:
| Allocator | Uses μ | Inverts Σ | Uses Σ |
|---|---|---|---|
| Equal weight (1/N) | no | no | no |
| Inverse volatility | no | no | diagonal only |
| Hierarchical risk parity | no | no | yes |
| Risk parity (ERC) | no | no | yes |
| Minimum variance | no | yes | yes |
| Maximum diversification | no | yes | yes |
| Maximum Sharpe | yes | yes | yes |
| Black-Litterman | yes | yes | yes |
That order is roughly increasing sensitivity to estimation error, and roughly the reverse of their out-of-sample ranking.
Equal risk contribution is the solution of
min ½w'Σw − Σ b_i log(w_i) over w > 0, no budget constraint
The omission is essential. The stationarity condition of the unconstrained
problem is exactly w_i(Σw)_i = b_i — the risk-budget condition — and risk
contributions are scale invariant, so normalizing afterwards lands on the answer.
The first implementation added sum(w) = 1 to the solve, which introduces a
multiplier that corrupts that condition; at annualized covariance scales the
barrier term then dominates and the weights barely move off equal weight. Risk
contributions came out spanning −0.06 to 0.41 instead of all being 0.10.
Each coordinate is now solved exactly: holding the others fixed, the condition in
w_i is the quadratic σ_ii w_i² + (Σw)_i^(−i) w_i − b_i = 0, whose positive root
is closed form. Cyclical coordinate descent on a strictly convex coercive objective
converges monotonically.
Result: risk contributions equal 1/N to 6e-15 for universes up to 80 assets, arbitrary risk budgets are honoured to 1e-14, and under an identity correlation matrix ERC equals inverse-volatility exactly — a structural identity the tests check.
Maximum diversification reuses the same machinery: maximizing (w'σ)/√(w'Σw) is
algebraically a Sharpe maximization with the volatility vector standing in for
expected returns.
HRP (López de Prado 2016) never inverts the covariance matrix, which is the whole point — inversion is where estimation error is amplified, so an allocator that avoids it cannot be an error maximizer. The tests assert this against the source, and check that HRP still produces valid weights on an exactly singular covariance where mean-variance cannot run at all.
Three steps, with SciPy's clustering written out by hand:
- Distance.
d_ij = sqrt(0.5(1 − ρ_ij))— a true metric, so the triangle inequality holds and clustering is meaningful. The common1 − ρis not a metric; the tests verify the triangle inequality on all n³ triples. - Tree. Agglomerative clustering via the Lance-Williams recurrence, covering
single, average, complete and Ward linkage in one formula, emitting SciPy's
(n−1) × 4format so the dendrogram can be drawn. All four recover a planted three-block correlation structure with the blocks contiguous in leaf order. - Recursive bisection. Split the quasi-diagonalized order in half repeatedly, sharing capital inversely to each half's inverse-variance-weighted variance.
Six estimators, all returning annualized, positive-definite matrices:
| Estimator | What it assumes |
|---|---|
sample |
nothing, and it shows |
ledoit_wolf |
optimal linear shrinkage toward constant correlation (LW 2004) |
oas |
Oracle Approximating Shrinkage toward a scaled identity |
ewma |
recent regime is the relevant one |
single_factor |
one market factor plus idiosyncratic noise — 2N+1 parameters instead of N(N+1)/2 |
denoised |
eigenvalues below the Marchenko-Pastur edge are noise |
The Ledoit-Wolf intensity is not a tuning knob: it is the value minimizing expected
squared Frobenius error, estimated from the asymptotic variances of the sample
entries. Implementing it correctly requires the third moment E[y_i³y_j] — the
first attempt used E[y_i²y_j], which is what the expression looks like before
you expand y_i² · y_i y_j.
Denoising fits the Marchenko-Pastur noise variance by grid search against a
hand-written Gaussian KDE of the eigenvalue spectrum, then replaces every
eigenvalue below the upper edge σ²(1 + √(N/T))² with their common average. That
flattens the noise subspace without changing the trace, so no asset's estimated
variance moves — the tests check the trace is preserved to 1e-8.
Scored against the true matrix over 40 trials:
| N, T | T/N | LW Frobenius | Sample Frobenius | LW min-var loss | Sample min-var loss |
|---|---|---|---|---|---|
| 25, 60 | 2.4 | 0.289 | 0.306 | 1.392 | 1.719 |
| 30, 90 | 3.0 | 0.245 | 0.249 | 1.236 | 1.473 |
| 15, 300 | 20 | 0.1241 | 0.1254 | 1.0560 | 1.0569 |
| 10, 504 | 50 | 0.0905 | 0.0917 | 1.0172 | 1.0167 |
The minimum-variance loss is the realized variance of each estimator's min-var portfolio divided by the truly attainable minimum, so 1.0 is perfect. Note how much wider the gap is in that column than in the Frobenius norm: shrinkage fixes exactly the directions an optimizer leans on. At T/N = 20 the two estimators agree to a tenth of a percent — which is why the test for this claim runs at T/N = 2.4, where it is not a coin flip.
Two moves. Reverse-optimize a prior: Π = δΣw_mkt asks what returns would make
the market portfolio optimal, replacing the noisiest input in finance with a
structural identity. Condition on views: a view is P E[R] = Q with uncertainty
Ω, and the posterior is the precision-weighted blend.
Ω follows He-Litterman — each view's variance is the prior variance of the view
portfolio itself, p'(τΣ)p, scaled by stated confidence. Tying Ω to the prior means
a view on a volatile combination is automatically treated as less precise than the
same numeric view on a stable one.
Properties the tests pin down:
- With no views, the posterior is the prior and the optimizer returns the market portfolio to 1e-6.
- A view moves the posterior strictly between prior and target — confidence below 1 never fully overrides.
- Raising confidence monotonically increases the move.
- Assets named in no view still move, by exactly the amount their correlation with the view assets implies. That is what stops one opinion producing a one-asset portfolio.
- A relative view's row sums to zero and widens the spread in the requested direction.
At each rebalance the estimators see only the trailing window and the weights are applied to the following period. The test for this is the one worth copying: it runs the backtest, then replaces every return after the final rebalance with garbage and re-runs. Every weight in the entire path must be bit-identical. An off-by-one in any estimation window fails it immediately.
Costs are charged on one-way turnover, because the optimizers that chase noise trade the most and a cost-free comparison systematically flatters them. On the built-in dataset maximum Sharpe turns over ~0.35/year against 0.17 for equal weight.
A backtest is one sample path. This draws fresh samples from a known distribution, estimates, allocates, and scores the weights under the true parameters — so the spread across trials is estimation error with no market luck in it. Thirty trials of two years each on the built-in dataset:
| Allocator | Oracle Sharpe | Mean achieved | Std across trials | Shortfall |
|---|---|---|---|---|
| Maximum diversification | 1.0814 | 1.0732 | 0.0119 | 0.0082 |
| Hierarchical risk parity | 1.0617 | 1.0712 | 0.0145 | −0.0095 |
| Minimum variance | 1.0624 | 1.0641 | 0.0138 | −0.0017 |
| Maximum Sharpe | 1.1405 | 1.0587 | 0.0275 | 0.0819 |
| Risk parity (ERC) | 1.0455 | 1.0417 | 0.0066 | 0.0038 |
| Inverse volatility | 0.9065 | 0.9062 | 0.0043 | 0.0003 |
| Equal weight (1/N) | 0.8172 | 0.8172 | 0.0000 | 0.0000 |
Maximum Sharpe has the highest ceiling and the largest shortfall, with four times the trial-to-trial spread of risk parity. Equal weight has a spread of exactly zero because it uses no estimates — the test asserts that identity.
One caveat stated plainly: the oracle is a strict ceiling only for maximum Sharpe, which optimizes the metric being scored. Minimum variance minimizes variance, so a sample-based version can post a higher true Sharpe than its own oracle by luck, and the small negative shortfalls above are legitimate rather than a bug. The standard deviation column is the clean cross-allocator comparison. The draws are also Gaussian even when the dataset has fat tails, which makes the table mildly optimistic.
Synthetic (default). A k-factor model with Student-t shocks, and it knows the true μ and Σ. That is what turns every question here from opinion into measurement: allocators can be scored against the genuinely optimal portfolio, covariance estimators against the real matrix. Verified against its own parameters — sample volatility within 0.01 of true at T = 40,000, relative Frobenius error under 0.05, excess kurtosis confirming the fat tails.
Live. Daily closes from Stooq's plain CSV endpoint (no key, no cookie handshake) with Yahoo's chart API as fallback. Tickers are inner-joined on common trading dates rather than forward-filled: filling gaps manufactures zero-return days, which biases volatility down and correlation up — exactly the quantities being estimated.
run.py launcher; loopback IPv4 + IPv6, opens a browser
portlib/
returns.py return panels, annualization, 5 expected-return estimators
covariance.py 6 estimators, Marchenko-Pastur denoising, diagnostics
optimize.py exact projection, FISTA QP with KKT certificates, frontier
portfolios.py 7 allocators behind one interface
hrp.py correlation distance, 4 linkages, recursive bisection
blacklitterman.py equilibrium prior, views, posterior
metrics.py risk contributions, effective breadth, performance stats
backtest.py walk-forward, in/out split, estimation-error Monte Carlo
market.py synthetic factor model with truth; Stooq/Yahoo fetch
server/app.py stdlib JSON API + static files
web/ ES modules, hand-rolled canvas charts
tests/test_portfolio.py 72 tests
The chart layer (web/js/chart.js, util.js, api.js) and the design tokens are
shared lineage with a companion options-pricing project — same palette, same
conventions — extended here with scatter, stacked-area and dendrogram renderers.
python3 -m unittest discover -s tests -v72 tests, ~16 seconds. Properties rather than golden numbers, because a golden number locks in whatever the code did on the day it was written. The ones that earned their keep by catching real defects:
- ERC risk contributions equal 1/N to machine precision (caught the budget-constraint bug).
- The KKT residual rejects a merely-feasible point (caught the broken CLA).
- Perturbing the future leaves the walk-forward weight path bit-identical.
- ERC equals inverse-volatility exactly under identity correlation.
- The correlation distance satisfies the triangle inequality on all n³ triples.
- Denoising preserves the trace, so no asset's variance moves.
- 1/N has exactly zero estimation-error spread.
- Every allocator respects weight bounds to 1e-9 and sums to 1 to 1e-9.
Two tests are deliberately shaped by what they can honestly claim: the shrinkage comparison runs at T/N = 2.4 rather than 20, and the oracle bound is asserted only for the allocator whose objective matches the metric.
- Expected returns are a single vector; there is no factor-based return model or regime switching. Given how the estimation study turns out, that is arguably the right scope — the marginal value of a better mean estimate is large, but so is the marginal risk of overfitting one.
- Constraints are box plus budget. No turnover constraints, sector limits, or transaction-cost-aware optimization; costs are charged in the backtest but not optimized against.
- The backtest rebalances on a fixed calendar, not on drift or signal.
- Stooq and Yahoo are undocumented endpoints that change without notice. Treat the live fetch as a demonstration, not a data pipeline.
- The estimation-error study assumes the true model is stationary, which is the most flattering possible assumption for every method in it.