The source files can be found in examples/.
Pareto surface
This example kicks up the complexity a couple of notches. We will introduce a new optimisation estimator, NearOptimalCentering optimiser.
Reach for a Pareto surface when a single risk/return frontier is not enough because you are trading off more than two competing criteria at once (e.g. two different risk measures alongside return). Where the efficient frontier sweeps a curve, this sweeps a surface/volume/hypervolume, exposing the full set of non-dominated portfolios. It builds directly on the efficient-frontier example — reach for that first if two criteria suffice.
using PortfolioOptimisers, PrettyTables# Format for pretty tables.tsfmt = (v, i, j) -> begin if j == 1 return Date(v) else return 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. ReturnsResult data
We will use the same data as the previous example.
using CSV, TimeSeries, DataFramesX = TimeArray(CSV.File(joinpath(@__DIR__, "..", "SP500.csv.gz")); timestamp = :Date)[(end - 252):end]pretty_table(X[(end - 5):end]; formatters = [tsfmt])# Compute the returnsrd = prices_to_returns(X)ReturnsResult
nx ┼ 20-element Vector{String}
X ┼ 252×20 Matrix{Float64}
nf ┼ nothing
F ┼ nothing
nb ┼ nothing
B ┼ nothing
ts ┼ 252-element Vector{Date}
iv ┼ nothing
ivpa ┼ nothing
nz ┼ nothing
Z ┴ nothing
2. Preparing solvers for pareto surface
The pareto surface is a generalisation of the efficient frontier, in fact, we can even think of hypersurfaces if we provide more parameters, but that would be difficult to visualise, so we will stick to a 2D surface in 3D space.
We'll provide a vector of solvers because the optimisation type we'll be using is more complex, and will contain various constraints.
using Clarabelslv = [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)), Solver(; name = :clarabel4, solver = Clarabel.Optimizer, settings = Dict("verbose" => false, "max_step_fraction" => 0.85), check_sol = (; allow_local = true, allow_almost = true)), Solver(; name = :clarabel5, solver = Clarabel.Optimizer, settings = Dict("verbose" => false, "max_step_fraction" => 0.8), check_sol = (; allow_local = true, allow_almost = true)), Solver(; name = :clarabel6, solver = Clarabel.Optimizer, settings = Dict("verbose" => false, "max_step_fraction" => 0.75), check_sol = (; allow_local = true, allow_almost = true)), Solver(; name = :clarabel7, solver = Clarabel.Optimizer, settings = Dict("verbose" => false, "max_step_fraction" => 0.70), check_sol = (; allow_local = true, allow_almost = true))];3. High order prior statistics
We will once again precompute the prior statistics because otherwise they'd have to be recomputed a few times.
We will be using high order risk measures, so we need to compute high order moments, we can do this with a HighOrderPriorEstimator, which needs a prior estimator that computes low order moments. Since we are only using a year of data, we will denoise our positive definite matrices by eliminating the eigenvalues corresponding to random noise. Denoising the non-positive definite matrix for the data we're using creates a negative square root, so we will not denoise it.
Note how many options this estimator contains.
dn = Denoise(; alg = SpectralDenoise(;))mp = MatrixProcessing(; dn = dn)pe = HighOrderPriorEstimator(; # Prior estimator for low order moments pe = EmpiricalPrior(; ce = PortfolioOptimisersCovariance(; mp = mp)), # Estimator for cokurtosis kte = Cokurtosis(; mp = mp), # Estimator for coskewness ske = Coskewness())HighOrderPriorEstimator
pe ┼ EmpiricalPrior
│ ce ┼ PortfolioOptimisersCovariance
│ │ ce ┼ Covariance
│ │ │ me ┼ SimpleExpectedReturns
│ │ │ │ w ┴ nothing
│ │ │ ce ┼ GeneralCovariance
│ │ │ │ ce ┼ SimpleCovariance: SimpleCovariance(true)
│ │ │ │ w ┴ nothing
│ │ │ alg ┴ FullMoment()
│ │ mp ┼ MatrixProcessing
│ │ │ pdm ┼ Posdef
│ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ dn ┼ Denoise
│ │ │ │ pdm ┼ Posdef
│ │ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ │ alg ┼ SpectralDenoise()
│ │ │ │ args ┼ Tuple{}: ()
│ │ │ │ kwargs ┼ @NamedTuple{}: NamedTuple()
│ │ │ │ kernel ┼ typeof(AverageShiftedHistograms.Kernels.gaussian): AverageShiftedHistograms.Kernels.gaussian
│ │ │ │ m ┼ Int64: 10
│ │ │ │ n ┴ Int64: 1000
│ │ │ dt ┼ nothing
│ │ │ alg ┼ nothing
│ │ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ me ┼ SimpleExpectedReturns
│ │ w ┴ nothing
│ horizon ┴ nothing
kte ┼ Cokurtosis
│ me ┼ SimpleExpectedReturns
│ │ w ┴ nothing
│ mp ┼ MatrixProcessing
│ │ pdm ┼ Posdef
│ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ dn ┼ Denoise
│ │ │ pdm ┼ Posdef
│ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ alg ┼ SpectralDenoise()
│ │ │ args ┼ Tuple{}: ()
│ │ │ kwargs ┼ @NamedTuple{}: NamedTuple()
│ │ │ kernel ┼ typeof(AverageShiftedHistograms.Kernels.gaussian): AverageShiftedHistograms.Kernels.gaussian
│ │ │ m ┼ Int64: 10
│ │ │ n ┴ Int64: 1000
│ │ dt ┼ nothing
│ │ alg ┼ nothing
│ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ alg ┼ FullMoment()
│ w ┴ nothing
ske ┼ Coskewness
│ me ┼ SimpleExpectedReturns
│ │ 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)
│ alg ┼ FullMoment()
│ w ┴ nothing
Let's compute the prior statistics.
pr = prior(pe, rd)HighOrderPrior
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
│ Z ┴ nothing
kt ┼ 400×400 Matrix{Float64}
D2 ┼ 400×210 SparseArrays.SparseMatrixCSC{Int64, Int64}
L2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
S2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
sk ┼ 20×400 Matrix{Float64}
V ┼ 20×20 Matrix{Float64}
skmp ┼ MatrixProcessing
│ pdm ┼ Posdef
│ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ dn ┼ nothing
│ dt ┼ nothing
│ alg ┼ nothing
│ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
fpr ┴ nothing
In order to generate a pareto surface/hyper-surface, we need more dimensions than we've previously explored. We can do this by adding more risk measure sweeps (and taking their product) to generate a mesh. PortfolioOptimisers does this internally and generally, but we will limit ourselves to two risk measures. This will generate a 2D surface which we can visualise in 3D.
We will use the square root NegativeSkewness and Kurtosis.
r1 = NegativeSkewness()r2 = Kurtosis()Kurtosis
settings ┼ RiskMeasureSettings
│ scale ┼ Float64: 1.0
│ ub ┼ nothing
│ rke ┴ Bool: true
w ┼ nothing
mu ┼ nothing
kt ┼ nothing
N ┼ nothing
alg1 ┼ FullMoment()
alg2 ┼ SOCRiskExpr()
pe ┴ nothing
4. Near optimal centering pareto surface
First we need to get the bounds of our pareto surface. We can do this in many different ways, the simplest are:
- Minimise the risk using both risk measures simultaneously subject to optional constraints.
- Maximise the return, utility or ratio subject to optional constraints.
We will simply maximise the risk-return ratio for both risk measures on their own with no added constraints. This will not give a complete surface, but it will give us a reasonable range of values.
The NearOptimalCentering estimator will not return the portfolio which satisfies the traditional MeanRisk constraints, but rather a portfolio which is at the centre of an analytical region (neighbourhood) around the optimal solution. The region is parametrised by binning the efficient frontier, we will use the automatic bins here, but it is possible to define them manually.
# Risk-free rate of 4.2/100/252rf = 4.2 / 100 / 252opt = JuMPOptimiser(; pe = pr, slv = slv)obj = MaximumRatio(; rf = rf)opt1 = NearOptimalCentering(; r = r1, obj = obj, opt = opt)opt2 = NearOptimalCentering(; r = r2, obj = obj, opt = opt)NearOptimalCentering
opt ┼ JuMPOptimiser
│ pe ┼ HighOrderPrior
│ │ 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
│ │ │ Z ┴ nothing
│ │ kt ┼ 400×400 Matrix{Float64}
│ │ D2 ┼ 400×210 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ L2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ S2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ sk ┼ 20×400 Matrix{Float64}
│ │ V ┼ 20×20 Matrix{Float64}
│ │ skmp ┼ MatrixProcessing
│ │ │ pdm ┼ Posdef
│ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ dn ┼ nothing
│ │ │ dt ┼ nothing
│ │ │ alg ┼ nothing
│ │ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ │ fpr ┴ nothing
│ slv ┼ 7-element Vector{Solver}
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ wb ┼ WeightBounds
│ │ lb ┼ Float64: 0.0
│ │ ub ┴ Float64: 1.0
│ bgt ┼ Float64: 1.0
│ sbgt ┼ nothing
│ gbgt ┼ nothing
│ xbgt ┼ Bool: false
│ 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
│ 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
│ linf ┼ nothing
│ lp ┼ nothing
│ brt ┼ Bool: false
│ x_src ┼ Symbol: :prior
│ z_src ┼ Symbol: :data
│ strict ┴ Bool: false
r ┼ Kurtosis
│ settings ┼ RiskMeasureSettings
│ │ scale ┼ Float64: 1.0
│ │ ub ┼ nothing
│ │ rke ┴ Bool: true
│ w ┼ nothing
│ mu ┼ nothing
│ kt ┼ nothing
│ N ┼ nothing
│ alg1 ┼ FullMoment()
│ alg2 ┼ SOCRiskExpr()
│ pe ┴ nothing
obj ┼ MaximumRatio
│ rf ┼ Float64: 0.0001666666666666667
│ ohf ┴ nothing
bins ┼ nothing
w_min ┼ nothing
w_min_ini ┼ nothing
w_opt ┼ nothing
w_opt_ini ┼ nothing
w_max ┼ nothing
w_max_ini ┼ nothing
ucs_flag ┼ Bool: true
alg ┼ UnconstrainedNearOptimalCentering()
fb ┴ nothing
Note the number of options in the estimator. In particular the alg property. Which in this case means the NearOptimalCentering alg will not have any external constraints applied to it.
Let's optimise the portfolios.
res1 = optimise(opt1)res2 = optimise(opt2)NearOptimalCenteringResult
jr ┼ JuMPOptimisationResult
│ pa ┼ ProcessedJuMPOptimiserAttributes
│ │ pr ┼ HighOrderPrior
│ │ │ 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
│ │ │ │ Z ┴ nothing
│ │ │ kt ┼ 400×400 Matrix{Float64}
│ │ │ D2 ┼ 400×210 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ │ L2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ │ S2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ │ sk ┼ 20×400 Matrix{Float64}
│ │ │ V ┼ 20×20 Matrix{Float64}
│ │ │ skmp ┼ MatrixProcessing
│ │ │ │ pdm ┼ Posdef
│ │ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ │ dn ┼ nothing
│ │ │ │ dt ┼ nothing
│ │ │ │ alg ┼ nothing
│ │ │ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ │ │ fpr ┴ 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}
│ │ sca ┴ SumScalariser()
│ retcode ┼ OptimisationSuccess
│ │ res ┴ nothing
│ sol ┼ JuMPOptimisationSolution
│ │ w ┴ 20-element Vector{Float64}
│ model ┼ A JuMP Model
│ │ ├ solver: Clarabel
│ │ ├ objective_sense: MIN_SENSE
│ │ │ └ objective_function_type: AffExpr
│ │ ├ num_variables: 273
│ │ ├ num_constraints: 47
│ │ │ ├ 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.ExponentialCone: 42
│ │ │ └ Vector{AffExpr} in MOI.PositiveSemidefiniteConeSquare: 1
│ │ └ Names registered in the model
│ │ └ :Gkt, :L2W, :M, :M_PSD, :W, :bgt, :ckurt_soc_1, :clog_delta_w, :clog_ret, :clog_risk, :clog_w, :k, :kurtosis_risk_1, :log_delta_w, :log_ret, :log_risk, :log_w, :lw, :noc_rk, :noc_rt, :obj_expr, :ret, :ret_1, :ret_vec, :risk, :risk_vec, :sc, :so, :w, :w_lb, :w_ub, :x_kurt_1
r ┼ Kurtosis
│ settings ┼ RiskMeasureSettings
│ │ scale ┼ Float64: 1.0
│ │ ub ┼ nothing
│ │ rke ┴ Bool: true
│ w ┼ nothing
│ mu ┼ 20-element Vector{Float64}
│ kt ┼ 400×400 Matrix{Float64}
│ N ┼ nothing
│ alg1 ┼ FullMoment()
│ alg2 ┼ SOCRiskExpr()
│ pe ┴ nothing
w_min_retcode ┼ OptimisationSuccess
│ res ┴ Dict{Any, Any}: Dict{Any, Any}()
w_opt_retcode ┼ OptimisationSuccess
│ res ┴ Dict{Any, Any}: Dict{Any, Any}()
w_max_retcode ┼ OptimisationSuccess
│ res ┴ Dict{Any, Any}: Dict{Any, Any}()
noc_retcode ┼ OptimisationSuccess
│ res ┴ Dict{Any, Any}: Dict{Any, Any}()
fb ┴ nothing
In order to allow for multiple risk measures in optimisations, certain measures can take different parameters. In this case, NegativeSkewness and Kurtosis take the moment matrices, which are used to compute the risk measures. We can use the factory function to create a new risk measure with the same parameters as the original, but with the moment matrices from the prior. Other risk measures require a solver, and this function is also used in those cases.
r1 = factory(r1, pr)r2 = factory(r2, pr)Kurtosis
settings ┼ RiskMeasureSettings
│ scale ┼ Float64: 1.0
│ ub ┼ nothing
│ rke ┴ Bool: true
w ┼ nothing
mu ┼ 20-element Vector{Float64}
kt ┼ 400×400 Matrix{Float64}
N ┼ nothing
alg1 ┼ FullMoment()
alg2 ┼ SOCRiskExpr()
pe ┴ nothing
Let's compute the risk bounds for the pareto surface. We need to compute four risks because we have two risk measures and two optimisations. This will let us pick the lower and upper bounds for each risk measure, as we explore the pareto surface from one optimisation to the other.
sk_rk1 = expected_risk(r1, res1.w, pr.X);kt_rk1 = expected_risk(r2, res1.w, pr.X);sk_rk2 = expected_risk(r1, res2.w, pr.X);kt_rk2 = expected_risk(r2, res2.w, pr.X);We will now create new risk measures bounded by these values. We will also use factories from the get-go. The optimisation procedure prioritises the parameters in the risk measures over the ones in the prior. This lets users provide the same risk measure with different parameters in the same optimisation. We will use two ranges of 5. The total number of points in the pareto surface will be the product of the points of each range.
Since we don't know which sk_rk1 or sk_r2, kt_rk1 or kt_rk2 is bigger or smaller, we need to use min, max.
r1 = factory(NegativeSkewness(; settings = RiskMeasureSettings(; # Risk upper bounds go from the minimum to maximum risk given the optimisations. ub = range(; start = min(sk_rk1, sk_rk2), stop = max(sk_rk1, sk_rk2), length = 5))), pr);r2 = factory(Kurtosis(; settings = RiskMeasureSettings(; ub = range(; start = min(kt_rk1, kt_rk2), stop = max(kt_rk1, kt_rk2), length = 5))), pr);Now we only need to maximise the return given both risk measures. Internally, the optimisation will generate the mesh as a product of the ranges in the order in which the risk measures were provided. This also works with the MeanRisk estimator, in fact, NearOptimalCentering uses it internally.
Since we are using an unconstrained NearOptimalCentering, the risk bound constraints will not be satisfied by the solution. If we wish to satisfy them, we can provide alg = ConstrainedNearOptimalCentering(), but would also make the optimisations harder, which may cause them to fail.
opt3 = NearOptimalCentering(; r = [r1, r2], obj = MaximumReturn(), opt = opt)NearOptimalCentering
opt ┼ JuMPOptimiser
│ pe ┼ HighOrderPrior
│ │ 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
│ │ │ Z ┴ nothing
│ │ kt ┼ 400×400 Matrix{Float64}
│ │ D2 ┼ 400×210 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ L2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ S2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ sk ┼ 20×400 Matrix{Float64}
│ │ V ┼ 20×20 Matrix{Float64}
│ │ skmp ┼ MatrixProcessing
│ │ │ pdm ┼ Posdef
│ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ dn ┼ nothing
│ │ │ dt ┼ nothing
│ │ │ alg ┼ nothing
│ │ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ │ fpr ┴ nothing
│ slv ┼ 7-element Vector{Solver}
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ │ Solver ⋯
│ wb ┼ WeightBounds
│ │ lb ┼ Float64: 0.0
│ │ ub ┴ Float64: 1.0
│ bgt ┼ Float64: 1.0
│ sbgt ┼ nothing
│ gbgt ┼ nothing
│ xbgt ┼ Bool: false
│ 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
│ 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
│ linf ┼ nothing
│ lp ┼ nothing
│ brt ┼ Bool: false
│ x_src ┼ Symbol: :prior
│ z_src ┼ Symbol: :data
│ strict ┴ Bool: false
r ┼ 2-element Vector{RiskMeasure}
│ NegativeSkewness ⋯
│ Kurtosis ⋯
obj ┼ MaximumReturn()
bins ┼ nothing
w_min ┼ nothing
w_min_ini ┼ nothing
w_opt ┼ nothing
w_opt_ini ┼ nothing
w_max ┼ nothing
w_max_ini ┼ nothing
ucs_flag ┼ Bool: true
alg ┼ UnconstrainedNearOptimalCentering()
fb ┴ nothing
See how r is a vector of risk measures with populated properties. We can now optimise the portfolios.
res3 = optimise(opt3)NearOptimalCenteringResult
jr ┼ JuMPOptimisationResult
│ pa ┼ ProcessedJuMPOptimiserAttributes
│ │ pr ┼ HighOrderPrior
│ │ │ 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
│ │ │ │ Z ┴ nothing
│ │ │ kt ┼ 400×400 Matrix{Float64}
│ │ │ D2 ┼ 400×210 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ │ L2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ │ S2 ┼ 210×400 SparseArrays.SparseMatrixCSC{Int64, Int64}
│ │ │ sk ┼ 20×400 Matrix{Float64}
│ │ │ V ┼ 20×20 Matrix{Float64}
│ │ │ skmp ┼ MatrixProcessing
│ │ │ │ pdm ┼ Posdef
│ │ │ │ │ alg ┼ UnionAll: NearestCorrelationMatrix.Newton
│ │ │ │ │ kwargs ┴ @NamedTuple{}: NamedTuple()
│ │ │ │ dn ┼ nothing
│ │ │ │ dt ┼ nothing
│ │ │ │ alg ┼ nothing
│ │ │ │ order ┴ NTuple{4, Symbol}: (:pdm, :dn, :dt, :alg)
│ │ │ fpr ┴ 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}
│ │ sca ┴ SumScalariser()
│ retcode ┼ OptimisationSuccess
│ │ res ┴ nothing
│ sol ┼ 25-element Vector{JuMPOptimisationSolution}
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ │ JuMPOptimisationSolution ⋯
│ model ┼ A JuMP Model
│ │ ├ solver: Clarabel
│ │ ├ objective_sense: MIN_SENSE
│ │ │ └ objective_function_type: AffExpr
│ │ ├ num_variables: 276
│ │ ├ num_constraints: 50
│ │ │ ├ AffExpr in MOI.EqualTo{Float64}: 1
│ │ │ ├ Vector{AffExpr} in MOI.Nonnegatives: 1
│ │ │ ├ Vector{AffExpr} in MOI.Nonpositives: 1
│ │ │ ├ Vector{AffExpr} in MOI.SecondOrderCone: 2
│ │ │ ├ Vector{AffExpr} in MOI.ExponentialCone: 42
│ │ │ ├ Vector{AffExpr} in MOI.PositiveSemidefiniteConeSquare: 1
│ │ │ └ VariableRef in MOI.Parameter{Float64}: 2
│ │ └ Names registered in the model
│ │ └ :L2W, :M, :M_PSD, :W, :bgt, :ckurt_soc_2, :clog_delta_w, :clog_ret, :clog_risk, :clog_w, :cnskew_soc_1, :k, :kurtosis_risk_2, :log_delta_w, :log_ret, :log_risk, :log_w, :lw, :noc_rk, :noc_rt, :nskew_risk_1, :obj_expr, :ret, :ret_1, :ret_vec, :risk, :risk_vec, :sc, :so, :w, :w_lb, :w_ub, :x_kurt_2
r ┼ 2-element Vector{RiskMeasure}
│ NegativeSkewness ⋯
│ Kurtosis ⋯
w_min_retcode ┼ OptimisationSuccess
│ res ┴ Dict{Any, Any}: Dict{Any, Any}()
w_opt_retcode ┼ 25-element Vector{OptimisationSuccess}
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
w_max_retcode ┼ OptimisationSuccess
│ res ┴ Dict{Any, Any}: Dict{Any, Any}()
noc_retcode ┼ 25-element Vector{OptimisationSuccess}
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
│ OptimisationSuccess ⋯
fb ┴ nothing
As expected, there are 5 × 5 = 25 solutions. Thankfully there are no warnings about failed optimisations, so there is no need to check the solutions.
The NearOptimalCentering estimator contains various return codes because it may need to compute some MeanRisk optimisations, it has a retcode which summarises whether all other optimisations succeeded. We can check this to make sure it was a success.
isa(res3.retcode, OptimisationSuccess)true5. Visualising the pareto surface
Let's view how the weights evolve along the pareto surface.
using StatsPlots, GraphRecipesplot_stacked_area_composition(res3.w, rd.nx)
Now we can view the pareto surface. For the z-axis and colourbar, we will use the conditional drawdown at risk to return ratio.
plot_measures(res3.w, pr; x = r1, y = r2, z = ExpectedReturnRiskRatio(; rk = ConditionalDrawdownatRisk(), rt = ArithmeticReturn(), rf = rf), c = ExpectedReturnRiskRatio(; rk = ConditionalDrawdownatRisk(), rt = ArithmeticReturn(), rf = rf), title = "Pareto Surface", xlabel = "Sqrt NSkew", ylabel = "Sqrt Kurt", zlabel = "CDaR/Return")
We can view it in 2D as well.
plot_measures(res3.w, pr; x = r1, y = r2, c = ExpectedReturnRiskRatio(; rk = ConditionalDrawdownatRisk(), rt = ArithmeticReturn(), rf = rf), title = "Pareto Front", xlabel = "Sqrt NSkew", ylabel = "Sqrt Kurt", colorbar_title = "\n\nCDaR/Return", right_margin = 8Plots.mm)
This page was generated using Literate.jl.