Explaining a trading P&L: what finance knows about performance attribution, and what is open

0. From diagnosing orders to explaining a P&L

A trading desk has lost money. It tracks spread, volatility, order flow, inventory and the signals used to trade, and wants to know which indicators explain the loss. Here, explain can mean finding an unusual market condition, assigning a share of the loss to each indicator, or checking that a relationship found in one sample still holds on new data.

The route of this page, from diagnosing orders to explaining a P&L.

Charles-Albert Lehalle’s talk describes an online diagnostic from Azencott et al. (2014), a paper on real-time transaction cost analysis. Every five minutes it labels as bad the orders with slippage among the worst 3% or 5% of active orders. It searches for threshold rules on market variables, called descriptors, and on pairs of variables that separate bad orders from the rest. These include volume, volatility and spread measures, rarity scores and anomaly detectors. A rule must reach a floor on both its hit rate on bad orders (sensitivity) and its true-negative rate on the rest (specificity). The paper uses 85%. A group’s influence score is the best hit rate found subject to that floor, or zero if no rule qualifies.

To explain portfolio P&L, we also need rules for sharing contributions between correlated indicators and testing relationships found by searching the data.

From orders to P&L. The diagnostic identifies orders that belong to a poorly performing group at a particular time. A P&L attribution starts with a continuous amount, positive or negative, accumulated over an interval. Finding where losses occurred does not tell us how much each indicator contributed.

The contributions must reconcile. Several descriptors can each give a useful alarm, but their scores are classification rates, not amounts of money, so they do not add up to the loss. To make the contributions add up, we must specify the total, use compatible units and show a residual for whatever remains unallocated.

Correlated indicators share information. In Azencott and colleagues’ example, volume, volatility and spread scores identify the same episode. Crediting each with the whole episode would count it repeatedly. We need a rule for sharing that information; choosing only the strongest score leaves the other contributions unresolved.

The search affects the statistical test. Azencott and colleagues test 406 candidate explanations (single variables and pairs) in each of 79 five-minute slices, with thresholds selected in sample (pp. 26–27). They restrict model complexity, but do not control false discoveries over that search. Thresholding a diagnostic score does not control errors across multiple tests. If we screen many indicators against past losses, we need to control errors across that search too.

We follow a roughly historical thread, starting with portfolio attribution against a benchmark in the 1980s, then moving to factor returns and risk, and to the signals that drive a portfolio’s positions. From there we look at correlated indicators and the market maker’s ledger, before examining the statistical pitfalls and possible research projects. The order follows the problems more closely than the dates: some of the tools come from older work, and we return to them when we need them.

↑ Sections

1. What it means to attribute a P&L

When we compare a portfolio with an index, we attribute its active return, $R_P-R_B$ , where $R_P$ is the portfolio return and $R_B$ the index return, measured over the same interval for both. For a trading book, we might compare P&L with a paper portfolio: the trades intended at decision time, valued without execution costs (section 5.1). Or we can fix a valuation rule, such as marking every position to the mid, and explain that P&L. We need to state this reference to know what the contributions mean.

Write the accounting requirement as

\Pi-\Pi^{\mathrm{ref}}=\sum_j C_j+\varepsilon.

Here $\Pi$ is the realised P&L, $\Pi^{\mathrm{ref}}$ is the reference P&L, $C_j$ is the amount assigned to component $j$ , and $\varepsilon$ is whatever the chosen model leaves unallocated. Keep this residual visible. Setting it equal to the unallocated difference makes the amounts add up; it does not show that the other terms explain behaviour, predict returns or identify causes.

1.1 What daily snapshots miss

In his book The Elements of Quantitative Investing, Paleologo (2025) splits a book’s P&L in three (p. 335; PDF p. 361). He records the holdings at fixed times, for example each close; these are the snapshots. Position P&L is what the snapshot holdings earn if held unchanged until the next snapshot. A factor model splits it into factor P&L and specific P&L (section 2.4). Trading P&L is the rest: the actual P&L minus position P&L, that is, everything gained or lost by trading between snapshots. If a book is flat at every snapshot used by a daily attribution, position P&L is zero and all realised P&L falls into trading P&L. Intraday alpha, spread earned and execution costs remain there, even though the book carried inventory and factor risk between snapshots. Explaining them requires fills and the inventory path, or a model of how signals produced those trades.

1.2 Information and accounting interval

Ex ante analysis uses the book’s exposures (its sensitivity to each factor) and a risk model, both recorded before returns are realised; ex post analysis combines them with realised returns. In Active Portfolio Management, Grinold & Kahn (2000) ask that skill be judged on the bets the manager declared beforehand. If we fit new exposures after a loss, we test an explanation chosen after seeing the result, not those bets. Holdings-based attribution uses observed positions. Returns-based attribution, such as Sharpe’s style analysis in section 2.3, estimates exposures by regressing the portfolio’s returns on index returns, so it may not recover the positions actually held.

Currency P&L adds across non-overlapping periods on a consistent accounting basis. Percentage returns compound: the period return is $\prod_t(1+R_t)-1$ , not $\sum_t R_t$ . Daily percentage contributions therefore do not add up to the period return, and a linking rule is needed to rescale them. Summing currency P&L avoids this, but we still need to choose reference prices and a rule for sharing interaction terms such as the one in section 2.1.

1.3 Two kinds of indicator

Signals used to take positions. A forecast enters the trading algorithm and may change its trades. If we can rerun the policy, we can remove signals, run them alone or combine subsets of them. We must state the market path and risk convention used for these runs. Section 3 shares the resulting P&L between signals; section 4 explains why this differs from sharing the predictive information that indicators carry.

State indicators. Spread, volatility, order flow and inventory describe the conditions in which a loss occurs; the trading algorithm need not use them as inputs. First split the fills’ P&L into spread earned and gains or losses on inventory (section 5). Then explain how each component varies with the recorded states (section 4), and test whether the findings could arise by chance, accounting for the indicators, horizons and thresholds tried (section 6). Section 7 turns that sequence into a research project.

The same variable can serve both roles. Removing a forecast changes the trades in a rerun of the policy. Comparing fills’ average outcomes between high- and low-volatility periods is a statistical comparison; it does not by itself establish a causal effect.

↑ Sections

2. The classical answers

Area diagram for segment A of the example in this section, with weight on the horizontal axis and return in basis points on the vertical axis. The benchmark area is 0.4 × 600 = 240. The portfolio rectangle, 0.6 × 800 = 480, adds allocation 48 above the benchmark total return of 360 bp, a pale part (w_A − b_A)R_B = 72 below it, selection 80 and interaction 40. — Segment A of the example below as areas: weight times return, in basis points.

2.1 Allocation, selection, interaction

A manager’s portfolio is compared with a benchmark index. Both split their money across segments: groups of assets such as countries, sectors or asset classes. In the paper we follow below, the segments are countries. The manager makes two kinds of decision: how much to put in each segment, and which stocks to hold inside it. The benchmark holds the stocks of each segment in index proportions, so the portfolio’s stocks inside a segment can earn more or less than the benchmark’s stocks in the same segment.

Take the synthetic example used in the figure below, with two segments, A and B. Over the period, segment A did well (the benchmark’s stocks in A returned 600 bp) and segment B did less well (200 bp). The portfolio puts 60% in A and 40% in B; the benchmark puts 40% in A and 60% in B. Inside A, the portfolio’s stocks returned 800 bp, better than the benchmark’s 600 bp; inside B, they returned 100 bp, worse than the benchmark’s 200 bp. Overall the portfolio earns 520 bp and the benchmark 360 bp: an active return of 160 bp. We want to know how much of it comes from the weights given to each segment and how much from the stocks chosen inside each segment.

In their paper on measuring non-US equity portfolio performance, Brinson & Fachler (1985) split active return into allocation, selection and interaction. For each segment $i$ , write $w_i$ for its weight in the portfolio and $b_i$ for its weight in the benchmark. Write $r_i$ for the return earned by the portfolio within that segment, and $r_i^B$ for the segment’s benchmark return. Both sets of weights sum to one. The total returns are

R_P=\sum_i w_i r_i, \qquad R_B=\sum_i b_i r_i^B.

Each total weights the segment returns by the amounts invested there. In our example, the portfolio return is $0.6\times800+0.4\times100=520$ bp; the benchmark return is $0.4\times600+0.6\times200=360$ bp.

Start with their difference:

R_P-R_B=\sum_i(w_i r_i-b_i r_i^B).

This compares the actual portfolio with the benchmark. To separate weights from returns, add and subtract $w_i r_i^B$ and $b_i r_i$ inside each segment’s contribution:

\begin{aligned} w_i r_i-b_i r_i^B &=(w_i r_i^B-b_i r_i^B)\\ &\quad +(b_i r_i-b_i r_i^B)\\ &\quad +(w_i r_i-w_i r_i^B\\ &\qquad -b_i r_i+b_i r_i^B). \end{aligned}

The first line changes only the weights, using benchmark segment returns. The second changes only the segment returns, using benchmark weights; the last two lines account for changing both together.

Because the weight differences sum to zero, we can subtract the total benchmark return from the segment return in the first term without changing its sum:

\sum_i(w_i-b_i)R_B=0.

An overweight in one segment must be offset by underweights elsewhere. This gives the allocation term of Brinson and Fachler (1985)

A_i=(w_i-b_i)(r_i^B-R_B).

Allocation measures the gain from holding more of a segment that beat the benchmark overall, or less of one that lagged it. It uses the segment’s benchmark return, so it does not yet credit the investments chosen within that segment.

S_i=b_i(r_i-r_i^B).

Selection measures the gain from beating the segment benchmark, applied to the benchmark’s weight in that segment.

I_i=(w_i-b_i)(r_i-r_i^B).

Interaction is the extra gain or loss from applying that same return difference to the overweight or underweight. Together the three terms give

R_P-R_B=\sum_i(A_i+S_i+I_i).

Every part of the active return is included. For segment A, the calculations in basis points are

\begin{aligned} A_A&=(0.6-0.4)(600-360)=48,\\ S_A&=0.4(800-600)=80,\\ I_A&=(0.6-0.4)(800-600)=40. \end{aligned}

The portfolio overweights A, whose benchmark return is above the overall benchmark, and also beats the benchmark within A. All three effects are positive. For B,

\begin{aligned} A_B&=(0.4-0.6)(200-360)=32,\\ S_B&=0.6(100-200)=-60,\\ I_B&=(0.4-0.6)(100-200)=20. \end{aligned}

Underweighting the lower-return segment helps allocation. Selection is negative because the investments in B lag its benchmark; interaction is positive because the underweight reduces that loss. Across both segments, allocation is 80 bp, selection 20 bp and interaction 60 bp: together 160 bp.

In their study Determinants of Portfolio Performance, Brinson et al. (1986) use the segment benchmark return directly for allocation:

A_i^{\mathrm{1986}}=(w_i-b_i)r_i^B.

Compared with the 1985 term, this adds $(w_i-b_i)R_B$ to each segment’s allocation. Those additions sum to zero, so total allocation stays the same even though its segment entries change.

A separate choice is to measure selection at portfolio weights:

S_i^{\mathrm{portfolio}}=w_i(r_i-r_i^B)=S_i+I_i.

This folds interaction into selection; it must no longer be reported as an additional contribution. Both conventions account for the full active return. Checking that the numbers add up cannot tell us which convention to use: we must decide which comparison we want each label to describe.

Use the figure to change a weight or return, then switch conventions while keeping those inputs fixed. The bars show the two segments and their totals in basis points. (Brinson & Fachler, 1985) (Brinson et al., 1986)

One active return, several attribution conventions

Synthetic two-segment portfolio. Returns and effects are in basis points (100 bp = 1 percentage point); both portfolios remain fully invested.

Synthetic, single-period returns. Each weight slider sets segment A; segment B gets the remainder.

Allocation convention

Interaction convention

Segment A · effects in basis points

Segment B · effects in basis points

Total · effects in basis points

Switching 1985 ↔ 1986 changes segment allocations; their sum is unchanged. Absorbing interaction into selection changes those totals by convention.

What to notice

Switch between 1985 and 1986 and compare the segment bars with the total allocation bar. Then absorb interaction into selection: this changes the total selection bar as well. Neither switch changes the active return.

Explaining variation in returns does not tell us the share of money earned through investment policy. Brinson, Hood and Beebower studied 91 pension plans over 40 quarters, 1974–1983. For each plan they built a policy portfolio: its long-run mix of stocks, bonds and cash, estimated from its 10-year average holdings, invested in passive indices. Its return is the policy return. The plan’s active return is its actual return minus that policy return.

The average $R^2$ of each plan’s quarterly return regressed on its policy return is 93.6% (Table VII, p. 6). A plan’s returns move with its policy returns, but that does not tell us what share of its earnings comes from policy or why one plan earns more than another. Across plans, the standard deviation of policy returns is 0.22 percentage points a year, compared with 1.45 for active returns (Table VI, p. 6). Adding a constant loss to a plan’s quarterly returns lowers its mean and leaves $R^2$ unchanged: the regression intercept absorbs the constant. (Brinson et al., 1986)

Increase the differences between plans’ average active returns. Their mean returns move further apart, while their time-series fits stay unchanged: adding a constant to a plan’s returns does not change its $R^2$ . (Brinson et al., 1986)

A high R² does not explain differences in return levels

Synthetic, fixed-seed additive panel: 91 plans, 40 quarters. The active mean dispersion changes; common policy shocks and quarterly active noise stay fixed.

Synthetic plans with the same quarterly policy shock. The slider changes permanent active-return differences across plans, holding time-series noise fixed.

Average time-series R² · actual versus policy

Cross-plan dispersion · annual arithmetic mean returns

Policy: blue · active: orange · actual: green. Each dot is a synthetic plan.

R² measures fluctuations within a plan, not the level of its mean return. Population dispersion excludes estimation noise; finite-sample dispersion includes it. Annualization is four times the quarterly mean, without compounding.

What to notice

This is the distinction between Tables VI and VII of Brinson, Hood and Beebower (1986, p. 6). The simulation uses an additive design with policy and active population mean dispersions of 0.22 and 1.45 percentage points per year at the default setting, with quarterly noise calibrated to a population $R^2$ of 93.6%. These are synthetic arithmetic returns, not a reproduction of Brinson, Hood and Beebower’s compound-return tables. They report 1.45 pp as a cross-plan realised dispersion; here it is a population parameter. The $R^2$ and dispersions measured in this simulated panel differ from those population values because the sample is finite.

2.3 A fitted exposure need not be a holding

In a paper on asset allocation and performance measurement, Sharpe (1992) fits a style: non-negative asset-class weights summing to one by minimising the variance of the difference from a fund’s returns. The mean difference is unconstrained. We can read the fitted weights as a long-only mix, though the fund need not have held it. To assess selection later, Sharpe estimates the style from earlier observations, keeping the data used to estimate the benchmark separate from those used to evaluate performance.

2.4 Factor return and factor risk

A linear factor model writes the vector of asset returns at time $t$ as

r_t=B_t f_t+u_t.

Here $B_t$ contains asset exposures to the factors, $f_t$ contains the factor returns, and $u_t$ contains the asset returns left unexplained by them. For dollar holdings $h_t$ (the amount held in each asset) over the return interval, position P&L is

h_t^\top r_t=(B_t^\top h_t)^\top f_t+h_t^\top u_t.

The vector $B_t^\top h_t$ gives the book’s factor exposures. Factor exposure times factor return gives each factor’s contribution; the remaining term is specific P&L. These terms add up even when factors are correlated. Correlation affects risk and how we interpret the factors, but adds no interaction term to this linear return identity. Replacing the factors by linear combinations of them (rotating the factor basis) can change the individual contributions while preserving the total. (Grinold & Kahn, 2000; Paleologo, 2025)

To attribute risk, take a fixed covariance matrix $\Sigma$ and position vector $x$ . Portfolio volatility is

\sigma(x)=\sqrt{x^\top\Sigma x}.

This combines the positions with their variances and covariances. Volatility scales with the book: multiplying every position by $c>0$ multiplies $\sigma$ by $c$ . For positive portfolio volatility, the derivative with respect to a position is $(\Sigma x)_i/\sigma(x)$ . Euler’s theorem then gives

\sigma(x)=\sum_i x_i\frac{\partial\sigma}{\partial x_i}, \qquad \mathrm{RC}_i=\frac{x_i(\Sigma x)_i}{\sigma(x)}.

In Hot Spots and Hedges, Litterman (1996) uses this property to write risk as a sum of position times marginal risk. Here $\mathrm{RC}_i$ is position $i$ ‘s contribution to volatility. Each position’s contribution depends on its covariance with the whole book, as well as its size and stand-alone volatility. A hedge can contribute negatively, and the largest position need not contribute the most risk.

In the next figure, change the size, volatility and correlation of two positions and watch their risk contributions. A larger position can contribute less risk than a smaller one. The decomposition allocates volatility; it does not allocate the realised loss on a particular day. (Litterman, 1996)

Positions and their contributions to risk

Synthetic two-position portfolio. Position sizes, annual return volatilities and return correlation are adjustable; all risk values use the same capital units.

Synthetic linear positions, in units of capital. Annual return volatilities and correlation define the covariance matrix. Positions may be short.

Position sizes · capital units

Euler contributions · capital units of annual volatility

Position	Marginal risk per unit	Risk share

Component = position × marginal risk. Components sum to portfolio volatility; a negative component hedges the portfolio. Removing a whole position is a different calculation.

