The source files can be found in examples/.

Asset pre-selection

Some assets in a universe add nothing to the optimisation. A constant column, such as a name that never traded, carries no information. Some assets rank too low, for example when you want the twenty names with the lowest risk and the universe has two hundred. Some are redundant, because they move so closely with another asset that keeping both gives the optimiser the same information twice. The optimiser then treats one common factor as two separate sources of risk.

An asset selector removes such assets from the universe, based on the returns. It is an ordinary preprocessing estimator and needs no pipeline. A Pipeline runs it through fit_preprocessing and apply_preprocessing, as it runs a prior estimator through prior. The library has three kinds:

• CompleteAssetSelector drops every asset with a non-finite return at any observation of the window.
• ScoreSelector scores every asset with a risk measure and keeps the assets a rule admits. ZeroVarianceFilter is a named case of it.
• RedundancySelector drops assets whose returns repeat the information of other assets.

A selector chooses its universe on the training window, and that universe is the result of the fit. Applying the fitted result to a later window keeps the same universe and does not choose again. A selector that chose again on each window would pick its assets with the returns it is about to be scored on, which is look-ahead bias.

using PortfolioOptimisers, PrettyTables, DataFrames, Statistics, StatsAPIresfmt = (v, i, j) -> begin    if isa(v, Number) && !isa(v, Integer)        return "$(round(v * 100, digits = 3)) %"    else        return v    endend;

1. The data

We take twenty S&P 500 names over the last 1000 trading days. Section 5 needs the correlations of real data, because on them the two redundancy algorithms keep different assets.

using CSV, TimeSeriesX = TimeArray(CSV.File(joinpath(@__DIR__, "..", "SP500.csv.gz")); timestamp = :Date)rd = prices_to_returns(X[(end - 1000):end])size(rd.X), rd.nx

((1000, 20), ["AAPL", "AMD", "BAC", "BBY", "CVX", "GE", "HD", "JNJ", "JPM", "KO", "LLY", "MRK", "MSFT", "PEP", "PFE", "PG", "RRC", "UNH", "WMT", "XOM"])

2. Score an asset with a risk measure

ScoreSelector scores asset i with a risk measure evaluated on the returns of that asset alone, score(view(X, :, i)). The library has no separate family of scores, because the risk measures already cover what you would rank assets on: variance and semi-variance, VaR and CVaR, the drawdown measures, and MeanReturn.

Two traits decide how a measure scores. supports_precomputed_returns says whether a measure can take a single return series, and bigger_is_better says which end of the ranking is the best. We print both traits for five measures.

measures = ["SCM()" => SCM(), "ConditionalValueatRisk()" => ConditionalValueatRisk(),            "MaximumDrawdown()" => MaximumDrawdown(), "MeanReturn()" => MeanReturn(),            "Variance()" => Variance()]pretty_table(DataFrame(; measure = first.(measures),                       scoreable = [PortfolioOptimisers.supports_precomputed_returns(r)                                    for (_, r) in measures],                       bigger_is_better = [PortfolioOptimisers.bigger_is_better(r)                                           for (_, r) in measures]))

┌──────────────────────────┬───────────┬──────────────────┐
│                  measure │ scoreable │ bigger_is_better │
│                   String │      Bool │             Bool │
├──────────────────────────┼───────────┼──────────────────┤
│                    SCM() │      true │            false │
│ ConditionalValueatRisk() │      true │            false │
│        MaximumDrawdown() │      true │            false │
│             MeanReturn() │      true │             true │
│               Variance() │     false │            false │
└──────────────────────────┴───────────┴──────────────────┘

A smaller ConditionalValueatRisk is better and a bigger MeanReturn is better. A rule that asks for the best assets takes the right end for each measure, and you do not state which end it is. We score every asset with both measures.

cvar = PortfolioOptimisers.asset_scores(ConditionalValueatRisk(), rd.X)mu = PortfolioOptimisers.asset_scores(MeanReturn(), rd.X)pretty_table(sort(DataFrame(; asset = rd.nx, cvar = cvar, mean_return = mu), :cvar);             formatters = [resfmt])

