The source files can be found in examples/.
Subset resampling and cross-validation
The meta-optimiser page builds the three meta-optimisers. This page asks two further questions:
- how far do the out-of-sample predictions of a plain optimiser and of a meta-optimiser move from one fold to the next?
- what does the efficient frontier look like when the optimiser draws random subsets of the universe and averages the weights?
MeanRisk is the benchmark here and SubsetResampling is the meta-optimiser. The prior and the clustering are the ones the meta-optimiser page builds, so you can read the allocations of the two pages against each other.
Reach for subset resampling, and for a meta-optimiser in general, when one fit over every asset moves too far for you to act on. A small change to the estimation window can move the whole allocation, and an average over many drawn universes moves less than one point estimate. Cross-validation on this page measures that movement. It does not build the portfolio you hold.
using PortfolioOptimisers, PrettyTables, StableRNGsresfmt = (v, i, j) -> begin if j == 1 return v else return isa(v, Number) ? "$(round(v * 100, digits = 3)) %" : v endend;1. Data and shared ingredients
We use the data of the meta-optimiser page, one year of daily prices for twenty S&P 500 stocks. We compute the prior and the clustering once, and every cell below uses them.
using CSV, TimeSeries, DataFrames, Clarabel, StatisticsX = TimeArray(CSV.File(joinpath(@__DIR__, "..", "SP500.csv.gz")); timestamp = :Date)[(end - 252):end]rd = prices_to_returns(X)slv = [Solver(; name = :clarabel1, solver = Clarabel.Optimizer, settings = Dict("verbose" => false), check_sol = (; allow_local = true, allow_almost = true)), Solver(; name = :clarabel2, solver = Clarabel.Optimizer, settings = Dict("verbose" => false, "max_step_fraction" => 0.95), check_sol = (; allow_local = true, allow_almost = true)), Solver(; name = :clarabel3, solver = Clarabel.Optimizer, settings = Dict("verbose" => false, "max_step_fraction" => 0.9), check_sol = (; allow_local = true, allow_almost = true))]pr = prior(EmpiricalPrior(), rd)clr = clusterise(ClustersEstimator(; alg = DBHT()), pr.X)jopti = JuMPOptimiser(; pe = pr, slv = slv)jopto = JuMPOptimiser(; slv = slv)JuMPOptimiser
pe ┼ EmpiricalPrior
│ ce ┼ PortfolioOptimisersCovariance
│ │ ce ┼ Covariance
│ │ │ me ┼ SimpleExpectedReturns
│ │ │ │ w ┴ nothing
│ │ │ ce ┼ GeneralCovariance
│ │ │ │ ce ┼ SimpleCovariance: SimpleCovariance(true)
│ │ │ │ w ┴ nothing
│ │ │ alg ┼ FullMoment()
│ │ │ w ┴ nothing
│ │ mp ┼ MatrixProcessing
│ │ │ pdm ┼ Posdef
│ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ dn ┼ nothing
│ │ │ dt ┼ nothing
│ │ │ alg ┼ nothing
│ │ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ me ┼ SimpleExpectedReturns
│ │ w ┴ nothing
│ horizon ┼ nothing
│ fill_limit ┼ nothing
│ x_src ┴ nothing
slv ┼ 3-element Vector{Solver}
│ Solver ⋯
│ Solver ⋯
│ Solver ⋯
wb ┼ WeightBounds
│ lb ┼ Float64: 0.0
│ ub ┴ Float64: 1.0
bgt ┼ Float64: 1.0
sbgt ┼ nothing
gbgt ┼ nothing
lt ┼ nothing
st ┼ nothing
lcse ┼ nothing
cte ┼ nothing
gcarde ┼ nothing
sgcarde ┼ nothing
smtx ┼ nothing
sgmtx ┼ nothing
slt ┼ nothing
sst ┼ nothing
sglt ┼ nothing
sgst ┼ nothing
tn ┼ nothing
fees ┼ nothing
sets ┼ nothing
tr ┼ nothing
ple ┼ nothing
ret ┼ ArithmeticReturn
│ settings ┼ JuMPReturnsSettings
│ │ scale ┼ Float64: 1.0
│ │ lb ┼ nothing
│ │ rte ┼ Bool: true
│ │ fee ┼ Bool: true
│ │ mic ┴ Bool: true
│ ucs ┼ nothing
│ mu ┼ nothing
│ mtx_sqrt ┴ EigenFallbackSquareRoot()
sca ┼ SumScalariser()
ccnt ┼ nothing
cobj ┼ nothing
sc ┼ Int64: 1
so ┼ Int64: 1
ss ┼ nothing
card ┼ nothing
scard ┼ nothing
l2c ┼ nothing
lpc ┼ nothing
linfc ┼ nothing
l1 ┼ nothing
l2 ┼ nothing
lp ┼ nothing
linf ┼ nothing
xbgt ┼ Bool: false
strict ┴ Bool: false
2. Reference allocations
We solve the plain minimum-variance portfolio and the three meta-optimisers. The meta-optimiser page builds the same four. Sections 3 and 4 take the plain optimiser and subset resampling further. Section 3 cross-validates both, and section 4 traces a frontier for each.
res_bench = optimise(MeanRisk(; opt = JuMPOptimiser(; pe = pr, slv = slv)))res_nco = optimise(NestedClustered(; pe = pr, cle = clr, opti = MeanRisk(; obj = MinimumRisk(), opt = jopti), opto = MeanRisk(; obj = MinimumRisk(), opt = jopto)), rd)res_stk = optimise(Stacking(; pe = pr, opti = [MeanRisk(; opt = jopti), HierarchicalRiskParity(; opt = HierarchicalOptimiser(; pe = pr)), InverseVolatility(; pe = pr)], opto = MeanRisk(; obj = MinimumRisk(), opt = jopto)), rd)res_ssr = optimise(SubsetResampling(; pe = pr, opt = MeanRisk(; obj = MinimumRisk(), opt = JuMPOptimiser(; slv = slv)), subset_size = 0.7, n_subsets = 10, rng = StableRNG(123), seed = 42), rd)pretty_table(DataFrame(; :assets => rd.nx, :MinVar => res_bench.w, :NCO => res_nco.w, :Stacking => res_stk.w, :SubsetResampling => res_ssr.w); formatters = [resfmt])┌────────┬──────────┬──────────┬──────────┬──────────────────┐
│ assets │ MinVar │ NCO │ Stacking │ SubsetResampling │
│ String │ Float64 │ Float64 │ Float64 │ Float64 │
├────────┼──────────┼──────────┼──────────┼──────────────────┤
│ AAPL │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ AMD │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ BAC │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ BBY │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ CVX │ 7.432 % │ 10.376 % │ 7.432 % │ 6.786 % │
│ GE │ 0.806 % │ 0.0 % │ 0.806 % │ 0.856 % │
│ HD │ 0.0 % │ 0.0 % │ 0.0 % │ 0.677 % │
│ JNJ │ 36.974 % │ 31.364 % │ 36.974 % │ 20.456 % │
│ JPM │ 0.749 % │ 0.0 % │ 0.749 % │ 1.915 % │
│ KO │ 11.161 % │ 10.028 % │ 11.161 % │ 13.75 % │
│ LLY │ 0.0 % │ 0.0 % │ 0.0 % │ 0.735 % │
│ MRK │ 17.467 % │ 15.876 % │ 17.467 % │ 19.192 % │
│ MSFT │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ PEP │ 8.978 % │ 10.115 % │ 8.978 % │ 10.81 % │
│ PFE │ 0.0 % │ 0.0 % │ 0.0 % │ 2.41 % │
│ PG │ 2.353 % │ 10.609 % │ 2.353 % │ 7.461 % │
│ RRC │ 0.0 % │ 0.0 % │ 0.0 % │ 0.0 % │
│ UNH │ 0.0 % │ 3.483 % │ 0.0 % │ 1.144 % │
│ WMT │ 9.355 % │ 4.863 % │ 9.355 % │ 7.507 % │
│ XOM │ 4.725 % │ 3.285 % │ 4.725 % │ 6.301 % │
└────────┴──────────┴──────────┴──────────┴──────────────────┘Read the three meta-optimiser columns against the plain fit, and read the largest weight in each. SubsetResampling averages over many smaller universes, so its largest weight is smaller than in the other three columns, and it holds more assets. The plot below stacks the four allocations.
using StatsPlots, GraphRecipesplot_stacked_bar_composition([res_bench, res_nco, res_stk, res_ssr], rd)
3. Cross-validation prediction
We score the benchmark and the resampled optimiser under cross-validation. cross_val_predict takes an estimator and returns one out-of-sample prediction per fold, so the two sets of predictions are comparable.
The two optimisers we pass it carry no computed prior. Their JuMPOptimiser holds a solver alone, so the prior is an estimator, the default EmpiricalPrior, and it is fitted again on each training fold. Cross-validation needs it that way. A prior computed once was fitted on every observation, including the ones each test fold holds, so cross-validation refuses it and asks for the estimator. The MeanRisk objectives page covers that distinction.
kfold = KFold(; n = 5)cv_bench = cross_val_predict(MeanRisk(; opt = JuMPOptimiser(; slv = slv)), rd, kfold)cv_ssr = cross_val_predict(SubsetResampling(; opt = MeanRisk(; opt = JuMPOptimiser(; slv = slv)), subset_size = 0.7, n_subsets = 8, rng = StableRNG(123), seed = 42), rd, kfold)scorer = NearestQuantilePrediction(; r = LowOrderMoment(; alg = SecondMoment()))pp_bench = PopulationPredictionResult(; pred = [cv_bench])pp_ssr = PopulationPredictionResult(; pred = [cv_ssr])median_bench = scorer(pp_bench)median_ssr = scorer(pp_ssr)println("MeanRisk cross-validated variance = $(expected_risk(LowOrderMoment(; alg = SecondMoment()), cv_bench))")println("SubsetResampling cross-validated variance = $(expected_risk(LowOrderMoment(; alg = SecondMoment()), cv_ssr))")MeanRisk cross-validated variance = 9.624610341649724e-5
SubsetResampling cross-validated variance = 9.75842830768231e-5Each bar of the first plot is the variance of one test fold of the benchmark.
plot_cv_scores(LowOrderMoment(; alg = SecondMoment()), cv_bench)
Each bar of the second plot is the variance of one test fold of the resampled optimiser.
plot_cv_scores(LowOrderMoment(; alg = SecondMoment()), cv_ssr)
The two printed lines give the variance of the out-of-sample returns of each optimiser, with the five test folds joined into one series. Each plot draws the variance of every fold on its own.
The scorer returns the path closest to the median of the population, so you do not choose a path by hand. The cell below prints the id of the path it returned for each optimiser. A population numbers its paths by position, and each population here holds one cross_val_predict stream, so both print 1. A population of many paths, such as the one CombinatorialCrossValidation returns, gives the id of the path the scorer selected.
println("Median benchmark path id = $(median_bench.id)")println("Median SubsetResampling path id = $(median_ssr.id)")Median benchmark path id = 1
Median SubsetResampling path id = 14. Efficient frontier of a meta-optimiser
The frontier page runs one MeanRisk problem per point. We run the same sweep over the resampled optimiser, so each point of the frontier is an average over drawn universes rather than one fit over every asset. The table prints the largest weight at each point of both frontiers, and the plot draws the two frontiers together.
frontier_ret = ArithmeticReturn(; settings = JuMPReturnsSettings(; lb = Frontier(; N = 15)))mr_front = MeanRisk(; opt = JuMPOptimiser(; pe = pr, slv = slv, ret = frontier_ret))ssr_front = SubsetResampling(; pe = pr, opt = MeanRisk(; opt = JuMPOptimiser(; slv = slv, ret = frontier_ret)), subset_size = 0.7, n_subsets = 8, rng = StableRNG(123), seed = 42)res_mf = optimise(mr_front)res_sf = optimise(ssr_front, rd)rf = factory(Variance(), pr)xs_m = [expected_risk(rf, w, pr.X) for w in res_mf.w]ys_m = [expected_return(ArithmeticReturn(), w, pr) for w in res_mf.w]xs_s = [expected_risk(rf, w, pr.X) for w in res_sf.w]ys_s = [expected_return(ArithmeticReturn(), w, pr) for w in res_sf.w]pretty_table(DataFrame(; :point => 1:length(res_mf.w), :MeanRisk_max_w => [maximum(w) for w in res_mf.w], :SubsetResampling_max_w => [maximum(w) for w in res_sf.w]); formatters = [resfmt])plot(xs_m, ys_m; seriestype = :scatter, marker = (:circle, 5), label = "MeanRisk", xlabel = "Variance", ylabel = "Arithmetic return", title = "Efficient frontiers of MeanRisk and SubsetResampling")plot!(xs_s, ys_s; seriestype = :scatter, marker = (:diamond, 6), label = "SubsetResampling")
5. A covariance fitted again on every subset and fold
Section 3 said that cross-validation refuses a prior computed once. It cannot refuse a matrix pasted into a risk measure the same way. Variance(; sigma = S) is a valid setting, and a matrix you measured somewhere else looks the same as one fitted on the very sample the portfolio is about to be scored on.
A field that a prior fills therefore takes a second form, the estimator that computes the value rather than the value. That is a DeferredQuantity. The optimiser fits it on the prior of each subset and each fold. What reaches the solver is still a plain matrix, and only the moment the matrix is computed has moved.
The two forms give the same answer unless the estimator depends on the universe or the window. Denoising depends on both. FixedDenoise fits a Marchenko-Pastur distribution to the eigenvalues of the correlation matrix, with the ratio of observations to assets, T / N, as its shape. It replaces every eigenvalue at or below the fitted edge with the mean of those eigenvalues. So the 14-asset block of a 20-asset fit is not a 14-asset fit.
ce_dn = PortfolioOptimisersCovariance(; mp = MatrixProcessing(; dn = Denoise(; alg = FixedDenoise())))sigma_full = cov(ce_dn, pr.X)idx = 1:14sigma_refit = cov(ce_dn, view(pr.X, :, idx))sigma_slice = view(sigma_full, idx, idx)println("Largest entry of the full-universe covariance = $(maximum(abs, sigma_full))")println("Largest gap, 14-asset refit against the cut block = $(maximum(abs, sigma_refit .- sigma_slice))")Largest entry of the full-universe covariance = 0.0015659674409100352
Largest gap, 14-asset refit against the cut block = 0.00011166025555061992The second number is the largest gap between the covariance fitted on the 14 assets and the block cut from the 20-asset fit. Read it against the first number, the largest entry of the full covariance.
Inside a resample
SubsetResampling cuts the inner optimiser down to the assets it drew before that optimiser computes its prior. A matrix you pass is cut to those assets. A deferred quantity is not yet a matrix, so it is fitted on the prior of the subset.
ssr_rm = r -> SubsetResampling(; pe = pr, opt = MeanRisk(; obj = MinimumRisk(), r = r, opt = JuMPOptimiser(; slv = slv)), subset_size = 0.7, n_subsets = 10, rng = StableRNG(123), seed = 42)res_pasted = optimise(ssr_rm(Variance(; sigma = sigma_full)), rd)res_deferred = optimise(ssr_rm(Variance(; sigma = ce_dn)), rd)pretty_table(DataFrame(; :assets => rd.nx, :pasted_matrix => res_pasted.w, :deferred_estimator => res_deferred.w); formatters = [resfmt])┌────────┬───────────────┬────────────────────┐
│ assets │ pasted_matrix │ deferred_estimator │
│ String │ Float64 │ Float64 │
├────────┼───────────────┼────────────────────┤
│ AAPL │ 0.0 % │ 0.0 % │
│ AMD │ 0.0 % │ 0.0 % │
│ BAC │ 0.194 % │ 0.229 % │
│ BBY │ 0.0 % │ 0.0 % │
│ CVX │ 7.342 % │ 7.005 % │
│ GE │ 0.044 % │ 0.036 % │
│ HD │ 0.719 % │ 0.691 % │
│ JNJ │ 18.746 % │ 19.666 % │
│ JPM │ 0.708 % │ 0.573 % │
│ KO │ 17.102 % │ 15.623 % │
│ LLY │ 0.039 % │ 0.163 % │
│ MRK │ 16.24 % │ 19.116 % │
│ MSFT │ 0.0 % │ 0.0 % │
│ PEP │ 14.483 % │ 13.508 % │
│ PFE │ 0.299 % │ 0.961 % │
│ PG │ 11.638 % │ 9.076 % │
│ RRC │ 0.0 % │ 0.067 % │
│ UNH │ 0.387 % │ 1.147 % │
│ WMT │ 3.732 % │ 3.512 % │
│ XOM │ 8.325 % │ 8.629 % │
└────────┴───────────────┴────────────────────┘The two columns differ only in how the covariance of each subset was made. The gap between them is the effect of fitting the denoised covariance on the subset instead of cutting it from the full fit.
Inside a fold
We run the same two forms under cross-validation. The pasted matrix was fitted on all 252 observations, so it contains every test fold. The estimator fits on the training fold alone.
cv_pasted = cross_val_predict(MeanRisk(; obj = MinimumRisk(), r = Variance(; sigma = sigma_full), opt = JuMPOptimiser(; slv = slv)), rd, kfold)cv_deferred = cross_val_predict(MeanRisk(; obj = MinimumRisk(), r = Variance(; sigma = ce_dn), opt = JuMPOptimiser(; slv = slv)), rd, kfold)sm = LowOrderMoment(; alg = SecondMoment())println("Pasted-matrix cross-validated variance = $(expected_risk(sm, cv_pasted))")println("Deferred-estimator cross-validated variance = $(expected_risk(sm, cv_deferred))")Pasted-matrix cross-validated variance = 8.908920388744428e-5
Deferred-estimator cross-validated variance = 9.295790278306191e-5The two runs differ in two ways. The pasted matrix was fitted on the observations each test fold holds, so its score is one no portfolio could have earned at the time. The deferred estimator also fits each training fold at a smaller T / N than the 252 rows of the pasted fit give, and T / N is the shape of the denoising fit. Read the deferred estimator's number as the one a portfolio could have earned.
The fields a prior fills, mu, sigma, kt and sk, all work this way. A measure with two or more such fields takes a prior estimator in pe instead, and one fit fills every field you left unstated. The fit fills mu and kt of a Kurtosis, so its pe must compute a cokurtosis, as a HighOrderPriorEstimator does. It fills mu, sigma and chol of a DistributionValueatRisk.
Kurtosis(; pe = HighOrderPriorEstimator())DistributionValueatRisk(; pe = EmpiricalPrior())A field you state by hand keeps the value you gave it. The constructor checks the shape of each field, not whether their values match, and the docstring of each such measure warns you of that. The library does not refuse the setting, because a matrix measured elsewhere is a valid thing to pass.
Summary
A meta-optimiser helps when one fit over every asset moves too far from window to window.
cross_val_predictscores the benchmark and the resampled optimiser out of sample, one prediction per fold.SubsetResamplingaverages many subset solves, and its weights sit closer together than those of one fit.- A frontier sweep runs over a meta-optimiser, so you can read its whole curve against the curve of the plain optimiser rather than one portfolio against another.
- A field a prior fills takes the estimator rather than the value, so it fits again on every subset and every fold instead of holding one answer from the whole sample.
This page was generated using Literate.jl.