What to notice

The default synthetic book has a larger position in A and a larger risk contribution from B. A negative contribution means that increasing that position slightly reduces risk. The derivative describes a small change in the position, so it does not measure how much risk would change if we deleted it. At zero portfolio volatility, the derivative is undefined.

2.5 Estimated components have uncertainty

In this model, factor returns are estimated each period by a cross-sectional regression of asset returns on the exposures $B_t$ . Paleologo treats these estimates as noisy and calculates the resulting uncertainty in factor P&L. Specific P&L is what remains of position P&L after factor P&L, so an error in factor P&L appears with the opposite sign in specific P&L, and the two cancel in the total. The report can add up exactly even when we are unsure how much P&L belongs to each component. His standard-error calculation assumes a fixed model, normal returns and factor returns estimated by cross-sectional weighted least squares; extending it to changing intraday exposures requires further work. (Paleologo, 2025)

↑ Sections

3. Attributing a P&L to signals

A portfolio receives several forecasts, but it trades one position. The construction algorithm combines the forecasts with risk, transaction costs and position limits. A signal that asks for a purchase may reinforce another signal, offset it, or leave the portfolio unchanged. To attribute the resulting P&L, we need to know how the signals become trades.

If we can rerun the algorithm, we can compare its outputs with different signals as inputs. If we also know its optimality equations, we can split the joint solution directly. Coalition methods need more reruns, but do not need those equations. Each method uses different information and gives a different meaning to a signal’s contribution.

Pipeline from signals through the optimiser, with its risk, costs and constraints, to positions and P&L. Above the pipeline, stand-alone proxies, leave-one-out and Shapley rerun the construction with subsets of signals and need a policy we can rerun. Below it, Durin 2014 and Grinold–Easton use the optimality equations and multipliers at the joint optimum and need a known optimiser. — What each method in this section needs: reruns of the construction, or the optimiser’s own equations.

3.1 Stand-alone portfolios as proxies

A stand-alone proxy method used in practice reruns the same construction algorithm with one signal active at a time. Each run keeps the construction rules, including costs and constraints, and gives that signal its own risk target, a chosen level of predicted portfolio volatility. In the simulation below, the optimiser adjusts risk aversion each period to reach that target when the constraints permit it. The resulting P&L series includes the trades generated by the signal: a fast signal’s proxy includes its turnover, and a slow signal’s proxy includes the positions it carries. We fit these P&Ls, not the raw forecasts.

Write the total P&L in period $t$ as $y_t$ , and the proxy for signal $k$ , at risk target $\sigma_k$ , as $p_{k,t}(\sigma_k)$ . Fit non-negative least squares without an intercept:

\widehat\beta =\arg\min_{\beta_k\geq0} \sum_t\left(y_t-\sum_{k=1}^{K}\beta_k p_{k,t}(\sigma_k)\right)^2.

Here $p_{k,t}$ abbreviates $p_{k,t}(\sigma_k)$ . The contribution assigned to signal $k$ is $\widehat\beta_k p_{k,t}(\sigma_k)$ . Keep the residual as a separate term:

e_t=y_t-\sum_k\widehat\beta_k p_{k,t}.

There is no constant term to absorb an average gain or loss. The non-negative coefficients prevent the fit from reversing a proxy’s sign. They can still assign a negative contribution to a proxy that loses money.

The method then updates each risk target by

\sigma_k\leftarrow\widehat\beta_k\sigma_k.

It multiplies each target by the fitted coefficient, regenerates the proxies and fits again. A coefficient below one lowers the proxy’s risk target; a coefficient above one raises it. At a fixed point, each surviving proxy enters the fit unscaled. A coefficient of zero drops the proxy. Scaling a completed proxy can give a different path from rerunning the nonlinear construction at a new risk target, so we need to repeat the regression.

The coefficients tell us how to weight these proxies. We cannot read them directly as independent risk budgets in the combined portfolio, because signals can be correlated, trades can net and constraints can bind. In his paper on signal-wise attribution for constrained optimisation, Durin (2014) describes regressions on stand-alone portfolio P&Ls as existing practice, but rejects them as a general solution: binding constraints can leave a large residual (p. 3). His attribution uses the optimiser’s internal equations to address this limitation.

The residual shows how much of the joint strategy these separate runs fail to reproduce. The construction algorithm can produce joint trades whose P&L no fixed non-negative mixture of stand-alone P&Ls reproduces. Assigning the residual to the signals would make the contributions add up without explaining any more of the joint P&L.

Even a small residual only tells us that the proxy P&Ls fit the observed total in this sample. Similar proxies can still have non-unique coefficients. A good fit neither establishes causality nor tells us what would happen if we removed a signal.

3.2 Removing one signal from the full portfolio

Leave-one-out attribution runs the construction with all signals and then with each signal omitted. In his book Quantitative Portfolio Management, Isichenko (2021) describes this method for forecast attribution: the deterioration of the historical simulation measures the forecast’s marginal contribution, and the contributions need not sum to the total because the combination is nonlinear (p. 134; PDF p. 151).

Let $N$ be the full set of signals, and let $v(S)$ be the P&L obtained by running only the signals in a subset $S$ , under a specified risk convention. The leave-one-out contribution is

L_k=v(N)-v(N\setminus\{k\}).

A stand-alone run measures what a signal does by itself; leave-one-out measures what disappears when we remove it from the existing combination. A signal that is redundant with another can have a substantial stand-alone P&L and a small leave-one-out contribution. A signal that helps the combined forecast clear a trading threshold can have the reverse pattern.

Neither method guarantees that contributions add up to the total. If several signals jointly enable a trade, removing any of them may remove the same gain, so adding the removal effects can count that gain more than once. For both methods, replay the same market path with the same cost model and treatment of initial positions. State which signals and risk allocations change between runs.

3.3 Using the optimiser’s internal equations

Durin’s method needs a portfolio built by a known optimiser. We need its optimal positions, active constraints and Lagrange multipliers, which measure how the optimum changes when a bound is relaxed. The setting below is convex, and its optimality conditions become linear in the forecasts once we handle the constraints at the joint optimum. A heuristic position rule, a black-box model or a learned policy does not supply these equations; sections 3.1, 3.2, 3.4 and 7 discuss approaches for those cases.

The calculation solves for the next trade. Durin also covers the dynamic model of Gârleanu and Pedersen through modified forecasts and matrices (pp. 4–5), but imposing constraints on that next trade while retaining the unconstrained value function is an approximation. He does not solve the full constrained dynamic programme (p. 7). Combinatorial and non-convex constraints, such as round lots or a minimum trade size, are outside the treatment (pp. 7, 14).

Let $x_{t-1}$ be the previous position, $\Delta x_t$ the trade, and $g_t^k$ the expected return supplied by signal $k$ for the next period; the combined forecast is $G_t=\sum_k g_t^k$ . After the trade the position is $x_t=x_{t-1}+\Delta x_t$ . In the simplest version, the static model, the optimiser chooses the trade that maximises the expected return of the new position, minus a risk penalty on that position, minus a quadratic trading cost:

x_t^\top G_t-\tfrac{\gamma}{2}\,x_t^\top\Sigma x_t-\tfrac12\,\Delta x_t^\top\Lambda\Delta x_t .

Here $\Sigma$ is the covariance of returns, $\gamma$ the risk aversion and $\Lambda$ the matrix of quadratic trading costs. Replace $x_t$ by $x_{t-1}+\Delta x_t$ and drop the terms that do not depend on the trade. What remains is

\begin{gathered} \Delta x_t^\top G_t-\tfrac12\,\Delta x_t^\top Q\Delta x_t-\Delta x_t^\top Px_{t-1},\\ Q=\Lambda+\gamma\Sigma,\quad P=\gamma\Sigma . \end{gathered}

$Q$ combines the cost of trading and the risk the trade adds; $P$ measures how the trade interacts with the risk of the position already held. Durin writes the problem in this form (his Eq. 5, p. 4; static matrices in Eq. 7, p. 5) because the dynamic trading model of Gârleanu and Pedersen (2013), on which his paper builds, leads to the same expression, with other matrices $Q$ and $P$ and forecasts adjusted for how fast each signal decays. Setting the derivative with respect to $\Delta x_t$ to zero gives the optimality equation (Eq. 8, p. 6)

Q\Delta x_t+Px_{t-1}=\sum_k g_t^k.

The equation is linear in the forecasts. Starting from $x_0^k=0$ for every signal, we compute its trade and then update its position (Eqs. 11–12, p. 6):

\Delta x_t^k=Q^{-1}(g_t^k-Px_{t-1}^k).

The trade responds to this signal’s forecast and the position it already holds.

x_t^k=x_{t-1}^k+\Delta x_t^k, \qquad x_t=\sum_k x_t^k.

The component positions add up to the actual position (Eq. 9, p. 6); a nonzero initial position must also be split consistently. Here and below, the inverses require nonsingular matrices; zero-bound cases use the limits described in the optional derivation.

Let $r_{t+1}$ be the vector of realised returns over the next holding period. The gross P&L assigned to signal $k$ is

\mathrm{PnL}_{t+1}^k=(x_t^k)^\top r_{t+1}.

Multiplying every component by the same returns gives an exact gross P&L split (the construction following Eq. 12, p. 6). For risk, write $R$ for portfolio variance and $R_k$ for the share assigned to signal $k$ (Eq. 13, p. 6):

R=x_t^\top\Sigma x_t =\sum_k (x_t^k)^\top\Sigma x_t =\sum_k R_k.

Each share measures the component’s covariance with the full portfolio, including its interaction with the other signals. For quadratic cost $C$ , the analogous allocation is (Eq. 14, p. 7)

\begin{aligned} C&=\tfrac12\Delta x_t^\top\Lambda\Delta x_t=\sum_k C_k,\\ C_k&=\tfrac12(\Delta x_t^k)^\top\Lambda\Delta x_t. \end{aligned}

A component receives its share of the cost of the joint trade, including the savings when component trades offset each other.

Two ways to handle constraints. Grinold and Easton’s attribution keeps each constraint as a separate source. In Durin’s notation, the optimality equation becomes (Eq. 17, p. 7)

Q\Delta x_t+Px_{t-1}=G_t+\sum_c\lambda_t^c.

Here $c$ indexes constraints and $\lambda_t^c$ is the signed vector contribution of a constraint’s multiplier to the equation, including its exposure direction. Feeding each of these sources through the same recursion builds a constraint portfolio alongside the signal portfolios. This tells us how much P&L is associated with the constraint, but leaves open which signals should bear its effect (pp. 7–8) (Durin, 2014).

An active position constraint contributes extra quadratic risk at the optimum; a trade constraint contributes extra quadratic cost. Write $\bar Q_t$ and $\bar P_t$ for the resulting matrices. The signal-wise recursion is then (Eq. 24, p. 10, with Eqs. 10–12, p. 6)

\begin{aligned} \Delta x_t^k&=\bar Q_t^{-1}(g_t^k-\bar P_t x_{t-1}^k),\\ x_t^k&=x_{t-1}^k+\Delta x_t^k. \end{aligned}

Every signal uses the same effective risk and cost, so its trades include the effect of the binding constraints. A binding trade limit slows trading in the constrained direction. A position limit reduces the aim in that direction, but can make adjustment towards that smaller aim faster. In Durin’s value/momentum example, even the signal pointing away from the forbidden short position is affected; the adjustment also depends on each signal’s position from previous trades (Eq. 34, p. 13).

We evaluate the matrices at the full portfolio’s optimum, so the component positions are shares of that joint solution. Running each signal alone and recomputing its multipliers can give different positions. Removing a signal can change the matrices and every remaining signal’s contribution.

Optional From active bounds to effective matrices: Squared constraints, solver multipliers and zero bounds

Let $v_c$ be the vector defining an exposure and $M_c$ its active bound. For example, a sector exposure uses a vector with an entry of one for each stock in that sector and zero elsewhere. At an active upper position bound, $v_c^\top x_t=M_c$ . Durin uses the squared constraint $(v_c^\top x_t)^2\leq M_c^2$ and a multiplier $\eta_c$ for it (Eqs. 26–29, p. 11). For a lower bound, use the bound that is active at the optimum. This construction matches the optimality equation there; squaring an arbitrary one-sided bound need not preserve its whole feasible set. The nonnegative-penalty interpretation applies to bounds that reduce the trade or position, in particular when the admissible set contains zero (p. 14).

Dropping terms constant in the trade, the augmented objective for one position constraint is (Eq. 29, p. 11)

\begin{aligned} F={}&-\tfrac12\Delta x_t^\top Q\Delta x_t -\Delta x_t^\top Px_{t-1}\\ &+\Delta x_t^\top G_t -\eta_c x_t^\top(v_c\otimes v_c)x_t. \end{aligned}

The last term penalises squared exposure; $v_c\otimes v_c=v_cv_c^\top$ is an outer product. With $A_c=2v_c\otimes v_c$ , differentiation gives (Eqs. 29–31, pp. 11–12)

\begin{gathered} (Q+\eta_c A_c)\Delta x_t\\ =G_t-(P+\eta_c A_c)x_{t-1}. \end{gathered}

Thus one position constraint adds $\eta_c A_c$ to both matrices; in the static model, it adds $\eta_c A_c/\gamma$ to the covariance matrix. A squared trade constraint instead adds a term $-\mu_c\Delta x_t^\top(v_c\otimes v_c)\Delta x_t$ to the objective, where $\mu_c$ is its multiplier, so it changes only $Q$ (p. 12).

For several constraints, let $\mathcal T$ index trade constraints and $\mathcal P$ position constraints. Including the time dependence of the multipliers, Eq. 31 (p. 12) gives

\begin{aligned} \bar Q_t&=Q+\sum_{c\in\mathcal T}\mu_{t,c}A_c\\ &\quad+\sum_{c\in\mathcal P}\eta_{t,c}A_c,\\ \bar P_t&=P+\sum_{c\in\mathcal P}\eta_{t,c}A_c. \end{aligned}

Only position constraints change the term involving the previous position; inactive constraints have zero multipliers. For individual stock bounds, these additions are diagonal, as in Durin’s long-only example (Eq. 32, p. 12).

The solver supplies multipliers for the original linear bounds. Denote one such scalar by $\lambda^c$ ; unlike $\lambda_t^c$ in Eq. 17 above, it does not include the exposure vector. Let $\epsilon$ equal $+1$ for an upper bound and $-1$ for a lower bound in the maximisation convention. The attribution multiplier is (Eqs. 35–37, p. 14)

\eta_c=\frac{\epsilon\lambda^c}{2M_c}.

This converts sensitivity to the original bound into sensitivity to its square; the same conversion gives $\mu_c$ for a trade bound.

At a zero bound this division cannot be used directly. For one trade constraint, define $h_t^k=g_t^k-Px_{t-1}^k$ and let $v$ be its exposure vector. Applying Durin’s rank-one inverse to each signal gives (Eqs. 39–41, p. 15)

\begin{aligned} \Delta x_t^k &=\left[Q^{-1} -\alpha\frac{Q^{-1}vv^\top Q^{-1}}{v^\top Q^{-1}v}\right]h_t^k,\\ \alpha&=\frac{2\eta v^\top Q^{-1}v} {1+2\eta v^\top Q^{-1}v}. \end{aligned}

As the active bound tends to zero and $\eta$ tends to infinity, $\alpha$ tends to one: the trade is projected onto $v^\top\Delta x_t^k=0$ , along $Q^{-1}v$ (Eq. 42, pp. 15–16). For a position bound, the trade expression has the additional term (Eqs. 43–44, p. 16)

-\alpha\frac{Q^{-1}vv^\top}{v^\top Q^{-1}v}\,x_{t-1}^k.

It removes the previous position’s exposure too, so the limiting component position satisfies $v^\top x_t^k=0$ . The infinite multiplier therefore has a finite, well-defined effect on the attributed trades and positions.

Spread costs and the no-trade band. Take a single asset in the static model. Its return variance is $\sigma^2$ , and $\Lambda$ now denotes its scalar quadratic cost coefficient (Durin calls this coefficient $\lambda$ in Eqs. 49–60). Charge a positive spread $\lambda_0$ per unit traded, adding $-\lambda_0|\Delta x_t|$ to the objective. For a nonzero joint trade, the effective quadratic-cost multiplier is (Eqs. 47–48, pp. 16–17)

\eta_t=\frac{\lambda_0}{2|\Delta x_t|}.

Smaller joint trades therefore face a larger effective quadratic cost. The optimality equation becomes (Eq. 52, p. 17)

(\gamma\sigma^2+\Lambda+2\eta_t)\Delta x_t =G_t-\gamma\sigma^2x_{t-1}.

Spread costs enter the trade coefficient, leaving the previous-position coefficient unchanged. Each signal uses the same $\eta_t$ , with $g_t^k$ and $x_{t-1}^k$ on the right. A nonzero trade exists only when (Eqs. 50–54, p. 17)

|G_t-\gamma\sigma^2x_{t-1}|>\lambda_0.

Inside this band, including its boundary, the combined forecast does not justify paying the spread. The limiting effective cost is infinite, and every signal’s attributed trade is zero. Durin calls this a voting rule: no signal trades unless the signals together push the forecast out of the band (p. 18). Two signals can cancel inside it, or together clear it. Stand-alone portfolios face their own thresholds and histories, so a regression coefficient cannot generally align their trading dates with the joint portfolio’s dates.