┌────────┬─────────┬─────────────┐
│  asset │    cvar │ mean_return │
│ String │ Float64 │     Float64 │
├────────┼─────────┼─────────────┤
│    JNJ │ 3.045 % │      0.05 % │
│    WMT │  3.25 % │     0.058 % │
│    PEP │ 3.276 % │     0.074 % │
│     PG │  3.38 % │     0.072 % │
│    MRK │ 3.562 % │     0.067 % │
│     KO │ 3.615 % │     0.054 % │
│    PFE │ 3.807 % │     0.051 % │
│    LLY │ 3.881 % │      0.14 % │
│    UNH │ 4.288 % │     0.103 % │
│   MSFT │ 4.509 % │     0.105 % │
│     HD │  4.51 % │     0.084 % │
│    JPM │  4.67 % │     0.062 % │
│   AAPL │ 4.785 % │     0.146 % │
│    XOM │ 5.122 % │     0.089 % │
│    BAC │ 5.236 % │     0.061 % │
│    CVX │ 5.273 % │     0.092 % │
│    BBY │ 6.213 % │      0.08 % │
│     GE │ 6.273 % │     0.063 % │
│    AMD │ 7.223 % │     0.172 % │
│    RRC │ 9.046 % │     0.187 % │
└────────┴─────────┴─────────────┘

2.1 Variance cannot score a single asset

One measure of the family does not work as a score. Variance is a WeightsInput measure. It takes the weights of a portfolio and a covariance matrix, not a return series. It has no meaning for one asset alone, so ScoreSelector throws an error when you construct it with Variance.

try    ScoreSelector(; score = Variance(), rule = ThresholdRule(; lo = 0.0))catch e    println(e.msg)end

`Variance` cannot score a single asset's return series: its `supports_precomputed_returns` is false, so its functor consumes portfolio weights rather than a precomputed return vector. `SCM()` computes the same quantity from a return series and is scoreable.

The error names the replacement, SCM(). It is the second central moment, an alias for LowOrderMoment with SecondMoment and FullMoment. It computes the same quantity from a return series, and it can score an asset. We print its scores for three assets beside the variances that var computes.

PortfolioOptimisers.asset_scores(SCM(), rd.X)[1:3],var(rd.X[:, 1:3]; dims = 1, corrected = false)

([0.0004644201994021774, 0.001189768278363829, 0.000565865613981831], [0.0004639557792027752 0.0011885785100854654 0.0005652997483678491])

3. Rules: what to do with the scores

A selection rule turns the scores into a mask of the assets to keep. There are two kinds of rule, and they treat the scores differently.

A literal rule compares the scores with fixed values. ThresholdRule compares each score with the bounds (lo, hi), both optional and both exclusive, and ignores bigger_is_better. A filter for zero variance must drop the assets of low variance, and a rule that kept the better assets by the trait would keep those assets instead.

An ordinal rule compares the scores with each other. RankRule and QuantileRule sort the assets from best to worst by bigger_is_better, and take a count or a fraction of assets from each end. A rule takes counts, not positions, and it can take the two ends at once. action = :drop inverts the selection, so the rule drops the assets it names and keeps the rest. With it you can drop the worst five without knowing how many assets there are.

rules = ["RankRule(best = 5)" => RankRule(; best = 5),         "RankRule(worst = 5)" => RankRule(; worst = 5),         "RankRule(worst = 3, action = :drop)" => RankRule(; worst = 3, action = :drop),         "RankRule(best = 2, worst = 2)" => RankRule(; best = 2, worst = 2),         "QuantileRule(best = 0.25)" => QuantileRule(; best = 0.25)]rows = map(rules) do (label, rule)    kept = fit_preprocessing(ScoreSelector(; score = ConditionalValueatRisk(), rule = rule),                             rd).nx    return (; rule = label, n = length(kept), kept = join(kept, ", "))endpretty_table(DataFrame(rows))

