The source files can be found in examples/.
Uncertainty sets
The view priors of the earlier pages, Black-Litterman, entropy pooling and opinion pooling, change the moments themselves. Robust optimisation keeps the estimated moments and assumes that they are wrong by some amount. You choose a region around the estimate that the true covariance or mean is likely to lie in, the uncertainty set. The optimiser then minimises the worst-case risk, or maximises the worst-case return, over that region. The risk it reports is an upper bound for every covariance inside the set, the estimate included.
An uncertainty set comes from an estimator, such as NormalUncertaintySet, and an algorithm. BoxUncertaintySetAlgorithm gives an interval for each entry of the moment, and EllipsoidalUncertaintySetAlgorithm gives one ellipsoid over all the entries. sigma_ucs builds a covariance set and mu_ucs builds a mean set. A robust risk measure such as UncertaintySetVariance takes the covariance set, and ArithmeticReturn takes the mean set to give a worst-case expected return. This page covers both sets, and the confidence level that sets their size.
Reach for uncertainty sets when the estimated moments are noisy, for example over a short window, and you want weights that allow for that noise more than you want to act on a forecast. A box set bounds each entry on its own. An ellipsoidal set bounds all the entries together. On this page the ellipsoidal set gives the less concentrated weights, for the covariance in section 4 and for the mean in section 5. If you hold views about where the moments are going, reach for a view prior instead. You can also use both, because a robust risk measure works with any prior.
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. The returns data
We use the same S&P 500 slice as the other examples.
using CSV, TimeSeries, DataFramesX = TimeArray(CSV.File(joinpath(@__DIR__, "..", "SP500.csv.gz")); timestamp = :Date)[(end - 252):end]rd = prices_to_returns(X)pr = prior(EmpiricalPrior(), rd)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
2. Building covariance uncertainty sets
We build two covariance uncertainty sets with a NormalUncertaintySet, a box and an ellipsoid. The box draws random covariance matrices from the Wishart law that normal returns imply, so we fix the random number generator to make the run repeatable. At its default the ellipsoid takes its radius from a chi-square quantile and draws nothing.
ucs_box = sigma_ucs(NormalUncertaintySet(; pe = EmpiricalPrior(), rng = StableRNG(1), alg = BoxUncertaintySetAlgorithm()), rd.X)ucs_ell = sigma_ucs(NormalUncertaintySet(; pe = EmpiricalPrior(), rng = StableRNG(1), alg = EllipsoidalUncertaintySetAlgorithm()), rd.X)EllipsoidalUncertaintySet
sigma ┼ 400×400 LinearAlgebra.Diagonal{Float64, Vector{Float64}}
k ┼ Float64: 21.15732657569969
class ┼ SigmaUncertaintySetClass()
val ┴ 20×20 Matrix{Float64}
3. The confidence level q sizes the set
NormalUncertaintySet takes a confidence level q, whose default is 0.05, and q sets the size of the set. A smaller q covers a more extreme worst case, so it gives a larger set. For a box set every interval widens, and for an ellipsoidal set the radius grows. We build a box set at four values of q and print the total width of each box, which is the sum of the lengths of its intervals.
qs = [0.01, 0.05, 0.10, 0.20]box_widths = [let u = sigma_ucs(NormalUncertaintySet(; rng = StableRNG(1), q = q, alg = BoxUncertaintySetAlgorithm()), rd.X) sum(abs, u.ub .- u.lb) end for q in qs]pretty_table(DataFrame(; q = qs, Symbol("box total width") => box_widths); title = "Total width of the box set at each q")Total width of the box set at each q
┌─────────┬─────────────────┐
│ q │ box total width │
│ Float64 │ Float64 │
├─────────┼─────────────────┤
│ 0.01 │ 0.0576656 │
│ 0.05 │ 0.0440054 │
│ 0.1 │ 0.0370207 │
│ 0.2 │ 0.028841 │
└─────────┴─────────────────┘4. Robust and nominal minimum variance
UncertaintySetVariance is the robust form of Variance. It minimises the worst-case variance over the uncertainty set instead of the variance at the point estimate. We solve the nominal minimum-variance portfolio, then the robust one with the box set and with the ellipsoidal set.
using Clarabelslv = Solver(; name = :clarabel1, solver = Clarabel.Optimizer, settings = Dict("verbose" => false), check_sol = (; allow_local = true, allow_almost = true))res_nom = optimise(MeanRisk(; r = Variance(), obj = MinimumRisk(), opt = JuMPOptimiser(; pe = pr, slv = slv)))res_box = optimise(MeanRisk(; r = UncertaintySetVariance(; ucs = ucs_box), obj = MinimumRisk(), opt = JuMPOptimiser(; pe = pr, slv = slv)))res_ell = optimise(MeanRisk(; r = UncertaintySetVariance(; ucs = ucs_ell), obj = MinimumRisk(), opt = JuMPOptimiser(; pe = pr, slv = slv)))MeanRiskResult
jr ┼ JuMPOptimisationResult
│ pa ┼ ProcessedJuMPOptimiserAttributes
│ │ pr ┼ 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
│ │ wb ┼ WeightBounds
│ │ │ lb ┼ 20-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
│ │ │ ub ┴ 20-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
│ │ lt ┼ nothing
│ │ st ┼ nothing
│ │ lcsr ┼ nothing
│ │ ctr ┼ nothing
│ │ gcardr ┼ nothing
│ │ sgcardr ┼ nothing
│ │ smtx ┼ nothing
│ │ sgmtx ┼ nothing
│ │ slt ┼ nothing
│ │ sst ┼ nothing
│ │ sglt ┼ nothing
│ │ sgst ┼ nothing
│ │ tn ┼ nothing
│ │ fees ┼ nothing
│ │ plr ┼ nothing
│ │ ret ┼ ArithmeticReturn
│ │ │ settings ┼ JuMPReturnsSettings
│ │ │ │ scale ┼ Float64: 1.0
│ │ │ │ lb ┼ nothing
│ │ │ │ rte ┼ Bool: true
│ │ │ │ fee ┼ Bool: true
│ │ │ │ mic ┴ Bool: true
│ │ │ ucs ┼ nothing
│ │ │ mu ┼ 20-element Vector{Float64}
│ │ │ mtx_sqrt ┴ EigenFallbackSquareRoot()
│ │ sca ┼ SumScalariser()
│ │ imsk ┴ nothing
│ retcode ┼ OptimisationSuccess
│ │ res ┴ Dict{Any, Any}: Dict{Any, Any}()
│ sol ┼ JuMPOptimisationSolution
│ │ w ┴ 20-element Vector{Float64}
│ model ┼ A JuMP Model
│ │ ├ solver: Clarabel
│ │ ├ objective_sense: MIN_SENSE
│ │ │ └ objective_function_type: AffExpr
│ │ ├ num_variables: 441
│ │ ├ num_constraints: 6
│ │ │ ├ AffExpr in MOI.EqualTo{Float64}: 1
│ │ │ ├ Vector{AffExpr} in MOI.Nonnegatives: 1
│ │ │ ├ Vector{AffExpr} in MOI.Nonpositives: 1
│ │ │ ├ Vector{AffExpr} in MOI.SecondOrderCone: 1
│ │ │ ├ Vector{AffExpr} in MOI.PositiveSemidefiniteConeTriangle: 1
│ │ │ └ Vector{AffExpr} in MOI.PositiveSemidefiniteConeSquare: 1
│ │ └ Names registered in the model
│ │ └ :E, :M, :M_PSD, :T, :W, :WpE, :bgt, :ceucs_variance, :eucs_variance_risk_1, :ge_soc1, :k, :lw, :obj_expr, :ret, :ret_1, :ret_vec, :risk, :risk_minimised, :risk_vec, :sc, :so, :t_eucs1, :variance_flag, :w, :w_gross_ub, :w_lb, :w_ub, :x_eucs1
r ┼ UncertaintySetVariance
│ settings ┼ RiskMeasureSettings
│ │ scale ┼ Float64: 1.0
│ │ ub ┼ nothing
│ │ rke ┴ Bool: true
│ ucs ┼ EllipsoidalUncertaintySet
│ │ sigma ┼ 400×400 LinearAlgebra.Diagonal{Float64, Vector{Float64}}
│ │ k ┼ Float64: 21.15732657569969
│ │ class ┼ SigmaUncertaintySetClass()
│ │ val ┴ 20×20 Matrix{Float64}
│ sigma ┼ 20×20 Matrix{Float64}
│ mtx_sqrt ┼ EigenFallbackSquareRoot()
│ ratio_alg ┼ SquareRootRatio()
│ corner ┴ nothing
fb ┴ nothing
The table puts the weights of the three portfolios side by side. On this data the ellipsoidal set gives the less concentrated weights.
pretty_table(DataFrame(["Assets" => rd.nx, "Nominal" => res_nom.w, "Box-robust" => res_box.w, "Ellipsoid-robust" => res_ell.w]); formatters = [resfmt], title = "Minimum-variance weights, nominal and robust") Minimum-variance weights, nominal and robust
┌────────┬──────────┬────────────┬──────────────────┐
│ Assets │ Nominal │ Box-robust │ Ellipsoid-robust │
│ String │ Float64 │ Float64 │ Float64 │
├────────┼──────────┼────────────┼──────────────────┤
│ AAPL │ 0.0 % │ 0.0 % │ 0.0 % │
│ AMD │ 0.0 % │ 0.0 % │ 0.0 % │
│ BAC │ 0.0 % │ 0.0 % │ 1.595 % │
│ BBY │ 0.0 % │ 0.0 % │ 0.0 % │
│ CVX │ 7.432 % │ 7.168 % │ 5.578 % │
│ GE │ 0.806 % │ 0.0 % │ 2.045 % │
│ HD │ 0.0 % │ 0.0 % │ 1.998 % │
│ JNJ │ 36.974 % │ 41.083 % │ 18.906 % │
│ JPM │ 0.749 % │ 0.0 % │ 2.822 % │
│ KO │ 11.161 % │ 14.02 % │ 10.551 % │
│ LLY │ 0.0 % │ 0.0 % │ 2.94 % │
│ MRK │ 17.467 % │ 17.04 % │ 14.114 % │
│ MSFT │ 0.0 % │ 0.0 % │ 0.13 % │
│ PEP │ 8.978 % │ 10.924 % │ 10.694 % │
│ PFE │ 0.0 % │ 0.0 % │ 3.76 % │
│ PG │ 2.353 % │ 0.001 % │ 7.888 % │
│ RRC │ 0.0 % │ 0.0 % │ 0.0 % │
│ UNH │ 0.0 % │ 0.0 % │ 4.343 % │
│ WMT │ 9.355 % │ 7.207 % │ 7.899 % │
│ XOM │ 4.725 % │ 2.557 % │ 4.736 % │
└────────┴──────────┴────────────┴──────────────────┘The composition of the nominal, box-robust and ellipsoid-robust portfolios.
using StatsPlots, GraphRecipesplot_stacked_bar_composition([res_nom, res_box, res_ell], rd; xticks = (1:3, ["Nominal", "Box", "Ellipsoid"]))
5. Robust expected returns with a worst-case mean
You can make the mean robust too. Build a mean uncertainty set with mu_ucs and pass it to ArithmeticReturn through its ucs keyword. The optimiser then maximises the worst-case expected return over the set instead of the return at the point estimate. A sample mean over one year of daily returns carries a large error, and an objective that seeks return acts on that error unless the mean is robust.
Pass ArithmeticReturn the mean set that mu_ucs returns, as UncertaintySetVariance takes the set that sigma_ucs returns. If you pass the estimator instead, the optimiser builds the set when it solves, and it then needs the returns data as well.
rf = 4.2 / 100 / 252mu_box = mu_ucs(NormalUncertaintySet(; pe = EmpiricalPrior(), rng = StableRNG(1), alg = BoxUncertaintySetAlgorithm()), rd)mu_ell = mu_ucs(NormalUncertaintySet(; pe = EmpiricalPrior(), rng = StableRNG(1), alg = EllipsoidalUncertaintySetAlgorithm()), rd)ret_nom = optimise(MeanRisk(; obj = MaximumRatio(; rf = rf), opt = JuMPOptimiser(; pe = pr, slv = slv)))ret_box = optimise(MeanRisk(; obj = MaximumRatio(; rf = rf), opt = JuMPOptimiser(; pe = pr, slv = slv, ret = ArithmeticReturn(; ucs = mu_box))))ret_ell = optimise(MeanRisk(; obj = MaximumRatio(; rf = rf), opt = JuMPOptimiser(; pe = pr, slv = slv, ret = ArithmeticReturn(; ucs = mu_ell))))MeanRiskResult
jr ┼ JuMPOptimisationResult
│ pa ┼ ProcessedJuMPOptimiserAttributes
│ │ pr ┼ 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
│ │ wb ┼ WeightBounds
│ │ │ lb ┼ 20-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
│ │ │ ub ┴ 20-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
│ │ lt ┼ nothing
│ │ st ┼ nothing
│ │ lcsr ┼ nothing
│ │ ctr ┼ nothing
│ │ gcardr ┼ nothing
│ │ sgcardr ┼ nothing
│ │ smtx ┼ nothing
│ │ sgmtx ┼ nothing
│ │ slt ┼ nothing
│ │ sst ┼ nothing
│ │ sglt ┼ nothing
│ │ sgst ┼ nothing
│ │ tn ┼ nothing
│ │ fees ┼ nothing
│ │ plr ┼ nothing
│ │ ret ┼ ArithmeticReturn
│ │ │ settings ┼ JuMPReturnsSettings
│ │ │ │ scale ┼ Float64: 1.0
│ │ │ │ lb ┼ nothing
│ │ │ │ rte ┼ Bool: true
│ │ │ │ fee ┼ Bool: true
│ │ │ │ mic ┴ Bool: true
│ │ │ ucs ┼ EllipsoidalUncertaintySet
│ │ │ │ sigma ┼ 20×20 LinearAlgebra.Diagonal{Float64, Vector{Float64}}
│ │ │ │ k ┼ Float64: 5.604501123581913
│ │ │ │ class ┼ MuUncertaintySetClass()
│ │ │ │ val ┴ 20-element Vector{Float64}
│ │ │ mu ┼ 20-element Vector{Float64}
│ │ │ mtx_sqrt ┴ EigenFallbackSquareRoot()
│ │ sca ┼ SumScalariser()
│ │ imsk ┴ nothing
│ retcode ┼ OptimisationSuccess
│ │ res ┴ Dict{Any, Any}: Dict{Any, Any}()
│ sol ┼ JuMPOptimisationSolution
│ │ w ┴ 20-element Vector{Float64}
│ model ┼ A JuMP Model
│ │ ├ solver: Clarabel
│ │ ├ objective_sense: MAX_SENSE
│ │ │ └ objective_function_type: AffExpr
│ │ ├ num_variables: 23
│ │ ├ num_constraints: 7
│ │ │ ├ AffExpr in MOI.EqualTo{Float64}: 1
│ │ │ ├ QuadExpr in MOI.LessThan{Float64}: 1
│ │ │ ├ Vector{AffExpr} in MOI.Nonnegatives: 1
│ │ │ ├ Vector{AffExpr} in MOI.Nonpositives: 1
│ │ │ ├ Vector{AffExpr} in MOI.SecondOrderCone: 2
│ │ │ └ VariableRef in MOI.GreaterThan{Float64}: 1
│ │ └ Names registered in the model
│ │ └ :G, :T, :bgt, :cdev_soc_1, :dev_1, :eucs_ret_1, :k, :lw, :obj_expr, :ohf, :ret, :ret_1, :ret_vec, :risk, :risk_vec, :sc, :sdp_corner, :so, :sr_risk, :t_eucs_gw_1, :variance_flag, :variance_risk_1, :w, :w_gross_ub, :w_lb, :w_ub, :x_eucs_w_1
r ┼ Variance
│ settings ┼ RiskMeasureSettings
│ │ scale ┼ Float64: 1.0
│ │ ub ┼ nothing
│ │ rke ┴ Bool: true
│ sigma ┼ 20×20 Matrix{Float64}
│ chol ┼ nothing
│ rc ┼ nothing
│ alg ┼ SquaredSOCRiskExpr()
│ mtx_sqrt ┼ EigenFallbackSquareRoot()
│ ratio_alg ┼ SquareRootRatio()
│ corner ┴ nothing
fb ┴ nothing
With a box mean set, the mean of each asset falls to the lower bound of its own interval. The worst-case portfolio then puts almost all its weight in the asset with the highest worst-case return. The ellipsoidal mean set bounds the means together, and its worst-case penalty grows as the weights concentrate. The worst-case portfolio spreads its weight over most of the assets.
The last cell of this section prints the largest worst-case return that a long-only portfolio can reach over each set, next to rf. When that return is below rf, the ratio is negative for every portfolio, and no tangency portfolio exists. MaximumRatio solves the ratio with a scale variable k, and here it returns the portfolio at the floor kmin of that variable, so the weights still meet their constraints. Read them as the best worst-case return at that scale, not as a worst-case tangency portfolio. The nominal problem has a tangency portfolio, so this note does not apply to the nominal column.
pretty_table(DataFrame(["Assets" => rd.nx, "Nominal" => ret_nom.w, "Box worst-case mean" => ret_box.w, "Ellipsoid worst-case mean" => ret_ell.w]); formatters = [resfmt], title = "Maximum-ratio weights with the nominal mean and the two worst-case means")plot_stacked_bar_composition([ret_nom, ret_box, ret_ell], rd; xticks = (1:3, ["Nominal", "Box μ", "Ellipsoid μ"]))
We find the long-only portfolio with the largest worst-case return over each set with MaximumReturn and NoRisk. Over a box, the worst-case return of the weights w is μᵀw - dᵀ|w|, where d is the vector of the half-widths of the intervals. Over an ellipsoid it is μᵀw - k‖Gw‖₂, where k is the radius of the set and G is the upper Cholesky factor of its covariance. We print both returns next to rf.
using LinearAlgebrawc_box(w) = dot(pr.mu, w) - dot((mu_box.ub .- mu_box.lb) ./ 2, abs.(w))wc_ell(w) = dot(pr.mu, w) - mu_ell.k * norm(cholesky(mu_ell.sigma).U * w)function best_wc(ucs) return optimise(MeanRisk(; r = NoRisk(), obj = MaximumReturn(), opt = JuMPOptimiser(; pe = pr, slv = slv, ret = ArithmeticReturn(; ucs = ucs)))).wend(rf = rf, box = wc_box(best_wc(mu_box)), ellipsoid = wc_ell(best_wc(mu_ell)))(rf = 0.0001666666666666667, box = 0.00012079562732887172, ellipsoid = -0.0006688770450720138)6. Other set estimators: Delta and the ARCH bootstrap
NormalUncertaintySet builds its sets from the sampling laws that normal returns imply, and other estimators exist. Two of them are DeltaUncertaintySet and ARCHUncertaintySet. DeltaUncertaintySet is the simpler, a box whose intervals are a fixed fraction of the point estimate on each side, with no sampling. It gives the same set on every run and costs almost nothing to build. You choose the fraction, and the data does not change it. We build a covariance box with it and print its total width next to that of the Normal box of section 2.
ucs_delta = sigma_ucs(DeltaUncertaintySet(), rd.X)set_width(u) = sum(abs, u.ub .- u.lb)pretty_table(DataFrame(; estimator = ["Delta (fixed)", "Normal (q=0.05)"], Symbol("box total width") => [set_width(ucs_delta), set_width(ucs_box)]); title = "Total width of the Delta box and the Normal box")Total width of the Delta box and the Normal box
┌─────────────────┬─────────────────┐
│ estimator │ box total width │
│ String │ Float64 │
├─────────────────┼─────────────────┤
│ Delta (fixed) │ 0.0131652 │
│ Normal (q=0.05) │ 0.0440054 │
└─────────────────┴─────────────────┘On this data the Delta box is narrower than the Normal box. UncertaintySetVariance takes it as it takes the Normal set, so we solve a robust minimum-variance portfolio with it and print its weights next to those of the Normal box of section 4.
res_delta = optimise(MeanRisk(; r = UncertaintySetVariance(; ucs = ucs_delta), obj = MinimumRisk(), opt = JuMPOptimiser(; pe = pr, slv = slv)))pretty_table(DataFrame(["Assets" => rd.nx, "Delta-robust" => res_delta.w, "Normal-box-robust" => res_box.w]); formatters = [resfmt], title = "Robust minimum-variance weights with the Delta box and the Normal box")plot_stacked_bar_composition([res_nom, res_delta, res_box], rd; xticks = (1:3, ["Nominal", "Delta", "Normal box"]))
ARCHUncertaintySet builds a set from the tails and the serial dependence of the returns themselves, not from a Gaussian or a fixed fraction. It resamples the returns in blocks, with StationaryBootstrap, MovingBootstrap or CircularBootstrap, and fits the moments again on each of the n_sim resamples. We build a covariance box with it and print its total width next to that of the Normal box.
ucs_arch = sigma_ucs(ARCHUncertaintySet(; alg = BoxUncertaintySetAlgorithm(), bootstrap = StationaryBootstrap(), n_sim = 100, seed = 1), rd.X)(arch = set_width(ucs_arch), normal = set_width(ucs_box))(arch = 0.04638754318257337, normal = 0.04400540957739831)On this data the ARCH box is a little wider than the Normal box. The block bootstrap draws from the observed returns, so their fat tails and autocorrelation enter the set, and the Normal estimator ignores both. You pass this set to UncertaintySetVariance in the same way as the other two sets.
This page was generated using Literate.jl.