For several assets, spread cost is assigned using the direction of each joint trade. With $i$ indexing assets, its allocation to signal $k$ is (Eq. 22, p. 9, written as a positive cost to subtract)

C_{\mathrm{spread},t}^k =\lambda_0\sum_i\Delta x_{t,i}^k\, \operatorname{sgn}(\Delta x_{t,i}).

Here $\operatorname{sgn}$ is the sign function, with value zero at zero. A component trading against the full trade can receive a negative cost allocation because it reduces the net trade to execute. Subtracting these spread allocations and the quadratic costs from each component’s gross P&L gives an exact net split. The effective costs determine the component trades; they are not additional cash costs to subtract.

The same construction handles a turnover constraint, with its multiplier in place of the spread coefficient. Financing costs and leverage constraints penalise absolute positions instead of absolute trades: they become effective risk terms and can force every component’s position to zero (§2.3, pp. 8–9; §3.4, Eqs. 55–60, p. 18) (Durin, 2014).

Move the signals and the spread in the next figure to see when the combined forecast produces a trade.

Spread costs create a no-trade band

Move the two forecasts. The combined forecast must overcome both the spread and the risk of the existing position before any signal trades.

Drag the point, click the plane, or move the signal sliders. All values use synthetic, dimensionless units.

Variance σ² = 1; quadratic cost Λ = 1; no position bound. Previous position is split equally between signals: x₋₁¹ = x₋₁² = x₋₁/2. This is one step with a prescribed prior split, not a full trading path.

What to notice

Inside the band, both attributed trades are zero, even when one signal alone would trade. Outside it, Durin applies the same effective quadratic cost to both signals. Their attributed trades add up to the actual trade, even if one goes in the opposite direction. Here we set the prior component positions equal. Along a full path, we would take them from the previous attribution step. This is the static model of Durin (2014), pp. 17–18, Eqs. 49–54 (Durin, 2014).

For a coalition game, rerun the algorithm with different subsets of signals. Define $v(S)$ as the P&L of the construction algorithm with only signals in $S$ , with $v(\emptyset)=0$ . A stand-alone proxy is its singleton value $v(\{k\})$ , provided the risk convention matches. If we update each proxy’s risk target, we need to check this correspondence again: changing the risk convention changes the game.

The Shapley value credits a signal with its marginal contribution averaged over all orders in which the signals could enter. Here $K$ is the number of signals in $N$ . Write the marginal contribution as

\Delta_k v(S)=v(S\cup\{k\})-v(S).

This is the extra P&L when signal $k$ joins the signals already in $S$ . The allocation is

\phi_k=\sum_{S\subseteq N\setminus\{k\}} \frac{|S|!\,(K-|S|-1)!}{K!} \Delta_k v(S).

The weight of $S$ is the probability that, when the $K$ signals enter in a uniformly random order, the signals already present when $k$ arrives are exactly those in $S$ : $|S|!$ orders of the signals before $k$ , times $(K-|S|-1)!$ orders of those after, out of $K!$ .

The singleton is the contribution when $k$ enters first; leave-one-out is its contribution when it enters last. Shapley includes the intermediate coalitions. For each order, the successive increments add up to the full P&L, so their average satisfies $\sum_k\phi_k=v(N)$ .

In an additive game the singleton values already sum and Shapley returns them. Otherwise $\phi_k-v(\{k\})$ records the interaction allocated to signal $k$ relative to its stand-alone value. This allocation uses the whole coalition game, whereas proxy regression only fits combinations of singleton series. Exact enumeration requires $2^K$ coalition values, so the extra information has a computational cost.

Follow a signal from its singleton value through its marginal contributions to different coalitions in the next figure.

Stand-alone, leave-one-out and Shapley

Edit an illustrative three-signal game. Select a signal to inspect each marginal contribution and its Shapley weight.

Synthetic teaching game, in integer P&L units. Edit any non-empty coalition. The empty coalition stays at zero. These values illustrate an allocation rule; they are not measured returns.

Inspect signal

Each edge adds one signal. Highlighted edges add the selected signal.

What to notice

Stand-alone measures what a signal earns by itself; leave-one-out measures what the full portfolio loses when we remove it. Shapley averages the signal’s marginal contribution over all orders in which signals could join. With three signals, each end coalition gets weight $1/3$ and each intermediate coalition gets weight $1/6$ . The shares sum to the full value. The gap between the full value and the sum of singleton values is the total interaction, not a regression residual.

You can edit these integer values, which were chosen for a teaching example. To apply this to trading, first decide how to rebuild each coalition’s portfolio, including its risk convention. The weighting follows Shapley’s allocation formula (Shapley, 1953); these P&L values are not a claim about a measured strategy.

The synthetic comparison uses two games. In the first, every coalition gets the same risk target. In the second, every coalition uses the risk aversion $\gamma_t$ that the full portfolio used at each date. This sequence is the full portfolio’s risk-aversion path, labelled “fixed γ” in the figure. With fixed risk aversion, no position bound and no spread, positions are additive. Gross P&L is linear and quadratic costs create pairwise cross terms; Shapley and Durin then coincide, including the allocation of quadratic costs.

3.5 What the synthetic comparison finds

All methods receive the same synthetic signals and returns. Durin gives an exact internal allocation for this optimiser, so we use it as a reference. Agreement with Durin does not make a contribution uniquely true.

The synthetic desk. A and B label different scores, not asset classes. The planted coefficients describe how the scores predict returns; they are not P&L shares after trading costs and constraints.

Base setup	Value
Assets; periods	50; 1,000
Asset covariance	Diagonal; volatility 1 per period
Portfolio risk target	1 per period
Scores	AR(1), unit marginal variance
Innovation correlation	0.3 between signals
Quadratic cost	$\Lambda=\gamma_{\rm ref}\Sigma$
Spread per unit traded	$0.5\,\mathrm{std}(G)$
Absolute position bound	$1.5/\sqrt{50}$ per asset
Study seeds	20 per configuration

Signal	Persistence	Planted coefficient
Fast A	0.5	0.06
Slow A	0.95	0.05
Fast B	0.5	0.04
Slow B	0.95	0.03
Mid A (six signals only)	0.8	0.05
Mid B (six signals only)	0.8	0.04

Returns are the sum of the planted signal components plus independent unit-volatility noise. Here $G$ is the combined forecast and $\gamma_{\rm ref}=\sqrt{50}\,\mathrm{std}(G)$ is the no-cost, no-bound risk aversion at the unit risk target. The quadratic cash cost is $\tfrac12\Delta x^\top\Lambda\Delta x$ . The two-signal configuration keeps fast A and slow A; the six-signal configuration adds the two mid signals. Each configuration below changes only the named parameter from the base.

The tables and the bar figure average 20 seeds per configuration. The P&L paths use one seed, 20261029: the same model and configurations, with different random realisations.

Iterated proxies compared with Durin. Distance is half the sum of the absolute differences between the methods’ cumulative net P&L shares, computed within each seed and then averaged. Zero means the same split. The residual stays separate. “Residual variance” is its variance divided by the book’s P&L variance.

Configuration	Distance	Residual variance
Base	0.1851	41.60%
Zero spread	0.0288	8.44%
Spread ×1	0.4344	87.52%
Spread ×2	0.4921	98.63%
Two signals	0.0587	18.87%
Six signals	0.3176	73.62%
Correlation 0.9	0.3929	60.81%

Configuration	Zero coefficients	Converged
Base	0.10 / 4	100%
Zero spread	0.00 / 4	100%
Spread ×1	1.20 / 4	35%
Spread ×2	3.00 / 4	0%
Two signals	0.00 / 2	100%
Six signals	1.50 / 6	65%
Correlation 0.9	2.95 / 4	100%

Spread ×1 and ×2 mean one and two times $\mathrm{std}(G)$ . Zero coefficients are the average number of dropped or zero-P&L proxies, out of the number of signals. Convergence requires every live proxy’s coefficient to be within $10^{-3}$ of one within 50 iterations. A run with no live proxy does not count as converged.

The iterated fit approaches Durin when spread is small and there are few signals. Even at zero spread, the risk target rescales the portfolio each period: removing a signal changes the risk aversion applied to the others. The stand-alone runs therefore need not add up to the joint run.

Large spreads leave most of the book’s variation in the residual. Strong correlation creates a different problem: every run converges, but usually with only one proxy left. Convergence does not recover the dropped signals’ contributions. Identical proxy columns can also give non-unique NNLS solutions, so the fit need not split identical signals equally.

The one-step regression’s residual tends to move with the coalition interaction in section 3.4, $\phi_k-v(\{k\})$ , but need not equal it. Regression absorbs part of the interaction into its coefficients, and the non-negativity constraint also affects the residual. The residual measures what the proxies leave unexplained; the distance measures how the signal allocations differ.

Choose a scenario, then switch allocation methods to follow each signal’s cumulative net P&L. Compare the endpoints below the paths: the book is the same for every method within a scenario.

Synthetic P&L paths under each allocation

Seed 20261029, four configurations; the tables and bar figure average 20 seeds per configuration. The thick line is the total book; dashed lines show residuals or the separate spread and position-limit portfolios.

Synthetic net P&L · one fixed seed

Final cumulative allocation by method

All panels use the same horizontal scale. The vertical black mark is the total book.

Units: asset volatility and the portfolio risk target are both 1 per period. These are simulation P&L units, not currency amounts. Paths show every fifth period. The path scale fits the selected method; the endpoint panels share one scale within each scenario.

Both proxy fits use NNLS without an intercept; the iterated fit updates risk targets without damping. Leave-one-out removes each signal and reruns the remaining book at the same risk target. Shapley averages marginal contributions over all coalitions, either at that target or at the full book’s risk aversion at each date. The Grinold–Easton spread and position-limit lines are attributed portfolios, each net of its allocated costs; the spread line is not the cash spread bill. All values are synthetic simulation units, with asset volatility and the portfolio risk target set to one per period.

What to notice

With the large spread, slow A, fast B and slow B never trade alone: their forecasts stay inside the no-trade band. Fast A trades alone on just four dates and ends with a small loss. Its P&L moves slightly against the book; giving it a positive weight would worsen the NNLS fit. Both proxy methods therefore assign zero to every signal and leave the entire book in the residual.

This is the voting rule of section 3.3 seen from the stand-alone side: signals that cannot trade alone can help clear the spread together. Durin allocates the resulting joint P&L to all four signals on this path.

With the tight position limit, Grinold–Easton’s negative spread and position-limit portfolios offset large positive signal portfolios. Their sum follows the same book as Durin, whose four signal paths include these effects.

At high innovation correlation, the iterated fit keeps only fast A and leaves a positive residual. Its convergence on that surviving proxy does not recover the other signals’ contributions. These are single-path examples; the following figure compares results across seeds.

Select a parameter to vary. The bars show each method’s share of cumulative net P&L for the four baseline signals described above, with the residual on its own row.

The same synthetic P&L, several signal allocations

Select one swept parameter and one of its calculated values. Compare the undamped stand-alone proxy method, one-step regression, Durin, and two Shapley games.

One axis varies at a time; other parameters stay at the baseline configuration. Bars: mean share of cumulative net P&L. Whiskers: ±1 sample standard deviation across synthetic seeds, not confidence intervals.

What to notice

The proxy methods use non-negative least squares without an intercept. The iterative method updates risk targets without damping. A small cumulative residual does not imply a small residual variance. The iteration can converge on the remaining proxies without recovering every signal.

“Shapley target” reruns every subset at the same risk target; “fixed γ” uses the full portfolio’s risk aversion at each date (section 3.4). Whiskers show dispersion across synthetic seeds, not estimation intervals for a real desk. The shares average per-seed ratios over 20 seeds per configuration, with the residual kept separate.

Costs vary little over time in this experiment. At the base configuration, the standard deviation of costs is 2.49% of that of gross P&L. This does not make costs deterministic or establish their relative size on a real desk. Fitting gross P&L and attaching each proxy’s scaled costs changes the distance from 0.185 to 0.181 and the residual variance ratio from 41.6% to 42.4%. Separating costs does not remove the construction interactions.

The fixed-risk-aversion equality in section 3.4 does not extend automatically to the nonlinear portfolio. At the base configuration, the net-game Shapley distances from Durin are 0.045 with a risk target per coalition and 0.037 with the full portfolio’s risk-aversion path. Both come close to Durin, but neither matches it exactly. Shapley’s contributions add up to the full-coalition value of its own game. To reconcile the original ledger, the numerical rerun of that coalition must also reproduce the original portfolio.

These methods need reruns with specified changes to the inputs. Section 7 keeps this comparison separate from later work on causal adjustment and learned policies.

Code: reproduce the synthetic comparison

Download the executed notebook. It runs with NumPy, SciPy and cvxpy (Clarabel). The code covers the optimiser, Durin’s allocation and the other methods in section 3, then checks the figure’s baseline and large-spread paths and runs a dense-covariance example.

Optimiser

Reproduce the synthetic signal comparison

Install numpy scipy cvxpy clarabel, then run the cells in order. The notebook is self-contained: it needs no checkout or download. NumPy handles the diagonal calculation; CVXPY with Clarabel handles a full covariance matrix. The final cell reproduces both published fixed-seed paths exactly. A separate small dense example checks what changes when assets are correlated.

The formulas are those in section 3.3 and Bruno Durin’s Signal-wise performance attribution for constrained optimisation (2014). Its equation numbers below make the code’s conventions explicit. All values are synthetic. The fixed-seed examples are not the comparison tables’ averages across 20 seeds.

from dataclasses import dataclass
from math import factorial
import numpy as np
from scipy.optimize import nnls
import cvxpy as cp
import warnings
solver_events = []

1. Optimiser

Desk holds the signal coefficients, persistence, cost scales and position bounds. covariance returns either the identity matrix or a one-factor covariance with unit diagonal. The static problem uses $\Lambda=\lambda\Sigma$ (section 3.3; Durin Eq. 7), so this compact example supports a general covariance matrix but does not expose an independent cost matrix. gamma_ref, lam and lam0 keep the published scaling conventions.

@dataclass(frozen=True)
class Desk:
    n_assets: int = 50
    n_periods: int = 1000
    names: tuple = ("fast A", "slow A", "fast B", "slow B")
    ic: tuple = (0.06, 0.05, 0.04, 0.03)
    rho: tuple = (0.5, 0.95, 0.5, 0.95)
    corr: float = 0.3
    risk_target: float = 1.0
    box: float = 1.5
    lam_rel: float = 1.0
    spread_rel: float = 0.5
    cov: str = "diag"
    long_only: bool = False
    factor_share: float = 0.3

    def sigma2(self):
        return np.ones(self.n_assets)

    def innovation_corr(self):
        C = np.full((len(self.ic), len(self.ic)), float(self.corr))
        np.fill_diagonal(C, 1.0)
        return C

    def g_std(self):
        rho = np.asarray(self.rho, float)
        a = np.sqrt(1 - rho**2)
        S = self.innovation_corr() * np.outer(a, a) / (1 - np.outer(rho, rho))
        b = np.asarray(self.ic, float)
        return float(np.sqrt(b @ S @ b))

    def gamma_ref(self):
        return np.sqrt(self.n_assets) * self.g_std() / self.risk_target

    def lam(self):
        return self.lam_rel * self.gamma_ref()

    def lam0(self):
        return self.spread_rel * self.g_std()

    def bounds(self):
        hi = np.full(self.n_assets, self.box * self.risk_target / np.sqrt(self.n_assets))
        return (np.zeros(self.n_assets) if self.long_only else -hi), hi

    def factor_model(self):
        beta = np.linspace(0.8, 1.2, self.n_assets)
        f = np.sqrt(self.factor_share) * beta
        return f, 1.0 - f**2

    def covariance(self):
        if self.cov == "diag":
            return np.diag(self.sigma2())
        f, s2 = self.factor_model()
        return np.outer(f, f) + np.diag(s2)

simulate generates all AR(1) innovations, then all return noise, in the published order. Returns equal the sum of forecast components plus noise (Durin Eqs. 1–4). The initial scores use the first innovations, leaving a small cross-correlation transient for mixed persistence. The dense example draws its factor noise before its idiosyncratic noise.

def simulate(desk, seed):
    rng = np.random.default_rng(seed)
    K, T, N = len(desk.ic), desk.n_periods, desk.n_assets
    rho = np.asarray(desk.rho, float)[:, None]
    L = np.linalg.cholesky(desk.innovation_corr() + 1e-12 * np.eye(K))
    e = np.einsum("kl,ltn->ktn", L, rng.standard_normal((K, T, N)))
    z = np.empty((K, T, N))
    z[:, 0] = e[:, 0]
    a = np.sqrt(1 - rho**2)
    for t in range(1, T):
        z[:, t] = rho * z[:, t - 1] + a * e[:, t]
    if desk.cov == "diag":
        u = rng.standard_normal((T, N))
    else:
        f, s2 = desk.factor_model()
        u = rng.standard_normal((T, 1)) * f + rng.standard_normal((T, N)) * np.sqrt(s2)
    g = np.asarray(desk.ic, float)[:, None, None] * z
    return g, g.sum(0) + u

Solution stores positions, trades, risk aversion, spread subgradients and signed box multipliers. For CVXPY runs, dual_upper, dual_lower and primal_raw retain the solver outputs. In particular, nu remains the raw upper-minus-lower dual; it is not overwritten to force attribution sums.