┌─────────────────────────────────────┬───────┬────────────────────────────────────────────────────────────────────────────────────┐
│                                rule │     n │                                                                               kept │
│                              String │ Int64 │                                                                             String │
├─────────────────────────────────────┼───────┼────────────────────────────────────────────────────────────────────────────────────┤
│                  RankRule(best = 5) │     5 │                                                             JNJ, MRK, PEP, PG, WMT │
│                 RankRule(worst = 5) │     5 │                                                             AMD, BBY, CVX, GE, RRC │
│ RankRule(worst = 3, action = :drop) │    17 │ AAPL, BAC, BBY, CVX, HD, JNJ, JPM, KO, LLY, MRK, MSFT, PEP, PFE, PG, UNH, WMT, XOM │
│       RankRule(best = 2, worst = 2) │     4 │                                                                 AMD, JNJ, RRC, WMT │
│           QuantileRule(best = 0.25) │     5 │                                                             JNJ, MRK, PEP, PG, WMT │
└─────────────────────────────────────┴───────┴────────────────────────────────────────────────────────────────────────────────────┘

The direction of the ranking changes with the measure. ConditionalValueatRisk is better when smaller, so best = 5 returns the five defensive names. MeanReturn is better when bigger, so the same rule returns the five growth names.

DataFrame(; rule = ["best = 5", "worst = 5"],          by_cvar = [join(fit_preprocessing(ScoreSelector(;                                                          score = ConditionalValueatRisk(),                                                          rule = RankRule(; best = 5)), rd).nx,                          ", "),                     join(fit_preprocessing(ScoreSelector(;                                                          score = ConditionalValueatRisk(),                                                          rule = RankRule(; worst = 5)),                                            rd).nx, ", ")],          by_mean_return = [join(fit_preprocessing(ScoreSelector(; score = MeanReturn(),                                                                 rule = RankRule(;                                                                                 best = 5)),                                                   rd).nx, ", "),                            join(fit_preprocessing(ScoreSelector(; score = MeanReturn(),                                                                 rule = RankRule(;                                                                                 worst = 5)),                                                   rd).nx, ", ")])

3.1 Degenerate columns

ZeroVarianceFilter is ScoreSelector(; score = SCM(), rule = ThresholdRule(; lo = tol)) under its own name. A constant column has zero variance, adds nothing to a portfolio, and makes a covariance matrix singular. We set the returns of one asset to zero, and the filter drops it.

CompleteAssetSelector drops the assets with a non-finite return. Use it when a pipeline receives a ReturnsResult directly, so the MissingDataFilter that works on prices never runs.

Xz = copy(rd.X)Xz[:, 4] .= 0.0rd_z = ReturnsResult(; nx = rd.nx, X = Xz)setdiff(rd.nx, fit_preprocessing(ZeroVarianceFilter(), rd_z).nx)

1-element Vector{String}:
 "BBY"

tol defaults to 1e-12 and not 0, because a column that moves by 1e-18 is constant for any use. The bound is exclusive, so tol = 0 still drops a column that is exactly constant.

4. Ties: if the scores cannot tell two assets apart, the rule keeps neither

When two assets have the same score, a rule keeps neither of them. A block of tied assets that spans a cut of the ranking is dropped whole and never split. To split it, the rule would have to break the tie by column index, and your portfolio would then depend on the order of the columns in your CSV file, which is not a property of the data.

So RankRule(; best = k) can return fewer than k assets. To show it, we copy the column of the asset with the second lowest variance over the column of the third, so that the two tie, and ask for the best one, two and three assets.

v = PortfolioOptimisers.asset_scores(SCM(), rd.X)ord = sortperm(v)Xs = copy(rd.X)Xs[:, ord[3]] = Xs[:, ord[2]]rd_tie = ReturnsResult(; nx = rd.nx, X = Xs)tie_rows = map([1, 2, 3]) do k    kept = fit_preprocessing(ScoreSelector(; score = SCM(), rule = RankRule(; best = k)),                             rd_tie).nx    return (; requested = k, returned = length(kept), kept = join(kept, ", "))endpretty_table(DataFrame(tie_rows))

