The source files can be found in examples/.
Covariance estimation
Almost every optimiser reads the covariance matrix, and a short window estimates it badly. With $N$ assets and a little more than $N$ observations, the sample covariance is noisy and ill-conditioned. Mean-variance optimisation inverts it in effect, and the inverse amplifies the noise. Two families of fixes help.
- Denoising separates the signal from the noise in the eigenvalues of the matrix.
Denoisetakes one of three algorithms for the eigenvalues below the noise threshold.FixedDenoisereplaces them with their mean.ShrunkDenoisekeeps the diagonal of their part of the matrix and shrinks its other entries toward zero.SpectralDenoisesets them to zero, which leaves a nearly singular matrix, as section 2 shows. - Sparsification imposes a structure on the inverse.
LoGokeeps only the entries that a network of the assets supports, and it builds the network from a similarity measure such asMaximumDistanceSimilarityorExponentialSimilarity.
You configure both through MatrixProcessing, the mp field of PortfolioOptimisersCovariance. A prior takes that estimator in its ce field.
Reach for denoising or sparsification when your window is short next to the number of assets and you run anything that reads the covariance, such as mean-variance optimisation, risk budgeting or clustering. A lower condition number gives a more stable inverse, and weights that move less when the data changes a little. Compare condition numbers before you settle on a technique.
using PortfolioOptimisers, PrettyTables, LinearAlgebrammtfmt = (v, i, j) -> begin if j == 1 return v else return isa(v, Number) ? "$(round(v, digits=6))" : v endend;resfmt = (v, i, j) -> begin if j == 1 return v else return isa(v, Number) ? "$(round(v*100, digits=3)) %" : v endend;1. The data
We load the daily prices of 20 assets over the last 253 trading days, and convert them to returns.
using CSV, TimeSeries, DataFramesX = TimeArray(CSV.File(joinpath(@__DIR__, "..", "SP500.csv.gz")); timestamp = :Date)[(end - 252):end]rd = prices_to_returns(X)ReturnsResult
nx ┼ 20-element Vector{String}
X ┼ 252×20 Matrix{Float64}
nf ┼ nothing
F ┼ nothing
nb ┼ nothing
B ┼ nothing
ne ┼ nothing
E ┼ nothing
ts ┼ 252-element Vector{Date}
iv ┼ nothing
ivpa ┼ nothing
pnl ┼ AssetPanel
│ pf ┼ Vector{AbstractPanelField}: AbstractPanelField[]
│ amsk ┼ 252×20 AllTrueMask
│ emsk ┴ 252×20 AllTrueMask
2. Covariance estimators
We build one prior per covariance estimator and change only the ce field. The expected returns stay at the sample mean. We compare the sample covariance with three denoisers and two LoGo sparsifications.
ces = ["Vanilla" => PortfolioOptimisersCovariance(), "FixedDenoise" => PortfolioOptimisersCovariance(; mp = MatrixProcessing(; dn = Denoise(; alg = FixedDenoise()))), "ShrunkDenoise" => PortfolioOptimisersCovariance(; mp = MatrixProcessing(; dn = Denoise(; alg = ShrunkDenoise(; alpha = 0.5)))), "SpectralDenoise" => PortfolioOptimisersCovariance(; mp = MatrixProcessing(; dn = Denoise(; alg = SpectralDenoise()))), "LoGo(MaxDist)" => PortfolioOptimisersCovariance(; mp = MatrixProcessing(; alg = LoGo())), "LoGo(ExpDist)" => PortfolioOptimisersCovariance(; mp = MatrixProcessing(; alg = LoGo(; sim = ExponentialSimilarity())))]prs = [k => prior(EmpiricalPrior(; ce = ce), rd) for (k, ce) in ces]6-element Vector{Pair{String, LowOrderPrior{Matrix{Float64}, Nothing, Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}}}:
"Vanilla" => LowOrderPrior
X ┼ 252×20 Matrix{Float64}
o_X ┼ nothing
mu ┼ 20-element Vector{Float64}
sigma ┼ 20×20 Matrix{Float64}
chol ┼ nothing
w ┼ nothing
ens ┼ nothing
kld ┼ nothing
ow ┼ nothing
rr ┼ nothing
fpr ┼ nothing
x_src ┼ nothing
ss_gaps ┼ nothing
pc ┴ nothing
"FixedDenoise" => LowOrderPrior
X ┼ 252×20 Matrix{Float64}
o_X ┼ nothing
mu ┼ 20-element Vector{Float64}
sigma ┼ 20×20 Matrix{Float64}
chol ┼ nothing
w ┼ nothing
ens ┼ nothing
kld ┼ nothing
ow ┼ nothing
rr ┼ nothing
fpr ┼ nothing
x_src ┼ nothing
ss_gaps ┼ nothing
pc ┴ nothing
"ShrunkDenoise" => LowOrderPrior
X ┼ 252×20 Matrix{Float64}
o_X ┼ nothing
mu ┼ 20-element Vector{Float64}
sigma ┼ 20×20 Matrix{Float64}
chol ┼ nothing
w ┼ nothing
ens ┼ nothing
kld ┼ nothing
ow ┼ nothing
rr ┼ nothing
fpr ┼ nothing
x_src ┼ nothing
ss_gaps ┼ nothing
pc ┴ nothing
"SpectralDenoise" => LowOrderPrior
X ┼ 252×20 Matrix{Float64}
o_X ┼ nothing
mu ┼ 20-element Vector{Float64}
sigma ┼ 20×20 Matrix{Float64}
chol ┼ nothing
w ┼ nothing
ens ┼ nothing
kld ┼ nothing
ow ┼ nothing
rr ┼ nothing
fpr ┼ nothing
x_src ┼ nothing
ss_gaps ┼ nothing
pc ┴ nothing
"LoGo(MaxDist)" => LowOrderPrior
X ┼ 252×20 Matrix{Float64}
o_X ┼ nothing
mu ┼ 20-element Vector{Float64}
sigma ┼ 20×20 Matrix{Float64}
chol ┼ nothing
w ┼ nothing
ens ┼ nothing
kld ┼ nothing
ow ┼ nothing
rr ┼ nothing
fpr ┼ nothing
x_src ┼ nothing
ss_gaps ┼ nothing
pc ┴ nothing
"LoGo(ExpDist)" => LowOrderPrior
X ┼ 252×20 Matrix{Float64}
o_X ┼ nothing
mu ┼ 20-element Vector{Float64}
sigma ┼ 20×20 Matrix{Float64}
chol ┼ nothing
w ┼ nothing
ens ┼ nothing
kld ┼ nothing
ow ┼ nothing
rr ┼ nothing
fpr ┼ nothing
x_src ┼ nothing
ss_gaps ┼ nothing
pc ┴ nothing
A lower condition number means a better-posed problem. On this data FixedDenoise lowers it the most, and SpectralDenoise raises it by many orders of magnitude. Measure it on your own data before you trust a technique.
pretty_table(DataFrame(; :estimator => [k for (k, _) in prs], Symbol("cond(sigma)") => [cond(p.sigma) for (_, p) in prs]); formatters = [mmtfmt], title = "Covariance conditioning by estimator") Covariance conditioning by estimator
┌─────────────────┬──────────────────────┐
│ estimator │ cond(sigma) │
│ String │ Float64 │
├─────────────────┼──────────────────────┤
│ Vanilla │ 177.278571 │
│ FixedDenoise │ 80.266747 │
│ ShrunkDenoise │ 118.981311 │
│ SpectralDenoise │ 5.006954228142141e13 │
│ LoGo(MaxDist) │ 157.814013 │
│ LoGo(ExpDist) │ 159.846418 │
└─────────────────┴──────────────────────┘3. Visualising the eigenspectrum
plot_eigenspectrum draws the eigenvalues of a covariance matrix as bars, with the Marchenko-Pastur upper bound $\lambda_+$ of that covariance. A bar above the line is larger than noise alone would give. The denoiser does not use this line. It fits its own threshold to the eigenvalues of the correlation matrix.
using StatsPlots, GraphRecipesThe eigenvalues of the sample covariance.
plot_eigenspectrum(prs[1].second, rd)
The eigenvalues after FixedDenoise.
plot_eigenspectrum(prs[2].second, rd)
The eigenvalues after ShrunkDenoise.
plot_eigenspectrum(prs[3].second, rd)
The eigenvalues after SpectralDenoise.
plot_eigenspectrum(prs[4].second, rd)
The eigenvalues after LoGo sparsification with the maximum distance similarity.
plot_eigenspectrum(prs[5].second, rd)
4. Minimum-variance portfolios
The minimum-variance portfolio reads only the covariance, so it isolates the estimator. We solve it with each prior and compare the weights.
using Clarabelslv = Solver(; name = :clarabel1, solver = Clarabel.Optimizer, settings = Dict("verbose" => false), check_sol = (; allow_local = true, allow_almost = true))ress = [k => optimise(MeanRisk(; r = Variance(), obj = MinimumRisk(), opt = JuMPOptimiser(; pe = p, slv = slv))) for (k, p) in prs]pretty_table(DataFrame(["Assets" => rd.nx; [k => r.w for (k, r) in ress]]); formatters = [resfmt], title = "Minimum-variance weights by covariance estimator") Minimum-variance weights by covariance estimator
┌────────┬──────────┬──────────────┬───────────────┬─────────────────┬───────────────┬───────────────┐
│ Assets │ Vanilla │ FixedDenoise │ ShrunkDenoise │ SpectralDenoise │ LoGo(MaxDist) │ LoGo(ExpDist) │
│ String │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├────────┼──────────┼──────────────┼───────────────┼─────────────────┼───────────────┼───────────────┤
│ AAPL │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ AMD │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ BAC │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ BBY │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ CVX │ 7.432 % │ 8.885 % │ 8.896 % │ 18.188 % │ 7.312 % │ 7.115 % │
│ GE │ 0.806 % │ 0.0 % │ 0.231 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ HD │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 5.834 % │ 4.76 % │
│ JNJ │ 36.974 % │ 33.938 % │ 37.46 % │ 81.812 % │ 32.147 % │ 32.179 % │
│ JPM │ 0.749 % │ 0.0 % │ 0.129 % │ 0.0 % │ 7.035 % │ 6.11 % │
│ KO │ 11.161 % │ 12.876 % │ 10.805 % │ 0.0 % │ 11.262 % │ 8.09 % │
│ LLY │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ MRK │ 17.467 % │ 12.975 % │ 14.07 % │ 0.0 % │ 19.716 % │ 19.314 % │
│ MSFT │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ PEP │ 8.978 % │ 14.197 % │ 9.854 % │ 0.0 % │ 1.468 % │ 0.001 % │
│ PFE │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.001 % │ 0.008 % │
│ PG │ 2.353 % │ 7.674 % │ 4.417 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ RRC │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ UNH │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │ 7.876 % │
│ WMT │ 9.355 % │ 2.528 % │ 9.226 % │ 0.0 % │ 11.158 % │ 10.551 % │
│ XOM │ 4.725 % │ 6.927 % │ 4.911 % │ 0.0 % │ 4.067 % │ 3.995 % │
└────────┴──────────┴──────────────┴───────────────┴─────────────────┴───────────────┴───────────────┘We stack the same weights into one bar per covariance estimator.
plot_stacked_bar_composition([r for (_, r) in ress], rd; xticks = (1:length(ress), [k for (k, _) in ress]))
This page was generated using Literate.jl.