@dataclass
class Solution:
    """Path of the algorithm for a batch of B runs."""
    x: np.ndarray
    dx: np.ndarray
    gamma: np.ndarray
    s: np.ndarray
    nu: np.ndarray
    at_bound: np.ndarray
    dual_upper: np.ndarray | None = None
    dual_lower: np.ndarray | None = None
    primal_raw: np.ndarray | None = None

_bcast_bounds repeats the same position bounds for each batched run (Durin Eq. 15).

def _bcast_bounds(lo, hi, B, N):
    lo = np.broadcast_to(np.asarray(lo, float), (B, N))
    hi = np.broadcast_to(np.asarray(hi, float), (B, N))
    return (lo, hi)

_trade_diag is the diagonal closed form: soft-threshold the trade and clip the position (section 3.3; Durin Eqs. 50–51).

def _trade_diag(G, xp, gamma, sig2, lam, lam0, lo, hi):
    """Closed-form optimal position for given gamma (B,) (Durin Eqs. 50-51, p. 17)."""
    g = gamma[:, None]
    q = (lam + g) * sig2
    z = G - g * sig2 * xp
    dstar = np.sign(z) * np.maximum(np.abs(z) - lam0, 0.0) / q
    xu = xp + dstar
    return (np.clip(xu, lo, hi), xu, z, q)

solve_diag uses 32 bisection steps on log risk aversion. It retains the published floor $10^{-8}\gamma_{ref}$ when the target is out of reach. A fixed risk-aversion path skips this search. This fast implementation is retained for exact parity with the published paths.

def solve_diag(desk: Desk, G: np.ndarray, sigma_target: np.ndarray | None=None, gamma_path: np.ndarray | None=None, lo: np.ndarray | None=None, hi: np.ndarray | None=None, n_bisect: int=32) -> Solution:
    """Run A on a batch of combined signals G (B, T, N), diagonal Sigma."""
    (B, T, N) = G.shape
    sig2 = desk.sigma2()
    (lam, lam0) = (desk.lam(), desk.lam0())
    if lo is None or hi is None:
        (lo, hi) = desk.bounds()
    (lo, hi) = _bcast_bounds(lo, hi, B, N)
    x = np.zeros((B, T, N))
    gam = np.zeros((B, T))
    xp = np.zeros((B, N))
    g_ref = desk.gamma_ref()
    (g_min, g_max) = (1e-08 * g_ref, 100000000.0 * g_ref)
    g_prev = np.full(B, g_ref)
    if sigma_target is not None:
        tgt2 = np.broadcast_to(np.asarray(sigma_target, float), (B,)) ** 2

    def risk2(gm):
        xx = _trade_diag(G[:, t], xp, gm, sig2, lam, lam0, lo, hi)[0]
        return (sig2 * xx ** 2).sum(1)
    for t in range(T):
        if gamma_path is not None:
            gm = gamma_path[:, t]
        else:
            a = np.maximum(g_prev / 30, g_min)
            b = np.minimum(g_prev * 30, g_max)
            for _ in range(10):
                low_a = risk2(a) < tgt2
                a = np.where(low_a & (a > g_min), np.maximum(a / 1000.0, g_min), a)
                high_b = risk2(b) > tgt2
                b = np.where(high_b & (b < g_max), np.minimum(b * 1000.0, g_max), b)
                if not ((low_a & (a > g_min)).any() or (high_b & (b < g_max)).any()):
                    break
            (la, lb) = (np.log(a), np.log(b))
            for _ in range(n_bisect):
                m = 0.5 * (la + lb)
                above = risk2(np.exp(m)) > tgt2
                la = np.where(above, m, la)
                lb = np.where(above, lb, m)
            gm = np.exp(0.5 * (la + lb))
            gm = np.where(risk2(a) < tgt2, a, gm)
            g_prev = gm
        x[:, t] = _trade_diag(G[:, t], xp, gm, sig2, lam, lam0, lo, hi)[0]
        gam[:, t] = gm
        xp = x[:, t]
    return _multipliers_diag(desk, G, x, gam, lo, hi)

_multipliers_diag recovers the diagonal spread subgradient and signed box multiplier from the optimality equation (section 3.3; Durin Eqs. 17, 21). These are closed-form multipliers, not outputs of CVXPY.

def _multipliers_diag(desk, G, x, gam, lo, hi) -> Solution:
    """Recover s and nu from the KKT equation given positions and gamma."""
    sig2 = desk.sigma2()
    (lam, lam0) = (desk.lam(), desk.lam0())
    xp = np.concatenate([np.zeros_like(x[:, :1]), x[:, :-1]], axis=1)
    dx = x - xp
    g = gam[:, :, None]
    q = (lam + g) * sig2
    z = G - g * sig2 * xp
    dstar = np.sign(z) * np.maximum(np.abs(z) - lam0, 0.0) / q
    xu = xp + dstar
    at_bound = (xu > hi[:, None]) | (xu < lo[:, None])
    if lam0 > 0:
        s = np.where(dx != 0, np.sign(dx), np.clip(z / lam0, -1.0, 1.0))
    else:
        s = np.zeros_like(dx)
    nu = np.where(at_bound, z - q * dx - lam0 * s, 0.0)
    return Solution(x=x, dx=dx, gamma=gam, s=s, nu=nu, at_bound=at_bound)

solve_cvxpy solves the same static problem with a dense covariance (Durin Eqs. 5, 7, 15, 19). In risk-target mode the risk constraint’s dual is $\mu$ , so $\gamma=2\mu$ . In fixed- $\gamma$ mode it enters the objective directly.

Clarabel tolerances match the study. Trades and positions within $10^{-6}\sigma^*/\sqrt N$ of zero or a bound are snapped, as in the study; the raw primal values remain available. The solver’s box duals remain untouched. Nonzero trades use their actual sign; at zero trades the KKT residual selects a spread subgradient in $[-1,1]$ . Numerical KKT and attribution discrepancies are measured below, rather than removed by reconstructing box multipliers. This compact solver requires finite box bounds.

If Clarabel fails while reusing solver state, retry the identical problem once from a cold start, with the same tolerances. Retry events and solver warnings are retained and printed; a second failure stops execution. Warning messages are recorded without local file paths.

Finite bounds also provide a rigorous risk bound. Put $m_i=\max(|lo_i|,|hi_i|)$ . Every feasible position satisfies $x^\top\Sigma x\leq m^\top|\Sigma|m$ . If the requested squared risk target exceeds this bound, its constraint is redundant. We solve the same problem with $\gamma=0$ and omit the redundant inequality, whose dual may be chosen as zero. This exact constraint elimination avoids a poorly scaled, enormous right-hand side; it does not cap targets or alter their updates.

def solve_cvxpy(desk, G, sigma_target=None, gamma_path=None, lo=None, hi=None,
                snap_tolerance=1e-6):
    """Dense static optimisation; retain raw primal values and box duals."""
    B, T, N = G.shape
    Sigma = desk.covariance()
    lam, spread = desk.lam(), desk.lam0()
    if lo is None or hi is None:
        lo, hi = desk.bounds()
    lo, hi = _bcast_bounds(lo, hi, B, N)
    if not (np.isfinite(lo).all() and np.isfinite(hi).all()):
        raise ValueError("This compact CVXPY example requires finite position bounds.")
    dxv, xv = cp.Variable(N), cp.Variable(N)
    linear, previous = cp.Parameter(N), cp.Parameter(N)
    qscale, risk_squared = cp.Parameter(nonneg=True), cp.Parameter(nonneg=True)
    lower, upper = cp.Parameter(N), cp.Parameter(N)
    upper_constraint, lower_constraint = xv <= upper, xv >= lower
    constraints = [xv - dxv == previous, upper_constraint, lower_constraint]
    risk_constraint = cp.quad_form(xv, Sigma) <= risk_squared
    objective = cp.Maximize(linear @ dxv - 0.5 * qscale * cp.quad_form(dxv, Sigma)
                            - spread * cp.norm1(dxv))
    risk_problem = cp.Problem(objective, constraints + [risk_constraint])
    gamma_problem = cp.Problem(objective, constraints)
    shape = (B, T, N)
    x, raw, trades = np.zeros(shape), np.zeros(shape), np.zeros(shape)
    upper_dual, lower_dual = np.zeros(shape), np.zeros(shape)
    gamma = np.zeros((B, T))
    at_bound = np.zeros(shape, dtype=bool)
    snap = snap_tolerance * desk.risk_target / np.sqrt(N)
    for b in range(B):
        maximum_position = np.maximum(np.abs(lo[b]), np.abs(hi[b]))
        variance_bound = maximum_position @ np.abs(Sigma) @ maximum_position
        risk_required = gamma_path is None and float(sigma_target[b])**2 < variance_bound
        xp = np.zeros(N)
        for t in range(T):
            previous.value, lower.value, upper.value = xp, lo[b], hi[b]
            if risk_required:
                linear.value, qscale.value = G[b, t], lam
                risk_squared.value = float(sigma_target[b])**2
                problem = risk_problem
            else:
                gm = 0.0 if gamma_path is None else float(gamma_path[b, t])
                linear.value = G[b, t] - gm * (Sigma @ xp)
                qscale.value = lam + gm
                problem = gamma_problem
            options = dict(solver="CLARABEL", tol_gap_abs=1e-10,
                           tol_gap_rel=1e-10, tol_feas=1e-10, tol_ktratio=1e-8)
            with warnings.catch_warnings(record=True) as caught:
                warnings.simplefilter("always")
                try:
                    problem.solve(**options)
                except cp.error.SolverError:
                    solver_events.append(dict(event="cold retry", run=b, period=t,
                                              risk_target=(None if gamma_path is not None
                                                           else float(sigma_target[b]))))
                    try:
                        problem.solve(**options, warm_start=False)
                    except cp.error.SolverError as error:
                        target = None if gamma_path is not None else float(sigma_target[b])
                        raise RuntimeError(
                            f"Clarabel failed after cold retry: run={b}, period={t}, "
                            f"risk_target={target}; {error}") from None
            for warning in caught:
                solver_events.append(dict(event="solver warning", message=str(warning.message)))
            if problem.status not in ("optimal", "optimal_inaccurate"):
                raise RuntimeError(f"Clarabel status: {problem.status}")
            if risk_required:
                gm = 2 * float(np.asarray(risk_constraint.dual_value).item())
            raw[b, t] = xv.value
            delta = np.where(np.abs(dxv.value) < snap, 0.0, dxv.value)
            xt = xp + delta
            on_hi, on_lo = np.abs(xt - hi[b]) < snap, np.abs(xt - lo[b]) < snap
            xt = np.where(on_hi, hi[b], np.where(on_lo, lo[b], xt))
            x[b, t], trades[b, t], gamma[b, t] = xt, xt - xp, gm
            upper_dual[b, t] = upper_constraint.dual_value
            lower_dual[b, t] = lower_constraint.dual_value
            at_bound[b, t] = on_hi | on_lo
            xp = xt
    nu = upper_dual - lower_dual
    # Use actual trade signs. At zero trades, choose a subgradient in [-1, 1]
    # from the KKT residual, leaving the solver's box duals unchanged.
    w = G - lam * trades @ Sigma - gamma[:, :, None] * (x @ Sigma) - nu
    s = (np.where(trades != 0, np.sign(trades), np.clip(w / spread, -1, 1))
         if spread > 0 else np.zeros_like(trades))
    return Solution(x, trades, gamma, s, nu, at_bound,
                    upper_dual, lower_dual, raw)

solve selects the diagonal closed form or CVXPY. All rerun methods below use this same entry point.

def solve(desk, G, **options):
    return (solve_diag if desk.cov == "diag" else solve_cvxpy)(desk, G, **options)

Diagonal oracle

This repeats the study’s diagonal-oracle configuration: 8 assets, 30 periods, box 1, spread 0.8, seed 3. The declared study tolerances are positions $10^{-5}$ absolute, box multipliers $10^{-4}$ absolute, and active risk aversion $10^{-3}$ relative. CVXPY’s inactive risk dual tends to zero; the closed form retains its floor, so inactive risk aversion is reported separately.

The first configuration has no active risk periods. A second case keeps the same seed and dimensions but uses box 1.5 to exercise active risk duals. At a zero trade on a position bound, spread and box multipliers are not separately identified: the raw solver box dual can differ from the closed form while their combined force agrees. We print the full raw discrepancy, compare raw duals where the trade is nonzero, and separately compare the study’s explicitly reconstructed spread-first convention. The raw duals used elsewhere remain unchanged.

Risk aversion is also non-unique when the portfolio remains on the risk sphere without a resolved trade: the spread subgradient can offset a range of risk multipliers. We print that disagreement too. The relative gamma assertion uses active periods with a trade larger than the solver’s declared snapping scale on an unbounded coordinate with nonzero position, and requires this mask to be nonempty in the box-1.5 case. This is an identifiability restriction, not a larger tolerance.

for box in (1.0, 1.5):
    oracle = Desk(n_assets=8, n_periods=30, box=box, spread_rel=0.8)
    g, r = simulate(oracle, 3)
    G = g.sum(0)[None]
    diagonal = solve_diag(oracle, G, sigma_target=np.array([1.0]))
    convex = solve_cvxpy(oracle, G, sigma_target=np.array([1.0]))
    active = diagonal.gamma[0] > 1e-6 * oracle.gamma_ref()
    resolved_trade = np.abs(diagonal.dx) > 1e-6 * oracle.risk_target / np.sqrt(oracle.n_assets)
    identifiable = resolved_trade
    risk_identifiable = active & np.any(
        resolved_trade[0] & ~diagonal.at_bound[0]
        & (np.abs(diagonal.x[0]) > 1e-6 * oracle.risk_target / np.sqrt(oracle.n_assets)), axis=1)
    assert identifiable.any()
    np.testing.assert_allclose(convex.x, diagonal.x, atol=1e-5)
    np.testing.assert_allclose(convex.nu[identifiable], diagonal.nu[identifiable], atol=1e-4)
    np.testing.assert_allclose(convex.gamma[0, risk_identifiable], diagonal.gamma[0, risk_identifiable], rtol=1e-3)
    if box > 1:
        assert risk_identifiable.any(), "Risk-dual comparison must be non-vacuous"
    # Explicitly reproduce the study's spread-first convention for comparison.
    # This diagnostic does not overwrite convex.nu or its raw dual arrays.
    w = G - oracle.lam() * convex.dx - convex.gamma[:, :, None] * convex.x
    convention_s = np.where(convex.dx != 0, np.sign(convex.dx),
                            np.clip(w / oracle.lam0(), -1, 1))
    convention_nu = np.where(convex.at_bound, w - oracle.lam0() * convention_s, 0)
    np.testing.assert_allclose(convention_nu, diagonal.nu, atol=1e-4)
    fixed = solve_cvxpy(oracle, G, gamma_path=diagonal.gamma)
    np.testing.assert_allclose(fixed.x, diagonal.x, atol=1e-5)
    print(f"Oracle box={box}; active risk periods={active.sum()}")
    print("Unidentified box coordinates:", int((~identifiable).sum()))
    print("Active risk periods with identifiable gamma:", int(risk_identifiable.sum()))
    if active.any():
        print("Raw active gamma max absolute difference (including periods not identified by this test):",
              np.max(np.abs(convex.gamma[0, active] - diagonal.gamma[0, active])))
    print("Maximum position error:", np.max(np.abs(convex.x - diagonal.x)))
    print("Raw box-dual error, all coordinates:", np.max(np.abs(convex.nu - diagonal.nu)))
    print("Raw box-dual error, identifiable coordinates:",
          np.max(np.abs(convex.nu[identifiable] - diagonal.nu[identifiable])))
    print("Study-convention box-dual error:", np.max(np.abs(convention_nu - diagonal.nu)))
    print("Maximum fixed-gamma position error:", np.max(np.abs(fixed.x - diagonal.x)))
    print("Maximum raw-to-snapped position change:", np.max(np.abs(convex.primal_raw - convex.x)))
    if risk_identifiable.any():
        print("Identifiable active gamma max relative error:", np.max(np.abs(
            convex.gamma[0, risk_identifiable] / diagonal.gamma[0, risk_identifiable] - 1)))
    if (~active).any():
        print("Inactive gamma, CVXPY max / diagonal floor:",
              convex.gamma[0, ~active].max(), 1e-8 * oracle.gamma_ref())
print("PASS: positions, identifiable active risk duals, identifiable raw box duals and study-convention box duals")
print("Solver events:", solver_events)
solver_events.clear()