┌───────────┬──────────┬─────────────┐
│ requested │ returned │        kept │
│     Int64 │    Int64 │      String │
├───────────┼──────────┼─────────────┤
│         1 │        1 │         JNJ │
│         2 │        1 │         JNJ │
│         3 │        3 │ JNJ, KO, PG │
└───────────┴──────────┴─────────────┘

best = 2 returns one asset, because ranks 2 and 3 tie and the cut falls between them, so the rule drops the pair. best = 3 returns three, because the pair now sits inside the cut and the rule keeps it whole. A window whose scores are all equal selects nothing, and fit_preprocessing throws an error rather than return an empty universe.

The redundancy selector treats ties in the same way. Two identical columns have a correlation of one and the same score, so it keeps neither. We add a copy of AAPL and check whether either stays.

rd_dup = ReturnsResult(; nx = [rd.nx; "AAPL_copy"], X = hcat(rd.X, rd.X[:, 1]))dup_kept = fit_preprocessing(RedundancySelector(; alg = PairwiseCorrelation(; t = 0.99),                                                score = SCM()), rd_dup).nx("AAPL" in dup_kept, "AAPL_copy" in dup_kept, length(dup_kept))

(false, false, 19)

5. Redundancy: pairwise pruning and connected components

RedundancySelector has two parts. alg decides which assets are redundant, and score decides which member of a redundant group stays. If you leave score as nothing, the correlation algorithms keep the asset with the lowest summary correlation to the rest of the universe, the least redundant one.

PairwiseCorrelation is greedy and is the default. It visits the pairs of assets from the most correlated to the least, and drops the worse asset of each pair, until no pair that remains has a correlation above t.

CorrelationComponents follows the correlations from one asset to the next. The assets are the nodes of a graph, and each correlation above the threshold is an edge. Each connected component of the graph keeps one asset. At t = 0.7, if ρ(A,B) = 0.80 and ρ(B,C) = 0.81 but ρ(A,C) = 0.32, the three assets form one component, and the algorithm drops two of them, although A and C are weakly correlated. We call such a sequence of edges a chain.

Chains form on real data too. We run both algorithms at the same threshold and print the largest correlation between two assets the greedy algorithm keeps.

greedy = RedundancySelector(; alg = PairwiseCorrelation(; t = 0.65, absolute = true),                            score = SCM())comps = RedundancySelector(; alg = CorrelationComponents(; t = 0.65, absolute = true),                           score = SCM())kept_g = fit_preprocessing(greedy, rd).nxkept_c = fit_preprocessing(comps, rd).nxgi = [findfirst(==(n), rd.nx) for n in kept_g]sub = abs.(cor(rd.X[:, gi]))max_surviving = maximum(sub[i, j] for j in axes(sub, 2) for i in (j + 1):size(sub, 1))DataFrame(; algorithm = ["PairwiseCorrelation", "CorrelationComponents"],          kept = [length(kept_g), length(kept_c)],          max_surviving_abs_cor = [round(max_surviving; digits = 3), missing],          extra_drops = ["", join(setdiff(kept_g, kept_c), ", ")])

absolute = true counts a correlation of -0.9 as redundant too. That is usually what you want, because two assets that move together carry the same information whatever the sign of their correlation.

The largest correlation the greedy algorithm keeps is below the threshold, and the last column names the assets that the components algorithm drops in addition. Choose the one whose rule you want.

5.1 Clustering as the grouping rule

ClusterGroups forms the groups with clusterise, so any clustering method of the library can define which assets are redundant: hierarchical linkage, the Direct Bubble Hierarchical Tree (DBHT), and the estimators of the optimal number of clusters. It has no default rule for which asset stays, so it needs a score, and the second cell prints the error you get without one.

clustered = RedundancySelector(; alg = ClusterGroups(), score = SCM())kept_cl = fit_preprocessing(clustered, rd).nxlength(kept_cl), join(kept_cl, ", ")try    RedundancySelector(; alg = ClusterGroups())catch e    println(e.msg)end

a ClusterGroups redundancy algorithm cannot choose the survivor of a group on its own; give the RedundancySelector a score

6. The universe is the result of the fit

Every selection above was fitted on one window, the whole sample. In a pipeline that window is the training window, and the fitted result holds the selected universe. A prediction on a later window keeps that universe. We fit a selector on the first 700 observations, apply it to the rest, and compare the result with the universe the selector would choose on the rest alone.

selector = ScoreSelector(; score = ConditionalValueatRisk(), rule = RankRule(; best = 10))train = ReturnsResult(; nx = rd.nx, X = rd.X[1:700, :])test = ReturnsResult(; nx = rd.nx, X = rd.X[701:end, :])fitted = fit_preprocessing(selector, train)replayed = apply_preprocessing(fitted, test)would_have_chosen = fit_preprocessing(selector, test).nxDataFrame(;          universe = ["fitted on the training window", "applied to the test window",                      "chosen on the test window"],          assets = [join(fitted.nx, ", "), join(replayed.nx, ", "),                    join(would_have_chosen, ", ")])

The second row, the universe applied to the test window, is the one fitted on the training window. The last row is the choice the test window would make from the returns it is scored on, and using it would be look-ahead bias.

6.1 The pipeline checks the order of its steps

A selector changes the assets of :returns. A step computed from :returns before the selector, such as a prior, a phylogeny, an uncertainty set or a set of constraints, then describes assets that are no longer in the universe. Pipeline throws an error for that order when you construct it, and the prior with the wrong number of assets never reaches the optimiser.

try    Pipeline(; steps = (EmpiricalPrior(), ZeroVarianceFilter(), EqualWeighted()))catch e    println(e.msg)end

a ScoreSelector step writes the :returns slot, invalidating the :prior slot written by an earlier EmpiricalPrior step; the stale :prior result would no longer match the assets of the new :returns data. Move the ScoreSelector step before the EmpiricalPrior step, or drop one of them.

With the selector first, the pipeline accepts the other steps with no change to them.

using Clarabelslv = Solver(; name = :clarabel, solver = Clarabel.Optimizer,             settings = Dict("verbose" => false),             check_sol = (; allow_local = true, allow_almost = true))pipe = Pipeline(;                steps = ("select" => ScoreSelector(; score = ConditionalValueatRisk(),                                                   rule = RankRule(; best = 10)),                         EmpiricalPrior(),                         "opt" => MeanRisk(; opt = JuMPOptimiser(; slv = slv))))res = StatsAPI.fit(pipe, rd)pretty_table(DataFrame(; asset = res.ctx.returns.nx, weight = res.w); formatters = [resfmt])

┌────────┬──────────┐
│  asset │   weight │
│ String │  Float64 │
├────────┼──────────┤
│    JNJ │ 25.982 % │
│     KO │ 17.616 % │
│    LLY │    0.0 % │
│    MRK │ 17.463 % │
│   MSFT │    0.0 % │
│    PEP │    0.0 % │
│    PFE │  6.215 % │
│     PG │  4.853 % │
│    UNH │    0.0 % │
│    WMT │ 27.871 % │
└────────┴──────────┘

The pipeline takes twenty assets and returns ten weights. predict on any window keeps the ten assets of the fit. We predict on the whole sample and print the CVaR of the portfolio of ten assets.

pred = StatsAPI.predict(res, rd)expected_risk(ConditionalValueatRisk(), pred)

0.024787275255167128

7. Selection is a parameter to tune

The number of assets to keep is a choice you cannot make by looking at the returns you will be scored on. search_cross_validation fits the whole pipeline, the selector with it, on each training window and scores it on the window held out, so no candidate sees the test window when it chooses its universe.

A key of the grid names a step and one of its fields. "select.rule" replaces the whole rule of the step named "select".