Oracle box=1.0; active risk periods=0
Unidentified box coordinates: 208
Active risk periods with identifiable gamma: 0
Maximum position error: 8.693970349193947e-09
Raw box-dual error, all coordinates: 0.0938566251796392
Raw box-dual error, identifiable coordinates: 4.006360382335927e-09
Study-convention box-dual error: 4.005961229402999e-09
Maximum fixed-gamma position error: 1.9367057402241983e-09
Maximum raw-to-snapped position change: 8.329894583113052e-09
Inactive gamma, CVXPY max / diagonal floor: 0.0 3.2756529752049163e-09
Oracle box=1.5; active risk periods=26
Unidentified box coordinates: 152
Active risk periods with identifiable gamma: 22
Raw active gamma max absolute difference (including periods not identified by this test): 0.07429351985869384
Maximum position error: 6.956109848826886e-07
Raw box-dual error, all coordinates: 0.10709899976812781
Raw box-dual error, identifiable coordinates: 8.28780995165701e-07
Study-convention box-dual error: 6.768218758224087e-07
Maximum fixed-gamma position error: 7.985581536917863e-06
Maximum raw-to-snapped position change: 1.2514183911349619e-08
Identifiable active gamma max relative error: 7.426598902293691e-06
Inactive gamma, CVXPY max / diagonal floor: 1.0608439388721236e-08 3.2756529752049163e-09
PASS: positions, identifiable active risk duals, identifiable raw box duals and study-convention box duals
Solver events: [{'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}]

Durin

2. Durin attribution

effective_position_matrix builds the general rank-one bound terms from vectors $v_c$ . Write a position constraint as $\epsilon_c v_c^\top x\leq\epsilon_c M_c$ , with signed bound $M_c$ , orientation $\epsilon_c\in\{-1,1\}$ and non-negative solver dual $\lambda^c$ . Then

$\eta_c=\frac{\epsilon_c\lambda^c}{2M_c},\qquad D_b=\sum_c2\eta_c v_cv_c^\top.$

This is section 3.3’s signed-bound form of Durin Eq. 37. For the coordinate box, $v_c$ is a coordinate vector and the oriented dual is nu. Zero bounds cannot use this division; their limiting condition is imposed in durin_matrix. The matrix helper accepts general constraint directions, while this notebook’s optimiser supplies coordinate bounds.

def effective_position_matrix(vectors, signed_bounds, oriented_duals):
    """Finite bound penalties: 2 eta_c v_c v_c.T; zero bounds handled separately."""
    eta = oriented_duals / (2 * signed_bounds)
    return vectors.T @ ((2 * eta)[:, None] * vectors)

durin_matrix forms the full matrices $\bar Q=(\lambda+\gamma)\Sigma+D_b+D_s$ and $\bar P=\gamma\Sigma+D_b$ (Durin Eqs. 30–31). On a nonzero joint trade, $\eta_s=\lambda_0/(2|\Delta x|)$ and $D_s=\operatorname{diag}(2\eta_s)$ (Eq. 48).

Each signal solves $\bar Q\Delta x^k=g^k-\bar P x^k_{t-1}$ . Zero joint trades impose $\Delta x^k=0$ (Eq. 54); a binding zero position imposes $x^k=0$ (Eq. 44). The remaining coordinates solve the reduced matrix equation, including their covariance with the fixed coordinates. This is a matrix solve, not an asset-by-asset diagonal formula.

Spread penalties act on trades and enter only $\bar Q$ ; position penalties act on the new position and enter both $\bar Q$ and $\bar P$ .

The zero-trade limit applies only with positive spread. With zero spread, opposite component trades may cancel; the ordinary matrix recursion must retain them.

def durin_matrix(desk, sol, g, b=0):
    K, T, N = g.shape
    Sigma, lam, spread = desk.covariance(), desk.lam(), desk.lam0()
    xk, xp = np.zeros((K, T, N)), np.zeros((K, N))
    coordinate_vectors = np.eye(N)
    for t in range(T):
        x, dx, nu = sol.x[b, t], sol.dx[b, t], sol.nu[b, t]
        bound, gamma = sol.at_bound[b, t], sol.gamma[b, t]
        zero_trade = (dx == 0) & (spread > 0)
        zero_position = bound & (x == 0) & ~zero_trade
        free = ~(zero_trade | zero_position)
        finite_bound = free & bound
        position_penalty = effective_position_matrix(
            coordinate_vectors[finite_bound], x[finite_bound], nu[finite_bound])
        eta_spread = np.zeros(N)
        moving = free & (dx != 0)
        eta_spread[moving] = spread / (2 * np.abs(dx[moving]))
        Qbar = (lam + gamma) * Sigma + position_penalty + np.diag(2 * eta_spread)
        Pbar = gamma * Sigma + position_penalty
        delta = np.where(zero_position, -xp, 0.0)
        indices = np.flatnonzero(free)
        rhs = g[:, t, indices] - xp @ Pbar[indices].T - delta @ Qbar[indices].T
        delta[:, indices] = np.linalg.solve(Qbar[np.ix_(indices, indices)], rhs.T).T
        xk[:, t] = xp + delta
        xp = xk[:, t]
    return xk, np.diff(xk, axis=1, prepend=0.0)

durin_positions retains the diagonal specialisation for the published-path check and calls the full matrix recursion for dense covariance. Its spread and bound penalties are the same effective multipliers (section 3.3; Durin Eqs. 37, 48).

def durin_positions(desk: Desk, sol: Solution, g: np.ndarray, b: int=0):
    """Signal-wise positions x^k and trades dx^k of run ``b`` of ``sol``."""
    if desk.cov != "diag":
        return durin_matrix(desk, sol, g, b)
    (K, T, N) = g.shape
    sig2 = desk.sigma2()
    (lam, lam0) = (desk.lam(), desk.lam0())
    xk = np.zeros((K, T, N))
    xp = np.zeros((K, N))
    for t in range(T):
        (x, dx, nu, ab, gm) = (sol.x[b, t], sol.dx[b, t], sol.nu[b, t], sol.at_bound[b, t], sol.gamma[b, t])
        zero_trade = (dx == 0) & (lam0 > 0)
        zero_pos = ab & (x == 0) & ~zero_trade
        free = ~zero_trade & ~zero_pos
        s = sol.s[b, t]
        with np.errstate(divide='ignore', invalid='ignore'):
            Ds = np.divide(lam0 * s, dx, out=np.zeros_like(dx), where=free & (dx != 0))
            Db = np.where(free & ab, nu / x, 0.0)
        Qbar = (lam + gm) * sig2 + Ds + Db
        Pbar = gm * sig2 + Db
        dk = (g[:, t] - Pbar * xp) / Qbar
        dk = np.where(zero_trade, 0.0, np.where(zero_pos, -xp, dk))
        xk[:, t] = xp + dk
        xp = xk[:, t]
    dxk = np.diff(xk, axis=1, prepend=0.0)
    return (xk, dxk)

allocations returns gross P&L, quadratic cost, spread cost, net P&L and, when given the total position, the contribution $x_k^\top\Sigma x$ to portfolio variance. Costs use the joint trade: $\tfrac12\Delta x_k^\top\Lambda\Delta x$ and $\lambda_0\Delta x_k^\top\operatorname{sgn}(\Delta x)$ (Durin Eqs. 13–14, 22). These signed contributions can be negative; their sums are checked separately.

def allocations(desk, xk, dxk, dx_total, returns, x_total=None):
    Sigma = desk.covariance()
    gross = (xk * returns).sum(-1)
    quadratic_cost = 0.5 * desk.lam() * (dxk * (dx_total @ Sigma)).sum(-1)
    spread_cost = desk.lam0() * (dxk * np.sign(dx_total)).sum(-1)
    result = dict(gross=gross, quadratic_cost=quadratic_cost, spread_cost=spread_cost,
                  net=gross - quadratic_cost - spread_cost)
    if x_total is not None:
        result["variance"] = (xk * (x_total @ Sigma)).sum(-1)
    return result

pnl_split applies that cost convention to component positions. The diagonal arithmetic is preserved to reproduce the published endpoints exactly.

def pnl_split(desk, xk, dxk, dx_total, r):
    """P&L of component positions xk with costs split as Durin Eqs. 14 and 22."""
    if desk.cov != "diag":
        return allocations(desk, xk, dxk, dx_total, r)["net"]
    (lam, lam0) = (desk.lam(), desk.lam0())
    sig2 = desk.sigma2()
    gross = (xk * r).sum(-1)
    quad = 0.5 * lam * (dxk * sig2 * dx_total).sum(-1)
    spread = lam0 * (dxk * np.sign(dx_total)).sum(-1)
    return gross - quad - spread

pnl evaluates the full run with the same position, quadratic-cost and spread-cost ledger (Durin Eq. 5 and section 3.3).

def pnl(desk: Desk, sol: Solution, r: np.ndarray) -> np.ndarray:
    """Net P&L per period (B, T): x_t . r_{t+1} minus quadratic and spread costs."""
    return pnl_split(desk, sol.x, sol.dx, sol.dx, r)

component_pnl evaluates signal or constraint portfolios using the joint run’s trades for the cost allocation (Durin Eqs. 14, 22).

def component_pnl(desk, xk, dxk, sol: Solution, r, b: int=0):
    """P&L per component with Durin's cost split (Eqs. 14, 22)."""
    return pnl_split(desk, xk, dxk, sol.dx[b][None], r[None])

Other methods

3. Other methods

grinold_easton_positions sends each signal, the spread force $-\lambda_0s$ and the box force $-\nu$ through $Q\Delta x+Px_{t-1}=\text{source}$ (section 3.3; Durin Eqs. 10–12, 17, 21). It keeps separate spread and position-limit portfolios and uses full matrices for dense covariance.

def grinold_easton_positions(desk, sol, g):
    K, T, N = g.shape
    lam, spread = desk.lam(), desk.lam0()
    sources = np.concatenate([g, -spread * sol.s[0][None], -sol.nu[0][None]])
    xk, xp = np.zeros((K + 2, T, N)), np.zeros((K + 2, N))
    Sigma = desk.covariance()
    for t in range(T):
        gamma = sol.gamma[0, t]
        if desk.cov == "diag":
            sig2 = desk.sigma2()
            delta = (sources[:, t] - gamma * sig2 * xp) / ((lam + gamma) * sig2)
        else:
            Q, P = (lam + gamma) * Sigma, gamma * Sigma
            delta = np.linalg.solve(Q, (sources[:, t] - xp @ P).T).T
        xk[:, t] = xp + delta
        xp = xk[:, t]
    return xk, np.diff(xk, axis=1, prepend=0.0)

regress fits non-negative least squares without an intercept to stand-alone net P&Ls (section 3.1; Durin’s discussion on pp. 2–3). The unexplained book P&L remains a separate residual.

def regress(total, proxies):
    beta = nnls(proxies.T, total)[0]
    attributed = beta[:, None] * proxies
    return beta, attributed, total - attributed.sum(0)

iterated_proxies reruns each signal at its updated risk target (section 3.1). Targets are multiplied by their NNLS coefficients without damping. Live zero coefficients drop proxies; zero-P&L proxies leave the convergence criterion. No live proxy gives no_live_proxy, not convergence. Convergence requires live coefficients within $10^{-3}$ of one, with at most 50 fits.

def iterated_proxies(desk, g, total, r, tol=1e-3, max_iter=50):
    K, T, N = g.shape
    sigma = np.full(K, desk.risk_target)
    status = "max_iter"
    for iteration in range(1, max_iter + 1):
        active = np.flatnonzero(sigma > 0)
        proxies = np.zeros((K, T))
        if len(active):
            sol = solve(desk, g[active], sigma_target=sigma[active])
            proxies[active] = pnl(desk, sol, r)
        active = active[np.any(proxies[active] != 0, axis=1)]
        beta = np.zeros(K)
        if len(active):
            beta[active], _, _ = regress(total, proxies[active])
        if not len(active):
            status = "no_live_proxy"
            break
        if np.max(np.abs(beta[active] - 1)) < tol:
            status = "converged"
            break
        if iteration == max_iter:
            break
        sigma[active] *= beta[active]
    return beta[:, None] * proxies, dict(
        status=status, n_iter=iteration, beta=beta, risk_targets=sigma)

coalition_signals sums forecasts for every non-empty subset, indexed by bitmask (section 3.4’s game). Every subset uses the same construction rules.

def coalition_signals(g: np.ndarray) -> np.ndarray:
    """Combined signal of every coalition, indexed by bitmask 1..2^K-1."""
    K = g.shape[0]
    masks = np.arange(1, 2 ** K)
    member = masks[:, None] >> np.arange(K) & 1
    return np.einsum('mk,ktn->mtn', member.astype(float), g)

shapley_weights computes the probability weight of each coalition size in a random ordering of signals (section 3.4).

def shapley_weights(K: int):
    return np.array([factorial(s) * factorial(K - s - 1) / factorial(K) for s in range(K)])

shapley averages marginal P&L over all orders. The empty coalition’s value is zero; this calculation is applied period by period (section 3.4).

def shapley(v: np.ndarray, K: int) -> np.ndarray:
    """Exact Shapley values per period."""
    w = shapley_weights(K)
    masks = np.arange(2 ** K)
    size = np.array([bin(m).count('1') for m in masks])
    phi = np.zeros((K, v.shape[1]))
    for k in range(K):
        S = masks[masks >> k & 1 == 0]
        phi[k] = (w[size[S]][:, None] * (v[S | 1 << k] - v[S])).sum(0)
    return phi

singleton_index identifies the stand-alone coalitions, so the first NNLS fit and the risk-target coalition game use exactly the same paths (sections 3.1 and 3.4).

def singleton_index(K: int) -> np.ndarray:
    """Bitmask of each singleton coalition."""
    return 1 << np.arange(K)

compare calculates all seven displayed methods. It evaluates both coalition risk conventions, computes leave-one-out as the full-coalition value minus each complementary coalition, and keeps proxy and leave-one-out residuals separate. The full-portfolio solver and the rerun solver are the same.

def compare(desk, seed=20261029, include_iterated=True):
    g, r = simulate(desk, seed)
    K, T, _ = g.shape
    full = solve(desk, g.sum(0)[None], sigma_target=np.array([desk.risk_target]))
    total = pnl(desk, full, r)[0]
    parts = {}
    for name, positions in [("durin", durin_positions),
                            ("grinold_easton", grinold_easton_positions)]:
        x, dx = positions(desk, full, g)
        parts[name] = component_pnl(desk, x, dx, full, r)
    coalitions = coalition_signals(g)
    for convention in ("target", "gamma"):
        options = (dict(sigma_target=np.full(len(coalitions), desk.risk_target))
                   if convention == "target" else
                   dict(gamma_path=np.broadcast_to(full.gamma, (len(coalitions), T))))
        sol = solve(desk, coalitions, **options)
        v = np.vstack([np.zeros(T), pnl(desk, sol, r)])
        parts["shapley_" + convention] = shapley(v, K)
        if convention == "target":
            _, parts["regression"], _ = regress(total, v[singleton_index(K)])
            parts["loo"] = v[-1] - v[(2**K - 1) ^ singleton_index(K)]
    iteration = None
    if include_iterated:
        parts["emmanuel"], iteration = iterated_proxies(desk, g, total, r)
    for method in ("regression", "loo", "emmanuel"):
        if method not in parts:
            continue
        parts[method] = np.vstack([parts[method], total - parts[method].sum(0)])
    return total, parts, iteration

Reproduction

4. Reproduction and checks

First run a small one-factor example (8 assets, 30 periods, seed 0), using CVXPY only. Long-only bounds exercise the zero-position limit as well as the no-trade limit. Raw duals are used, so the declared absolute tolerance is $10^{-5}$ for dense attribution sums, consistent with the diagonal oracle’s position tolerance. This is a numerical reconciliation check, not the exact endpoint check below. The output gives the errors in positions, trades, gross P&L, each cost, net P&L and variance, plus the KKT residual.

The separate four-asset/four-period example also runs the methods through compare. Its iterated targets can exceed the maximum risk permitted by the box. The solver removes such redundant risk inequalities using the rigorous bound above; a separate check compares this branch with an explicit zero-risk-aversion run. Iteration status is printed as returned, including max_iter when applicable. Default compare raises any remaining solver failure.

dense = Desk(n_assets=8, n_periods=30, cov="factor", long_only=True,
             box=1.0, spread_rel=0.8)
g, r = simulate(dense, 0)
full = solve_cvxpy(dense, g.sum(0)[None], sigma_target=np.array([1.0]))
Sigma = dense.covariance()
kkt = (g.sum(0) - dense.lam() * full.dx[0] @ Sigma
       - full.gamma[0, :, None] * (full.x[0] @ Sigma)
       - dense.lam0() * full.s[0] - full.nu[0])
print("Dense KKT max residual:", np.max(np.abs(kkt)))
print("Raw box-dual identity error:",
      np.max(np.abs(full.nu - (full.dual_upper - full.dual_lower))))
zero_trades = full.dx[0] == 0
zero_bounds = full.at_bound[0] & (full.x[0] == 0)
previous = np.concatenate([np.zeros_like(full.x[0, :1]), full.x[0, :-1]])
liquidations = zero_bounds & (previous > 0) & ~zero_trades
assert zero_trades.any() and zero_bounds.any() and liquidations.any()
print("Actual liquidations from positive position to zero:", liquidations.sum())
print("Zero joint trades / active zero bounds:", zero_trades.sum(), zero_bounds.sum())
ledger = allocations(dense, full.x[0], full.dx[0], full.dx[0], r, full.x[0])
for name, positions in [("Durin", durin_positions), ("Grinold–Easton", grinold_easton_positions)]:
    xk, dxk = positions(dense, full, g)
    np.testing.assert_allclose(component_pnl(dense, xk, dxk, full, r),
                               allocations(dense, xk, dxk, full.dx[0], r)["net"],
                               atol=1e-12, rtol=0)
    a = allocations(dense, xk, dxk, full.dx[0], r, full.x[0])
    errors = dict(positions=np.max(np.abs(xk.sum(0) - full.x[0])),
                  trades=np.max(np.abs(dxk.sum(0) - full.dx[0])))
    errors.update({key: np.max(np.abs(values.sum(0) - ledger[key]))
                   for key, values in a.items()})
    if name == "Durin":
        np.testing.assert_array_equal(dxk[:, zero_trades], 0)
        np.testing.assert_array_equal(xk[:, zero_bounds], 0)
    for key, error in errors.items():
        assert error <= 1e-5, (name, key, error)
    print(name, {key: float(value) for key, value in errors.items()})
print("PASS: dense allocations, variance sums and zero-trade/zero-bound limits")

# Zero-spread cancellation: zero joint trade must allow opposite component trades.
for covariance in ("diag", "factor"):
    cancellation = Desk(n_assets=2, n_periods=2, cov=covariance, spread_rel=0.0)
    opposing = np.array([[[0.05, -0.03], [0.04, -0.02]],
                         [[-0.05, 0.03], [-0.04, 0.02]]])
    joint = solve(cancellation, opposing.sum(0)[None], gamma_path=np.ones((1, 2)))
    xk, dxk = durin_positions(cancellation, joint, opposing)
    assert np.any(np.abs(dxk) > 1e-3)
    np.testing.assert_allclose(dxk.sum(0), joint.dx[0], atol=1e-8, rtol=0)
    np.testing.assert_allclose(xk.sum(0), joint.x[0], atol=1e-8, rtol=0)
print("PASS: zero-spread cancelling trades, diagonal and dense")

# Check the redundant-risk branch against the same problem at gamma=0.
small_dense = Desk(n_assets=4, n_periods=4, cov="factor")
g_small, r_small = simulate(small_dense, 0)
large_target = solve_cvxpy(small_dense, g_small.sum(0)[None],
                          sigma_target=np.array([11881.909026551642]))
zero_gamma = solve_cvxpy(small_dense, g_small.sum(0)[None],
                        gamma_path=np.zeros((1, small_dense.n_periods)))
np.testing.assert_array_equal(large_target.gamma, 0)
np.testing.assert_allclose(large_target.x, zero_gamma.x, atol=1e-8, rtol=0)
np.testing.assert_allclose(large_target.nu, zero_gamma.nu, atol=1e-8, rtol=0)
print("PASS: redundant risk target equals gamma=0; max position error",
      float(np.max(np.abs(large_target.x - zero_gamma.x))))
total_dense, methods_dense, _ = compare(small_dense, seed=0, include_iterated=False)
for method, values in methods_dense.items():
    assert np.isfinite(values).all()
    error = np.max(np.abs(values.sum(0) - total_dense))
    assert error <= 1e-5, (method, error)
    print("Dense method check:", method, "max period reconciliation error", float(error))
g_small, r_small = simulate(small_dense, 0)
try:
    attributed_dense, iteration_dense = iterated_proxies(small_dense, g_small, total_dense, r_small)
    residual_dense = total_dense - attributed_dense.sum(0)
    assert np.isfinite(attributed_dense).all()
    np.testing.assert_allclose(attributed_dense.sum(0) + residual_dense,
                               total_dense, atol=1e-5, rtol=0)
    print("Dense iteration status:", iteration_dense["status"], "fits:", iteration_dense["n_iter"])
except RuntimeError as error:
    print("DENSE ITERATION NOT VERIFIED:", str(error))
print("Solver events:", solver_events)
solver_events.clear()

Dense KKT max residual: 4.401681470917254e-09
Raw box-dual identity error: 0.0
Actual liquidations from positive position to zero: 8
Zero joint trades / active zero bounds: 203 118
Durin {'positions': 1.9557196417885336e-09, 'trades': 1.9557196417885336e-09, 'gross': 3.085415589687557e-09, 'quadratic_cost': 5.83463121567318e-11, 'spread_cost': 1.6965026605753053e-10, 'net': 3.085415589687557e-09, 'variance': 1.947450312123067e-09}
Grinold–Easton {'positions': 2.2032357271939418e-08, 'trades': 1.762616990905186e-08, 'gross': 2.7894036414011225e-08, 'quadratic_cost': 9.632457158292862e-11, 'spread_cost': 1.2902979899809264e-09, 'net': 2.7888603204573315e-08, 'variance': 6.229143978409013e-09}
PASS: dense allocations, variance sums and zero-trade/zero-bound limits
PASS: zero-spread cancelling trades, diagonal and dense
PASS: redundant risk target equals gamma=0; max position error 0.0

Dense method check: durin max period reconciliation error 3.3075678536331754e-07
Dense method check: grinold_easton max period reconciliation error 1.5779620150813134e-07
Dense method check: shapley_target max period reconciliation error 2.6728619317850644e-14
Dense method check: regression max period reconciliation error 1.3877787807814457e-17
Dense method check: loo max period reconciliation error 0.0
Dense method check: shapley_gamma max period reconciliation error 3.6393206892526564e-08

Dense iteration status: max_iter fits: 50
Solver events: [{'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}]

Check the two published paths

The reference below contains the final cumulative values copied from the figure’s f13-pnl-paths.json, including residuals and the two Grinold–Easton portfolios. It is embedded so this notebook works offline. The source JSON SHA-256 is 91b750176b9263260aa13a85063dbf2c1f3b12bb561139e752bf4f9cc8f53544.

The check uses absolute tolerance $10^{-8}$ simulation P&L units and zero relative tolerance. It compares every component endpoint for every method and the book endpoint. It also checks that each method, including any residual, adds back to the book in every period. Neither check changes the reference. The examples may take a few minutes.

reference = {'base': {'total': 519.0910938220723,
          'methods': {'regression': {'fast A': 100.53242064658039,
                                     'slow A': 70.38954858667016,
                                     'fast B': 14.304523929024574,
                                     'slow B': 5.152252732430301,
                                     'residual': 328.7123479273673},
                      'emmanuel': {'fast A': 152.764529210917,
                                   'slow A': 112.63119504250166,
                                   'fast B': 58.46371903282879,
                                   'slow B': 9.128433518895262,
                                   'residual': 186.10321701692942},
                      'loo': {'fast A': 111.60697161328982,
                              'slow A': 119.72681697759026,
                              'fast B': 45.54721780588971,
                              'slow B': 46.361535410303155,
                              'residual': 195.84855201499892},
                      'grinold_easton': {'fast A': 256.72848474852657,
                                         'slow A': 366.92552784896986,
                                         'fast B': 130.11725884186922,
                                         'slow B': 152.89402503899424,
                                         'spread cost': -241.89397508564036,
                                         'position limit': -145.68022757064776},
                      'durin': {'fast A': 175.06247412657459,
                                'slow A': 179.2643826576884,
                                'fast B': 101.85999505010044,
                                'slow B': 62.904241987708474},
                      'shapley_target': {'fast A': 177.48994085845976,
                                         'slow A': 176.86366093064652,
                                         'fast B': 98.66062755250604,
                                         'slow B': 66.07686448045925},
                      'shapley_gamma': {'fast A': 165.110443254377,
                                        'slow A': 201.34140232216058,
                                        'fast B': 84.36564359163518,
                                        'slow B': 68.27360465389944}}},
 'spread_2': {'total': 104.32833661086099,
              'methods': {'regression': {'fast A': 0.0,
                                         'slow A': 0.0,
                                         'fast B': 0.0,
                                         'slow B': 0.0,
                                         'residual': 104.32833661086099},
                          'emmanuel': {'fast A': 0.0,
                                       'slow A': 0.0,
                                       'fast B': 0.0,
                                       'slow B': 0.0,
                                       'residual': 104.32833661086099},
                          'loo': {'fast A': 76.94178438231884,
                                  'slow A': 79.71850185927482,
                                  'fast B': 83.2415101226865,
                                  'slow B': 38.198069494507166,
                                  'residual': -173.77152924792617},
                          'grinold_easton': {'fast A': 584.582758851705,
                                             'slow A': 1411.9782224503138,
                                             'fast B': 427.665383299356,
                                             'slow B': 414.98905688713205,
                                             'spread cost': -2717.716466794506,
                                             'position limit': -17.1706180831368},
                          'durin': {'fast A': 39.74723533758032,
                                    'slow A': 25.44001817571785,
                                    'fast B': 27.503588408741983,
                                    'slow B': 11.637494688821034},
                          'shapley_target': {'fast A': 30.38494084386784,
                                             'slow A': 30.230508340967734,
                                             'fast B': 29.720064717153566,
                                             'slow B': 13.992822708871934},
                          'shapley_gamma': {'fast A': 31.427851827138475,
                                            'slow A': 29.636409594471246,
                                            'fast B': 30.21533987280228,
                                            'slow B': 13.048735316449124}}}}

tolerance = 1e-8
for scenario, spread in [("base", 0.5), ("spread_2", 2.0)]:
    desk = Desk(spread_rel=spread)
    total, parts, iteration = compare(desk)
    expected = reference[scenario]
    book_endpoint = np.cumsum(total)[-1]
    np.testing.assert_allclose(book_endpoint, expected["total"], atol=tolerance, rtol=0)
    print(f"{scenario}: book={book_endpoint:.12f}; "
          f"iteration={iteration['status']} ({iteration['n_iter']} fits)")
    largest_error = abs(book_endpoint - expected["total"])
    for method, rows in expected["methods"].items():
        endpoints = np.cumsum(parts[method], axis=1)[:, -1]
        np.testing.assert_allclose(endpoints, list(rows.values()), atol=tolerance, rtol=0)
        np.testing.assert_allclose(parts[method].sum(0), total, atol=tolerance, rtol=0)
        largest_error = max(largest_error, np.max(np.abs(endpoints - list(rows.values()))))
        print(f"  {method}: " + ", ".join(
            f"{label}={value:.9f}" for label, value in zip(rows, endpoints)))
    print(f"PASS: all component endpoints and period sums; max endpoint error={largest_error:.3g}")

base: book=519.091093822072; iteration=converged (10 fits)
  regression: fast A=100.532420647, slow A=70.389548587, fast B=14.304523929, slow B=5.152252732, residual=328.712347927
  emmanuel: fast A=152.764529211, slow A=112.631195043, fast B=58.463719033, slow B=9.128433519, residual=186.103217017
  loo: fast A=111.606971613, slow A=119.726816978, fast B=45.547217806, slow B=46.361535410, residual=195.848552015
  grinold_easton: fast A=256.728484749, slow A=366.925527849, fast B=130.117258842, slow B=152.894025039, spread cost=-241.893975086, position limit=-145.680227571
  durin: fast A=175.062474127, slow A=179.264382658, fast B=101.859995050, slow B=62.904241988
  shapley_target: fast A=177.489940858, slow A=176.863660931, fast B=98.660627553, slow B=66.076864480
  shapley_gamma: fast A=165.110443254, slow A=201.341402322, fast B=84.365643592, slow B=68.273604654
PASS: all component endpoints and period sums; max endpoint error=0

spread_2: book=104.328336610861; iteration=no_live_proxy (2 fits)
  regression: fast A=0.000000000, slow A=0.000000000, fast B=0.000000000, slow B=0.000000000, residual=104.328336611
  emmanuel: fast A=0.000000000, slow A=0.000000000, fast B=0.000000000, slow B=0.000000000, residual=104.328336611
  loo: fast A=76.941784382, slow A=79.718501859, fast B=83.241510123, slow B=38.198069495, residual=-173.771529248
  grinold_easton: fast A=584.582758852, slow A=1411.978222450, fast B=427.665383299, slow B=414.989056887, spread cost=-2717.716466795, position limit=-17.170618083
  durin: fast A=39.747235338, slow A=25.440018176, fast B=27.503588409, slow B=11.637494689
  shapley_target: fast A=30.384940844, slow A=30.230508341, fast B=29.720064717, slow B=13.992822709
  shapley_gamma: fast A=31.427851827, slow A=29.636409594, fast B=30.215339873, slow B=13.048735316
PASS: all component endpoints and period sums; max endpoint error=0

↑ Sections

4. Correlated indicators: which value?

Venn diagram of an outcome and two overlapping indicators; the part of the overlap that is also about the outcome is hatched as shared information. Beside it, four games that can value a subset of indicators: refitted regression R², conditional variance, local conditional SHAP and SAGE, each with what its value measures. — Correlated indicators share information about the outcome; the chosen game sets what their shares mean.

4.1 Shared information is not a regression coefficient

An indicator can explain a P&L without having generated it. A variable that tracks another signal may predict the same losses even if the trading rule never reads it. Statistical attribution shares the information that indicators provide about an outcome. The coalition P&L game of the previous section shares the P&L obtained by rerunning a trading algorithm.

When indicators move together, a regression can predict their combined effect accurately while estimating their separate coefficients poorly. Small changes in the sample can move a large coefficient from one indicator to another, or reverse its sign. The coefficient describes a change holding the other inputs fixed, which may be a rare change in the observed data. Its size also depends on the input’s units.

In their game-theoretic analysis of regression, Lipovetsky & Conklin (2001) divide the regression’s explained variance. Each coalition is a subset of predictors, each subset gets its own fitted regression, and its value is its $R^2$ . Shapley averages the gain from adding an indicator over all orders of entry. A redundant indicator receives credit when it arrives before its substitute, even if it adds little when it arrives last.

Change the correlation between the paired indicators. The same synthetic model is fitted again across independent seeds; only the first pair is displayed, with true coefficients equal to 1 for both indicators. The other predictors stay in each regression.

Coefficients and predictive shares under collinearity

Synthetic correlated-indicator study. The coefficient scatter shows individual seeded fits; the other panels show means ±1 sample SD, not confidence intervals.

Within-pair correlation ρ

Indicator 1 (blue)Indicator 2 (orange)

When the inputs are nearly duplicates, the OLS coefficients vary more across samples. The Shapley shares of refitted $R^2$ vary less, because each averages marginal predictive gains over all predictor subsets. These are different estimands and their axes have different units. Stability does not establish causality or significance. (Lipovetsky & Conklin, 2001)

4.2 Choose the value before choosing the allocation

Let $N$ be the set of indicators, $d$ its size, and $v(S)$ the value of the subset $S$ . The allocation rule is

\begin{aligned} \phi_j(v)&=\sum_{S\subseteq N\setminus\{j\}}w(S)\,\Delta_j(S),\\ w(S)&=\frac{|S|!\,(d-|S|-1)!}{d!},\\ \Delta_j(S)&=v(S\cup\{j\})-v(S). \end{aligned}

The weights average marginal contributions over orders of arrival. The shares sum to $v(N)-v(\varnothing)$ . To interpret them, we still need to specify the units, baseline and meaning of $v$ .

For refitted regression $R^2$ , every subset gets a newly estimated model. Its value measures how much of the observed response variation that subset can fit. Training $R^2$ cannot fall when an ordinary least-squares model with an intercept gains a predictor. Its Shapley shares are therefore nonnegative, including when some predictors contain only noise.

For conditional variance, the value is the variance of the conditional mean prediction given the subset. The fitted model stays fixed. The value measures how much of its output variability those inputs can recover. In his paper on Sobol’ indices and Shapley values, Owen (2014) shows that, for independent inputs, the Shapley share of variance has an exact formula in terms of the variance decomposition (ANOVA) of the model, and lies between two Sobol’ indices: the variance an input explains alone and the variance it adds when it comes last. With correlated indicators, conditional variance is still a possible game, but this formula and these bounds no longer hold as such.

In their paper on interpreting model predictions, Lundberg & Lee (2017) introduce SHAP. For local conditional SHAP, the subset’s value is the conditional expected prediction at the particular observation being explained. The contributions are signed and measured in the prediction’s units. They add to that prediction minus its baseline, not to a fraction of variance. An indicator can push one prediction down and another up.

In their paper on global feature contributions, Covert et al. (2020) introduce SAGE. Its value is the reduction in expected prediction loss relative to using no features. A restricted model averages over missing inputs conditional on the retained ones. This makes the answer depend on both the fitted model and the loss: information useful for squared-error prediction need not have the same value for another objective.

In the synthetic linear experiment, three conditions hold: the linear model is the true conditional mean, the loss is squared error, and the predictors are jointly Gaussian, so the conditional mean given any subset is also linear. Under these conditions, refitted explained variance, conditional variance and SAGE give the same population shares once their normalisations match. Differences between their estimated global shares can come from training versus evaluation samples, estimation error or misspecification; this example does not establish distinct population rankings for those games.

What to notice

Even a stable share need not measure a causal effect, be statistically significant or represent a share of realised trading profit. State the target, the subset operation and the baseline beside the ranking.

4.3 Importance needs a comparison under the null

Simpler importance measures make different comparisons. A single-indicator fit can give high scores to several substitutes. Dropping an indicator and refitting may barely hurt because another replaces it. Permuting a column without refitting can create combinations absent from the joint distribution. In the synthetic linear example, near-duplicates have small drop-column importance, but permutation importance does not vanish for both copies. These quantities need not add to a common total.

Shuffling an indicator to create a null variable also breaks its correlations with the other indicators, so the shuffled copy is not a fair competitor. In their paper on false discovery control, Barber & Candès (2015) build a knockoff of each indicator: an artificial variable with the same correlations with the other indicators as the original, constructed without using the outcome. They then compare each original with its knockoff. Knockoff+ controls the false discovery rate for tests of each indicator’s coefficient, under their assumptions: a Gaussian linear model with the indicators treated as fixed (fixed design). It selects variables; it does not allocate profit.

That guarantee does not automatically cover serially dependent intraday errors. The synthetic knockoff comparison deliberately uses a separate design with independent errors. Its comparator that shuffles rows does not have the correlation structure knockoffs require, so it gives no FDR guarantee.

A model fitted to unrelated inputs and outcomes still produces shares, and one of them will be largest. So an attribution needs a null comparison even when we do not use it to select variables. In The Dead Salmons of AI Interpretability, Méloux & Dirupo (2025) motivate explicit null comparisons for interpretability; their experiments concern neural-network probing, not a validation or refutation of this Shapley estimator.

The synthetic attribution experiment uses a response generated independently of the indicators, so no indicator truly predicts it. It compares its shares with those from shuffled versions of that response, one observation at a time or in blocks. Blocks preserve some temporal dependence but cut it at their boundaries, so calibration is approximate. If the finding is that an indicator has the largest share, compare it with the largest share in each shuffled run. A band for one indicator chosen in advance does not account for choosing the top of the ranking. Shuffling the whole response tests whether any indicator predicts the outcome; it does not test whether one adds information when others truly predict it.

↑ Sections

5. The market maker’s ledger

A passive fill earns a spread relative to the current midprice. It also changes inventory, which remains exposed to the next price move. If the counterpart trades just before that move, spread capture can coexist with an inventory loss. We need to record both the spread earned and the inventory loss.

A buy fill on a time axis. The fill price is below the pre-fill mid, and the mid then falls through two markout horizons to the final mid. The gap from the fill price up to the pre-fill mid is the spread; the move from the pre-fill mid down to the final mid is positioning. A schematic waterfall then adds spread, positioning, hedge P&L, fees and rebates to give the P&L. — One fill split into spread and positioning, and the ledger lines that add up to the P&L (schematic).

In his study of high-frequency trading and the new market makers, Menkveld (2013) measures this for a trader in Dutch equities. He splits P&L into net spread, what fills earn against the mid after exchange and clearing fees, and positioning, gains or losses on inventory as the mid moves. The working-paper version reports €1.55 per trade of net spread and a positioning loss of €0.68 per trade. A frequency decomposition splits inventory into components that vary over different horizons. The component at horizons up to five seconds earns positive positioning profit; every longer horizon band loses. This is evidence about the trader and sample studied, not a law that inventory becomes unprofitable after five seconds.

5.1 Markouts and the reference price

For a book starting with zero inventory, let $d_i$ be the signed quantity of fill $i$ , positive for a purchase, $p_i$ its price, and $m_T$ the final midprice. Marked at $m_T$ , fill $i$ adds $d_i(m_T-p_i)$ to the P&L. Choose a reference price $r_i$ for each fill, for example the mid just before it, and split that amount in two:

d_i(m_T-p_i)=d_i(r_i-p_i)+d_i(m_T-r_i).

The first part is what the fill earned against the reference: the spread. The second is what the position earned as the price moved from the reference to the final mid: positioning. Summing over fills,

\begin{aligned} S(r)&=\sum_i d_i(r_i-p_i),\\ I(r)&=\sum_i d_i(m_T-r_i), \end{aligned}

and the ledger is

\mathrm{PnL}=S(r)+I(r)+H-\mathrm{fees}+\mathrm{rebates}.

Here $H$ is the gross marked P&L of a separately booked hedge instrument; fees and rebates are separate cash lines covering both quoting and hedging. The fill set below excludes those hedge trades, so they are not counted twice. If all instruments are instead consolidated into one fill ledger, the hedge line is already included and must not be added again.

The terminal inventory is marked at $m_T$ ; this does not include a cost of liquidating it. Changing $r_i$ cannot change the total. Moving the reference to a later mid adds $d_i(r_i^{\text{later}}-r_i^{\text{pre}})$ to spread and subtracts the same amount from positioning.

In their paper on dealing with inventory risk, Guéant et al. (2013) model quotes that respond to inventory risk and extend the model to adverse price jumps at fills. We can use this model to introduce a known adverse-selection effect and see where it appears in the accounts. With the pre-jump mid as reference, the jump enters inventory revaluation. With a later mark, part of the loss enters measured spread. A later mark also includes unrelated price moves and subsequent trading: the transfer is not a causal estimate of the current fill’s adverse selection.

For a fill at time $t_i$ , let $m(t_i+h)$ be the midprice at horizon $h$ after the fill. Its signed execution-to-mid markout is

d_i[m(t_i+h)-p_i].

This measures the price change from execution to the later mid, multiplied by the signed fill quantity. Some reports call only the mid-to-mid move a markout; here the execution price is the starting point. Take as reference the mid $h$ seconds after each fill, $r_i=m(t_i+h)$ , and write $r(h)$ for these references. Averaging the execution-to-mid amounts over $n$ fills gives the markout at horizon $h$ , $S(r(h))/n$ ; plotted against $h$ , it is the markout curve. $I(r(h))/n$ is the corresponding positioning line. For unequal sizes, also report the amount per unit traded, dividing by $\sum_i |d_i|$ . Always retain the fill count and traded quantity.

Compare the execution markouts at all four reference horizons. The second curve compares the change from the pre-fill reference with the adverse jump planted at the current fill.

Markouts and the market maker’s ledger

Synthetic market making: 32 fixed-seed paths of 600 seconds, quotes from the approximate formulas of Guéant, Lehalle and Fernandez-Tapia (GLFT), unit-size fills. Values are in ticks per fill, averaged over all fills of all paths. Positive horizons are capped at the terminal time; 0 s uses the pre-fill mid. Jump sizes range from 0 to 2 ticks per fill in steps of 0.5.

Planted adverse jump

Signed execution markout = side × (reference mid − fill price). At 0 s the reference is the pre-fill mid; later marks include the current jump, later fills and price moves. Positive means favourable to the market maker.

Blue: execution markout · Orange: change from pre-fill markout · Green dashed: planted current-fill jump cost (−ξ), not a prediction of a later mark.

Horizon	Markout	Positioning	Total

What to notice

Changing the reference transfers P&L between execution markout and positioning; their sum remains terminal marked wealth. Later marks contain the current jump, later fills’ jumps and Brownian moves. A positive positioning entry at a later reference does not establish that inventory management is profitable. The dashed line is the planted current-fill jump cost, not an expectation for the full later markout change.

When you change the jump, the simulation recomputes quotes and fills, so total P&L can change. Inventory starts at zero and is marked at the terminal mid. Hedge P&L, fees, rebates and liquidation costs are not modelled. This accounting illustration is not a calibrated trading result.

In his paper on implementation shortfall, Perold (1988) measures the shortfall by fixing a paper portfolio at the decision time and comparing it with what was actually implemented. The difference includes execution costs and the opportunity cost of intended trades that were never made. A market maker continuously revises quotes rather than executing a fixed investment list, so defining its paper portfolio is itself a research problem. The ledger above accounts for fills; it does not value missed opportunities.

5.2 From state indicators to a component explanation

With the ledger fixed, relate each fill’s markout or positioning component to state indicators available at the chosen decision time: spread, volatility, flow and inventory. Use the conditional predictive allocations of section 4, with the null comparisons and search corrections of section 6. A state observed after the fill may be an outcome of the trade, so its role must be justified before using it as a control.

The Azencott diagnostic described in section 0 scores threshold rules for bad-order classification. These scores are neither significance levels nor contributions in currency. Summing them cannot reconcile them with daily P&L. (Azencott et al., 2014)

In the synthetic experiment, the rarity rule adds a loss to orders whose volume rarity score exceeds 0.97, in six slices of the day. The score is the fraction of historical volume observations below the current volume. The detection rate counts how often volume alone passes the sensitivity and specificity floor in those planted slices. On matched days without the added loss, we count each variable-or-pair × slice cell that still passes: these are the null alarms.

For a planted loss of two spreads (select −2 spreads and the rarity rule below), a floor of 0.85 detects volume in only 2.8% of planted slices. Lowering the floor to 0.65 raises detection to 81.9%, but retains about 669 null cells per day. At 0.85, the simulated null retains none.

Lower the floor on sensitivity and specificity to keep more rules. The same change makes the planted volume effect easier to detect and produces more alarms on days without any planted effect.

Detection and false alarms

Synthetic sensitivity sweep over 24 matched day seeds per configuration. Detection is the fraction of six planted slices retaining the volume-only rule; null alarms count all retained variable-or-pair × slice cells over a 79-slice day. Error bars show ±1 across-day standard deviation, not confidence intervals.

Effect on slippage Affected orders per planted slice Sensitivity and specificity floor

━ Selected effect · ┄ Volume under null
● Selected floor · ☆ Original synthetic setting, paper floor

Volume singleton detection rate (%)

Loading synthetic sensitivity sweep…

What to notice

The figure opens at −2 spreads with the rarity rule. The marked 85% floor is the setting used by Azencott and coauthors. The −8-spread effect and the rarity-above-0.97 intervention are choices of this synthetic experiment, not parameters estimated by Azencott and coauthors. “Rarity rule” selects a varying number of orders; it does not mean zero affected orders. The other settings select an exact number of orders by their volume score, without using slippage.

Each retained variable-or-pair × slice cell counts as an alarm; it is not an independent discovery or an order-level error. The dotted curve shows how often the volume-only rule is retained under the null. Detection here is an in-sample association, not causal recovery or a forecast. Lowering the floor adds no false-discovery control; zero alarms in these simulated days do not establish a zero population rate. The lower two floors fall below the range explored by Azencott and coauthors.

Test the component explanation on later slices, keeping those results separate from the fit to the orders used to select it. The null comparison and search correction in the next section account for that selection.

↑ Sections

6. Statistical pitfalls

Two ways a statistical explanation can mislead: selecting the best of a search, and controlling for a collider.

6.1 Testing the best of many candidates

A research process can select indicators, horizons, thresholds and trading rules on the same history. Even if every candidate has zero expected excess performance, the winner tends to look profitable. Testing that winner as if we had chosen it before seeing the data ignores the other candidates we tried. The test needs to account for that search.

The familywise error rate, or FWER, is the probability of making any false rejection within a specified family. Bonferroni divides the allowed error probability by the number of tests; Holm uses a sequential correction that can reject more hypotheses while retaining FWER control. Both require valid individual p-values. Neither repairs a p-value computed under an incorrect model of temporal dependence.

The false discovery rate, or FDR, is the expected fraction of false rejections among all rejections, with that fraction defined as zero when nothing is selected. In their paper introducing FDR control, Benjamini & Hochberg (1995) give the Benjamini–Hochberg rule for independent tests. In their study of multiple testing under dependency, Benjamini & Yekutieli (2001) establish control under certain positive-dependence conditions and supply the more conservative Benjamini–Yekutieli correction for arbitrary dependence. Correlated indicators do not by themselves establish the positive-dependence condition required by BH.

FWER limits the chance of any false discovery in the family; FDR allows some false discoveries while controlling their expected share. Neither says that every selected indicator has a particular probability of being false. When every tested null hypothesis is true, every discovery is false, so FDR and FWER coincide.

In their study of the cross-section of expected returns, Harvey et al. (2016) apply these corrections to the growing set of tested return predictors. The t-statistic needed to claim significance depends on the search, so the same threshold need not suit every project.

Increase the number of indicators while keeping the same seeded sequence of null tests. The largest absolute statistic cannot decrease in this nested experiment. Every rejection here is a false discovery.

Searching more indicators raises the hurdle

Synthetic independent two-sided standard Normal tests at 5%. The z-statistics illustrate the large-sample t-statistic setting; no finite-sample t distribution is assumed.

Number of indicators tested

Bonferroni uses one cutoff for the whole family. Holm tests in order and stops at the first failure. Benjamini–Hochberg uses the last qualifying rank, controlling the expected false-discovery proportion under independence. The coloured comparison curves depend on rank; they cannot be used interchangeably as single cutoffs. (Benjamini & Hochberg, 1995) (Harvey et al., 2016)

This experiment uses independent Normal test statistics. It does not reproduce the bootstrap of White’s Reality Check or Romano and Wolf’s StepM (section 6.2), which needs the joint performance series of the strategies.

6.2 Test the strategy family

Trading strategies often share dates, positions and shocks. White’s Reality Check tests the joint null that no candidate beats a specified benchmark in expected performance. It compares the best observed performance with a bootstrap distribution of the maximum under the null. The resampling keeps the strategies together and, with an appropriate time-series bootstrap, retains temporal dependence.

Rejecting that joint null says that some candidate improves on the benchmark. In their paper on stepwise multiple testing, Romano & Wolf (2005) introduce StepM, which tests which candidates can be declared better than the benchmark, with asymptotic FWER control. It tests the strongest remaining statistics against a joint critical value, removes rejected nulls, and repeats on the remaining family. Dividing each candidate’s performance by its standard error (studentization) can help when candidates have different estimation uncertainties.

The synthetic comparison uses common bootstrap indices across strategies. Its finite-sample rejection frequencies illustrate the procedure; they do not prove that the asymptotic guarantee is exact in that sample. We also need a record of the candidates tried. Neither procedure corrects for searches over datasets or specifications left out of that record.

6.3 More controls can damage the explanation

In their paper on specification errors in factor investing, Lopez de Prado & Zoonekynd (2024) examine how controls can bias an estimate. A control variable can remove confounding, but a variable jointly affected by the proposed cause and the outcome is a collider. For example, $X$ could be a signal, $Y$ a markout and $Z$ a score used after the fact to select trades, which depends on both. Holding that score fixed restricts which signal and markout pairs we compare. In the model below, conditioning on $Z$ lowers the estimated effect of $X$ on $Y$ and can even reverse its sign, although the trading process has not changed. Choosing controls only because they improve predictions cannot tell these roles apart.

An execution outcome or a selected-trade flag may result from the process we want to explain. Check when and how it was generated before using it to explain P&L. Including it may improve the fit to observed losses while making an earlier indicator’s coefficient misleading.

The synthetic model starts with independent unit-variance Normal variables $X,u,v$ . The outcome is

Y=2X+u.

Here $X$ affects $Y$ , with noise $u$ . The control variable is

Z=\gamma Y+2X+v.

Here $\gamma$ controls how strongly $Y$ affects $Z$ , and $v$ is separate noise. This continuous synthetic model illustrates selection; it is not a fitted model of binary fills. The variable $Z$ is a consequence of both $X$ and $Y$ : the graph contains $X\to Z\leftarrow Y$ .

Conditioning on a collider

Exact population equations, with and without Z. Change γ while the causal coefficient of X stays fixed.

Collider coefficient γ (unitless)

Without $Z$ , the coefficient of $X$ is $2$ and the explained variance is

R^2=4/5.

With $Z$ included, the coefficient of $X$ becomes

(2-2\gamma)/(1+\gamma^2).

This is the fitted coefficient after controlling for $Z$ ; the causal effect of $X$ on $Y$ has not changed. The explained variance becomes

R^2=1-1/[5(1+\gamma^2)].

This measures the fit using both $X$ and $Z$ . Fit improves while the coefficient crosses zero. These equations hold in the population, so adding observations cannot remove the bias from controlling for the collider. (Lopez de Prado & Zoonekynd, 2024)

The parameter labels follow López de Prado’s equation (10): $\gamma$ multiplies $Y$ . The text describing Exhibit 6 (p. 14) swaps the two coefficients.

6.4 A Sharpe ratio is estimated, then often selected

A Sharpe ratio estimates a mean relative to volatility. Both quantities are uncertain, and the observation interval matters. In his paper on the statistics of Sharpe ratios, Lo (2002) shows that serial correlation changes Sharpe inference and time aggregation. Scaling a Sharpe ratio by the square root of the number of periods is not generally valid for dependent returns. More intraday observations do not provide the same information as equally many independent returns.

In their paper on the deflated Sharpe ratio, Bailey & López de Prado (2014) account for selection by comparing the selected estimate with a benchmark based on the expected maximum of zero-skill trials. The calculation incorporates sample length, skewness and kurtosis as well as the dispersion and effective number of trials. This model-based correction cannot recover trials that were never recorded. Correlated variants cannot be counted as independent attempts.

In their paper on the probability of backtest overfitting, Bailey et al. (2016) define PBO using the matrix of candidate returns. Combinatorially symmetric cross-validation repeatedly assigns blocks to an in-sample and complementary out-of-sample set, selects the in-sample winner, then records its relative out-of-sample rank. PBO measures how often that winner falls below the out-of-sample median of the candidates.

Falling in rank does not necessarily mean losing money, so PBO is not a probability of future loss. The splits overlap: they are neither independent replications nor an untouched future period. The synthetic search over monthly calendar rules (entry day, holding period, stop and side) lets us compare selection on noise with selection when a recurring effect is present. Its deflated-Sharpe illustration uses a separate hypothetical independent-trial model.

Choose an observed annual Sharpe ratio and increase the number of independent trials that preceded it. Even under zero skill, the expected maximum rises. The deflated Sharpe ratio compares the reported result with that selection benchmark. (Bailey & López de Prado, 2014)

Trial count changes the evidence

Hypothetical independent-trial example: default annual Sharpe 2.5, annualised trial-Sharpe variance 0.5, 1,250 daily observations, 250 days per year, skewness −3 and non-excess kurtosis 10.

Independent trials (count) Observed annual Sharpe ratio

The expected maximum uses the Bailey–López de Prado approximation; it is not exact in a finite sample. The deflated Sharpe ratio (DSR) uses per-period Sharpe units internally; the maximum is displayed in annual units. The model assumes independent trials. Counting correlated parameter combinations does not give that number.

The separate PBO panel uses the synthetic seasonal-rule experiment. Combinatorially symmetric cross-validation (CSCV) selects the best rule on half the blocks and ranks that same rule on the complementary blocks, over every symmetric split. PBO measures how often its rank falls below the median. It neither estimates the independent trial count for DSR nor gives a probability of future losses. (Bailey et al., 2016)

6.5 Losses do not identify their own cause

Long drawdowns can occur under an unchanged positive expected return. In their paper on when to worry about a drawdown, Rej et al. (2017) derive drawdown distributions for Brownian log-P&L with constant drift and volatility. Their duration is the time from the sample’s global maximum to its endpoint, not the longest completed drawdown. An unusually long episode can reflect bad luck, an overestimated Sharpe ratio, or a poor model of the return process.

Performance can also deteriorate out of sample while the underlying effects stay unchanged. In the synthetic selection experiment, each candidate’s true mean stays fixed. Choosing the best historical estimate selects some favourable noise along with any real effect; that noise need not recur after selection. Performance can therefore appear to decay without the signal weakening or the market responding to its publication.

In his study Does Empirical Finance Replicate?, Muravyev (2026) reports deterioration beyond the original sample and examines its timing. Its empirical evidence concerns a selected, replicable subset of the literature. The synthetic experiment isolates one possible mechanism for a drop; it does not identify the cause of the drop in those empirical results.

What to notice

To assess an indicator ranking, specify what it measures, compare it with a null and account for the search that produced it. If losses follow, check whether sampling variation or changes in the model’s relationships could explain them. Dividing the loss into additive contributions cannot answer these statistical questions on its own.

↑ Sections

7. Research directions

Start by reconciling the ledger and computing markout curves from fills and mids. Then use state indicators recorded at decision time to explain the components conditionally. To allocate P&L to signals, we also need a policy we can rerun. The no-trade band in section 3 explains why separate signal backtests cannot simply be added. The proposals below specify a first deliverable, success criteria and possible negative results; they make no claims about a particular trading operation.

The projects of this section in order, with their prerequisites and the sections they build on.

7.1 Ledger, benchmark and state indicators

Open problem. A fill ledger can reconcile without uniquely identifying spread capture, adverse selection and inventory risk. Menkveld’s accounting and Perold’s decision-time paper portfolio use different references. Neither tells us which paper portfolio to use when quotes are continuously revised. (Menkveld, 2013; Perold, 1988)

First increment: a ledger with explicit accounting conventions. Use simulated market making with a known adverse-selection mechanism. Reconcile identical fills against the own-venue mid, a consolidated mid and future marks. Keep total wealth, hedge P&L, fees and rebates fixed; report markout and positioning per fill and per unit traded. Measure how the components move when the reference changes. Before estimating a planted parameter from a component, state the estimator and the assumptions that make it identifiable.

Deliver the reconciled ledger, a table showing how the reference changes each component, and the estimator’s uncertainty across repeated simulations. Success means that the ledger reconciles at every reference and the declared estimator recovers the planted parameter within a tolerance set before the simulations, with coverage checked against its declared confidence level. A negative result is that no tested reference identifies the parameter, or that apparently stable component labels change under plausible reference choices. Report this result even if the ledger reconciles; reconciliation alone cannot identify a cause.

Second increment: explain markout components by state indicators with search control. Fix the ledger convention, the markout horizons and the decision-time information set. Fit conditional models of each component using spread, volatility, flow and inventory, then assess predictive value on later days. Compare indicator allocations with null distributions that preserve the relevant temporal and cross-indicator dependence. Record every indicator, threshold and horizon tried, and control the declared error rate over that family rather than only the final ranking.

Success means recovering planted conditional relationships in controlled simulations with calibrated uncertainty and error control, then reporting which relationships improve held-out component predictions. A negative result is that no relationship passes the null comparison and search correction, or that apparent gains disappear on new data. Under this design, we would have no evidence that these indicators improve predictions of the components.

Data and dependency. Fills with time, side, price, size and execution status; own and market quotes; mids; inventory and hedge paths; fee and rebate records; indicators as available at each decision. The conditional step depends on the first ledger and on a research log started before model selection. Test instrument-level results as well as pooled results.

Later increments. In their paper on incorporating signals into optimal trading, Lehalle & Neuman (2019) write the optimal execution schedule as a liquidation part plus a signal part; its impact cost then contains a cross term between the two. Compare the two orders of attribution (liquidation first, then signal, and the reverse), then share the cross term by a stated coalition rule. They supply the execution model; the allocation of its cross term is a proposed experiment, not their established attribution result. A separate experiment can test Menkveld’s split of mid price changes into a permanent part (adverse selection) and a transitory part (price pressure), in a simulation where both parts are known. Neither extension is needed to deliver the first ledger.

7.2 Proxy, coalition and exact signal allocations, with uncertainty

Open problem. Even with the same inputs and trading algorithm, stand-alone proxies, coalition reruns and an internal exact split need not answer the same counterfactual question. Removing a signal can change risk scaling, trading dates and active constraints. Before applying Shapley’s rule, we need to choose which game to allocate. (Shapley, 1953)

First increment. On a controlled simulator with correlated signals, compare the undamped NNLS proxy method, coalition Shapley under each stated risk convention, and Durin’s internal allocation where its equations apply. Keep the coalition baseline explicit: a policy with all signals off may still hedge or unwind initial inventory. Allocate the increment above that baseline, not an assumed zero. Compare residuals, reconciliation, distances between shares and sensitivity to correlation and binding constraints. (Durin, 2014)

In this first comparison, separate the numerical error of coalition sampling from uncertainty due to the observed trading history. Calibrate block resampling, placebo indicators and null comparisons on repeated simulated histories with serial dependence and seasonality before using them on a single history. Even a stable-looking attribution needs a null comparison. (Méloux & Dirupo, 2025)

Acceptance and negative result. Deliver a comparison that reproduces the fixed-risk-aversion, unconstrained, no-spread equality, reconciles every method to its stated target including residuals, and reports uncertainty calibration against known simulated allocations. Outside that equality case, Durin is a reference convention, not ground truth for every method. Success means finding the conditions under which proxy error stays within a predeclared tolerance and uncertainty estimates are calibrated. A negative result is finding no such conditions under realistic constraints, or finding unstable proxy shares despite a small residual.

Data and scope. Decision-time signals, positions, costs, risk settings, initial inventory and a policy that can be rerun; the exact comparison also needs the optimiser and its multipliers. Aggregate P&L alone cannot recover the coalition experiment. Predictive games, causal adjustment and learned surrogates are later increments, not requirements of this first comparison.

7.3 Later extensions: predictive games, causal structure and policies

Once the ledger targets and rerun conventions are established, compare regression- $R^2$ sharing and SAGE on the same component histories. Distinguish reductions in prediction loss from allocations in currency. Success would be explaining ranking differences by the target or estimation error; a negative result would be rankings too uncertain to distinguish. This requires held-out component outcomes and fitted predictive models. (Covert et al., 2020)

Causal adjustment needs a separate experiment with simulated event sequences whose causal structure is known. Compare a fit-selected control set with an adjustment set justified by that structure. Inventory and fills may be consequences of a signal; timestamps alone do not identify the graph. Success means recovering the planted causal effect under the declared identification assumptions. A negative result is that the available states cannot identify it, even with good predictive fit. This requires decision, fill and inventory timestamps in addition to the first project’s data. (Lopez de Prado & Zoonekynd, 2024)

Learned policies and surrogates are another later experiment. Practitioner “black-box” trading includes rule-based automation; it need not mean a learned position-taking policy. Guéant’s chapter in the volume Machine Learning and Data Sciences for Financial Markets, edited by Capponi & Lehalle (2023), characterises published RL trading work as proof of concept (ch. 11, p. 200; PDF p. 226), while the sources provide no adoption share. In his book Algorithmic Trading Methods, Kissell (2021) describes a network that gives an optimiser a starting point; the optimiser still runs to find the solution. With rerun access, train a surrogate, a model that predicts the optimiser’s positions from its inputs, on solved cases and measure its position and P&L error before comparing proxy and coalition allocations; SHAP of its position predictions has a different target. Success requires accurate policy predictions on held-out cases and attribution discrepancies within predeclared tolerances. If the surrogate fails to reproduce the policy, report this before interpreting its attributions. The experiment’s conclusions are limited by what the sources report about practice (Isichenko, 2021; Paleologo, 2025) and by the target of the explanation (Lundberg & Lee, 2017). The practitioner overview of graph neural networks by Antulov-Fantulin & Kolm (2026) discusses interpretability concerns. The book Trades, Quotes and Prices by Bouchaud et al. (2018) discusses market impact, which historical replay omits.

7.4 Statistics of a recorded research process

Prerequisite. Start recording every tried variant, rejected specification, parameter choice, evaluation date and P&L history before the research begins. Keep this log as the research proceeds: no correction can recover experiments that were never recorded. The log is also a prerequisite for the state-indicator search in project 7.1.

Open problem. Assessing a selected indicator requires accounting for all the variants tried. Correlated variants make an effective trial count hard to estimate, while backtest-overfitting diagnostics depend on the candidate set supplied. White and Romano–Wolf resample the family together. We still need to test whether that resampling is calibrated for intraday dependence and seasonality. (Bailey et al., 2016; Bailey & López de Prado, 2014; Romano & Wolf, 2005; White, 2000)

First increment. Generate families of strategies and component models with known nulls, planted effects, redundant candidates, seasonality and inventory-driven autocorrelation. Freeze the search and evaluation protocol, then measure rejection rates, missed effects and interval coverage across independent simulated days or histories. Compare effective-count corrections with family-level resampling. Separate a test that any strategy beats a benchmark from an indicator’s conditional contribution test.

Acceptance and negative result. Deliver a table of measured error rates and power for each declared family and dependence design, with Monte Carlo uncertainty. Success requires the stated error criterion to be met within that uncertainty, while retaining useful power at effect sizes chosen in advance. A negative result is excess false discoveries, unreliable coverage, or negligible power at those effects. Missing research history prevents a claim about the unrecorded search; it must remain an explicit limitation.

Later checks. Assess Sharpe aggregation and uncertainty against repeated simulated days rather than treating intraday observations as independent (Lo, 2002). Freeze models at research-window ends and measure changes in later windows, and repeat the measurement at placebo boundaries, dates where no model changed, to distinguish selection effects from ordinary variation. Muravyev’s findings motivate checking whether performance holds up on later data; they supply no intraday decay rate (Muravyev, 2026). Compare forecast decay with realised information coefficients, following the analysis of signal horizons in Breadth, Skill, and Time by Grinold & Kahn (2011), before interpreting short-lived shares as failed research. Finally test error control for the diagnostic in section 0 across order descriptors and time slices, separating continuous slippage prediction from binary bad-order detection (Azencott et al., 2014).

Data. The prospective research log, aligned indicator and return histories, recorded model introduction dates and untouched later evaluation windows. The diagnostic extension additionally needs order-level outcomes and market states; it depends on the ledger and state-indicator work in project 7.1.

↑ Sections

References

1953 L. S. Shapley. A value for n-person games. Princeton University Press. 3, 7
1985 Gary P. Brinson; Nimrod Fachler. Measuring non-U.S. equity portfolio performance. The Journal of Portfolio Management. 2
1986 Gary P. Brinson; L. Randolph Hood; Gilbert L. Beebower. Determinants of Portfolio Performance. Financial Analysts Journal. 2
1988 André F. Perold. The implementation shortfall: Paper versus reality. The Journal of Portfolio Management. 5, 7
1992 William F. Sharpe. Asset allocation: Management style and performance measurement. The Journal of Portfolio Management. 2
1995 Yoav Benjamini; Yosef Hochberg. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society Series B: Statistical Methodology. 6
1996 Robert Litterman. Hot Spots and Hedges. The Journal of Portfolio Management. 2
2000 Halbert White. A Reality Check for Data Snooping. Econometrica. 6, 7
2000 Richard C. Grinold; Ronald N. Kahn. Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Controlling Risk. McGraw-Hill. 1, 2
2001 Stan Lipovetsky; Michael Conklin. Analysis of regression in game theory approach. Applied Stochastic Models in Business and Industry. 4
2001 Yoav Benjamini; Daniel Yekutieli. The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics. 6
2002 Andrew W. Lo. The Statistics of Sharpe Ratios. Financial Analysts Journal. 6, 7
2005 Joseph P. Romano; Michael Wolf. Stepwise Multiple Testing as Formalized Data Snooping. Econometrica. 6, 7
2011 Richard C. Grinold; Ronald N. Kahn. Breadth, Skill, and Time. The Journal of Portfolio Management. 7
2013 Albert J. Menkveld. High frequency trading and the new market makers. Journal of Financial Markets. 5, 7
2013 Olivier Guéant; Charles-Albert Lehalle; Joaquin Fernandez-Tapia. Dealing with the inventory risk: a solution to the market making problem. Mathematics and Financial Economics. 5
2014 Art B. Owen. Sobol' Indices and Shapley Value. SIAM/ASA Journal on Uncertainty Quantification. 4
2014 Bruno Durin. Signal-wise performance attribution for constrained portfolio optimisation. arXiv. 3, 7
2014 David H. Bailey; Marcos López de Prado. The Deflated Sharpe Ratio: Correcting for Selection Bias, Backtest Overfitting, and Non-Normality. The Journal of Portfolio Management. 6, 7
2014 R. Azencott; A. Beri; Y. Gadhyan; N. Joseph; C.-A. Lehalle; M. Rowley. Real-time market microstructure analysis: online transaction cost analysis. Quantitative Finance. 0, 5, 7
2015 Rina Foygel Barber; Emmanuel J. Candès. Controlling the false discovery rate via knockoffs. The Annals of Statistics. 4
2016 Campbell R. Harvey; Yan Liu; Heqing Zhu. ... and the Cross-Section of Expected Returns. The Review of Financial Studies. 6
2016 David H. Bailey; Jonathan M. Borwein; Marcos López de Prado; Qiji Jim Zhu. The Probability of Backtest Overfitting. Journal of Computational Finance. 6, 7
2017 Adam Rej; Philip Seager; Jean-Philippe Bouchaud. You are in a drawdown. When should you start worrying?. arXiv. 6
2017 Scott Lundberg; Su-In Lee. A Unified Approach to Interpreting Model Predictions. arXiv. 4, 7
2018 Jean-Philippe Bouchaud; Julius Bonart; Jonathan Donier; Martin Gould. Trades, quotes and prices: financial markets under the microscope. Cambridge University Press. 7
2019 Charles-Albert Lehalle; Eyal Neuman. Incorporating signals into optimal trading. Finance and Stochastics. 7
2020 Ian Covert; Scott Lundberg; Su-In Lee. Understanding Global Feature Contributions With Additive Importance Measures. arXiv. 4, 7
2021 Michael Isichenko. Quantitative portfolio management: The art and science of statistical arbitrage. John Wiley & Sons. 3, 7
2021 Robert Kissell. Algorithmic Trading Methods: Applications Using Advanced Statistics, Optimization, and Machine Learning Techniques. Academic Press. 7
2023 Agostino Capponi; Charles-Albert Lehalle. Machine Learning and Data Sciences for Financial Markets: A Guide to Contemporary Practices. Cambridge University Press. 7
2024 Marcos Lopez de Prado; Vincent Zoonekynd. Why Has Factor Investing Failed?: The Role of Specification Errors. theoretical. 6, 7
2025 Giuseppe A. Paleologo. The Elements of Quantitative Investing. John Wiley & Sons. 1, 2, 7
2025 Maxime Méloux; Giada Dirupo. The Dead Salmons of AI Interpretability. arXiv. 4, 7
2026 Dmitriy Muravyev. Does Empirical Finance Replicate?. empirical. 6, 7
2026 Nino Antulov-Fantulin; Petter N. Kolm. Graph Neural Networks in Asset and Investment Management: An Overview for Practitioners. The Journal of Portfolio Management. 7

0. From diagnosing orders to explaining a P&L

1. What it means to attribute a P&L

1.1 What daily snapshots miss

1.2 Information and accounting interval

1.3 Two kinds of indicator

2. The classical answers

2.1 Allocation, selection, interaction

Segment A · effects in basis points

Segment B · effects in basis points

Total · effects in basis points

2.2 Explained variance is not a share of money

Average time-series R² · actual versus policy

Cross-plan dispersion · annual arithmetic mean returns

2.3 A fitted exposure need not be a holding

2.4 Factor return and factor risk

Position sizes · capital units

Euler contributions · capital units of annual volatility

2.5 Estimated components have uncertainty

3. Attributing a P&L to signals

3.1 Stand-alone portfolios as proxies

3.2 Removing one signal from the full portfolio

3.3 Using the optimiser’s internal equations

3.4 Sharing value across coalitions

3.5 What the synthetic comparison finds

Final cumulative allocation by method

Reproduce the synthetic signal comparison

1. Optimiser

Diagonal oracle

2. Durin attribution

3. Other methods

4. Reproduction and checks

Check the two published paths

4. Correlated indicators: which value?

4.1 Shared information is not a regression coefficient

4.2 Choose the value before choosing the allocation

4.3 Importance needs a comparison under the null

5. The market maker’s ledger

5.1 Markouts and the reference price

5.2 From state indicators to a component explanation

6. Statistical pitfalls

6.1 Testing the best of many candidates

6.2 Test the strategy family

6.3 More controls can damage the explanation

6.4 A Sharpe ratio is estimated, then often selected

6.5 Losses do not identify their own cause

7. Research directions

7.1 Ledger, benchmark and state indicators

7.2 Proxy, coalition and exact signal allocations, with uncertainty

7.3 Later extensions: predictive games, causal structure and policies

7.4 Statistics of a recorded research process

References