p = ["select.rule" => [RankRule(; best = 5), RankRule(; best = 10), RankRule(; best = 15)]]gscv = GridSearchCrossValidation(p; cv = IndexWalkForward(500, 250),                                 r = ConditionalValueatRisk())tuned = search_cross_validation(pipe, gscv, rd)pretty_table(DataFrame(; k = [v[1].best for v in tuned.val_grid],                       mean_score = vec(mean(tuned.test_scores; dims = 1))))

┌───────┬────────────┐
│     k │ mean_score │
│ Int64 │    Float64 │
├───────┼────────────┤
│     5 │ -0.0214013 │
│    10 │ -0.0204401 │
│    15 │ -0.0204521 │
└───────┴────────────┘

tuned.opt is the pipeline with the best score, selector included, and you can fit it on the full sample. Its selector keeps this many assets:

tuned.opt.steps[1].rule.best

10

A count larger than the number of assets keeps every asset, so best = 50 on a universe of 20 assets keeps 20 and throws no error, and the largest value of a grid never stops a search. Every other degenerate case throws an error: a measure that cannot score an asset, a rule that selects nothing, a score that is not finite, and a fitted asset that is missing from a test window.

8. Plot the two decisions

A score and a rule sort the assets along one axis and cut them at one place. The bar chart shows the CVaR of each asset, and the blue bars are the ten assets that RankRule(; best = 10) keeps.

using StatsPlotsperm = sortperm(cvar)kept10 = fit_preprocessing(ScoreSelector(; score = ConditionalValueatRisk(),                                         rule = RankRule(; best = 10)), rd).nxcolours = [n in kept10 ? :steelblue : :lightgray for n in rd.nx[perm]]bar(1:length(perm), cvar[perm]; color = colours, legend = false,    xticks = (1:length(perm), rd.nx[perm]), xrotation = 60, ylabel = "CVaR",    title = "CVaR per asset, with the ten that RankRule(best = 10) keeps in blue")

The redundancy decision depends on the threshold. We count the assets each algorithm keeps over a range of thresholds. The gap between the two lines is the assets the components algorithm drops through chains, and it closes at the thresholds where the graph has no chains.

thrs = 0.5:0.025:0.85n_greedy = [length(fit_preprocessing(RedundancySelector(;                                                        alg = PairwiseCorrelation(; t = t,                                                                                  absolute = true),                                                        score = SCM()), rd).nx)            for t in thrs]n_comps = [length(fit_preprocessing(RedundancySelector(;                                                       alg = CorrelationComponents(; t = t,                                                                                   absolute = true),                                                       score = SCM()), rd).nx)           for t in thrs]plot(thrs, n_greedy; label = "PairwiseCorrelation (greedy)", marker = :circle, lw = 2)plot!(thrs, n_comps; label = "CorrelationComponents (transitive)", marker = :square, lw = 2)hline!([length(rd.nx)]; label = "full universe", ls = :dash, color = :gray)plot!(; xlabel = "|correlation| threshold", ylabel = "assets kept",      title = "Assets kept at each correlation threshold", legend = :bottomright)

Summary

• An asset selector is a preprocessing estimator that works on returns. It chooses the universe on the training window, and applies that universe to later windows, so the selection uses no returns from the window it is scored on.
• ScoreSelector scores each asset with any risk measure whose supports_precomputed_returns is true. bigger_is_better sets the direction of best and worst, so you pass no flag for it. Variance cannot score an asset, and SCM() replaces it.
• ThresholdRule compares each score with fixed bounds and ignores the direction. RankRule and QuantileRule rank the assets and take counts or fractions from each end, so you can drop the worst five without knowing how many assets there are.
• A rule keeps neither of two tied assets and never splits a tied block. RankRule(; best = k) can return fewer than k assets, and two identical columns both leave the universe.
• PairwiseCorrelation keeps no pair above the threshold. CorrelationComponents follows chains of correlations and drops more assets. On the same returns the two keep different assets.
• A pipeline throws an error at construction if a step changes :returns after a prior, a phylogeny, an uncertainty set or a constraint step, so a step never computes on a universe that no longer exists.
• The thresholds and counts of a selection are parameters. Tune them with search_cross_validation, never on the full sample.

This page was generated using Literate.jl.