By Islam BarakaBuying $100 of one asset and selling $100 of another looks self-financing. Add a 0.1% buying fee and a 0.2% selling fee, and the trade needs another 30 cents. That small amount exposes the whole modeling problem: where does the cost come from, and does it change which allocation is best?
In portfolios containing Egyptian, Japanese, and US equities expressed in USD, I have had to reconcile units and FX before comparing values. I bring the same habit to a cost-aware rebalance: the proposed holdings, cash, and fee must reconcile in one currency. A precise fee applied to an invented trade is still wrong.
This tutorial is for readers who know portfolio weights, signed trades, and the mean-variance objective. We will fund a small trade, derive a cost-aware optimum, and explain why doing nothing can be the best modeled decision. All example amounts, returns, and fee rates are synthetic.
The decision and its units
Let x⁰ be current holdings in US-dollar value at one valuation time, and let
x be the holdings after one rebalance. The signed trade is t=x−x⁰. A
buying rate c⁺ᵢ and selling rate c⁻ᵢ are dollars of cost per dollar traded
for this rebalance:
C(t) = Σᵢ c⁺ᵢ max(tᵢ,0) + c⁻ᵢ max(−tᵢ,0).
The budget is not optional bookkeeping:
Σᵢ xᵢ + C(t) = Σᵢ x⁰ᵢ, xᵢ ≥ 0.
The explicit cash coordinate is the residual after risky holdings and their cost have been funded. The cost rates apply once to this rebalance. Expected returns and covariance describe the chosen investment period after it. Compare all terms over that same decision period: an annual trading-cost estimate cannot be used as a per-trade rate, and a one-time commission should not be annualized again.
The MOSEK Portfolio Optimization Cookbook’s transaction-cost chapter introduces the positive/negative trade parts, explains that costs must be aggregated over the investment period, and writes the self-financing equality. The convex portfolio reference by Boyd and collaborators places return, risk, and trading costs in a convex portfolio framework. These references describe model structure; they do not supply a broker’s fee table, tax policy, settlement rule, or FX conversion.
Use current, base-currency position values in this equation. Historical total-return prices are useful for return analysis, but they are not a substitute for current marks plus an explicit cash account. The split-adjustment and dividend-adjustment lessons explain why those data bases need to stay identifiable.
Figure 1. The fee is a dollar outflow in the same equation as the holdings. Open the ledger visual at full size.
A utility that charges the fee once
For the reference utility, μ is expected return per dollar over one declared
investment period, Σ is a positive-semidefinite covariance matrix of those returns (it produces nonnegative portfolio variance), and λ
is a nonnegative risk-aversion coefficient with the reciprocal units needed to
combine dollar return and dollar quadratic risk. We maximize
μᵀx − λ xᵀΣx − C(t).
It can look as though the fee has been charged twice because it appears in both equations. Write the expected terminal wealth before applying the risk penalty:
Expected terminal wealth = Σx + μᵀx
= W₀ − C(t) + μᵀx
Expected gain over W₀ = μᵀx − C(t).
The budget determines how many dollars remain invested. The expected-gain objective measures their gain relative to the original wealth W₀. These are consistent. If your objective already uses Σx+μᵀx, subtracting C(t) again would double-count the same fee. The risk penalty then expresses the model's tradeoff between gain and variance; its coefficient is a chosen input, not a universal setting.
The reference scope is long-only, one-rebalance, one base currency, separable linear variable costs, and a zero-return, zero-risk, zero-cost cash account. It excludes fixed per-trade fees, nonlinear market impact, taxes, borrowing, shorting, and multi-period settlement. A fixed fee is not a renamed linear rate, and a market-impact epigraph is not trusted unless its tightness is proved or audited.
The synthetic ledger
Consider current holdings
x⁰ = [600, 400, 100] USD
and the controlled trade
t = [100, −100, −0.30] USD
c⁺ = c⁻ = [0.001, 0.002, 0].
The direct fee is
C(t) = .001×100 + .002×100 = .30 USD.
The post-trade holdings are [700,300,99.70]. Therefore
Σx = 1,099.70 USD
Σx + C(t) = 1,100.00 USD.
The cash line is nonnegative and the complete equation is checkable without trusting an optimizer’s hidden slack. I want this ledger beside the allocation, because a fee that is merely reported after the solve may already have made the proposed purchase impossible.
Why trade parts need an audit
It is convenient to introduce nonnegative t⁺ and t⁻ with
t=t⁺−t⁻. But recombination alone does not prove that they are the true
positive and negative parts. For the same signed trade, these two descriptions
both recombine:
| representation | buy part | sell part | reported cost |
|---|---|---|---|
| exact | [100,0,0] | [0,100,0.30] | $0.30000 |
| slack | [100.25,0,0] | [0.25,100,0.30] | $0.30050 |
The second pair has unnecessary opposing parts. The package reports direct cost from the signed trade and separately flags nonminimal supplied parts. This is the same kind of contract discipline needed when upstream units, marks, splits, dividends, or FX are reconciled: an invalid input can create a trade or a fee that the economic state did not create.
Figure 2. Recombination is necessary but not sufficient for cost tightness. Open the auxiliary audit at full size.
What the optimizer returns
Use the same current holdings and costs, with μ=[.08,.04,0], Σ=diag(.04,.01,0), and λ=.0005. The last coordinate is cash. The resulting dollar holdings are approximately:
| Quantity | Cost-aware result | Zero-cost comparison |
|---|---|---|
| Asset A | $957.39 | $1,020.00 |
| Asset B | $141.74 | $80.00 |
| Cash | $0.00 | $0.00 |
| Fee | $0.873904 | $0.00 |
| Risk-adjusted expected-gain score | $62.954368 | $63.960000 |
The table rounds for reading; the solver and fixture retain full precision. You should not round dollar allocations before checking their budget.
Here is how to reproduce the candidate. In the region that buys A and sells B, with amounts x_A≥600 and x_B≤400,
C = .001(x_A−600) + .002(400−x_B)
= .001x_A − .002x_B + .20.
At zero remaining cash, the budget becomes 1.001x_A+.998x_B=1099.8. One extra dollar of A consumes $1.001 including its fee; retaining one extra dollar of B reduces the sale fee by .002, so it consumes only $.998 of this regional budget. Substitute the same cost into the objective:
U = .079x_A + .042x_B − .00002x_A² − .000005x_B² − .20.
For an interior optimum in this buy/sell region, marginal utility per dollar of budget must match. The Karush–Kuhn–Tucker (KKT) conditions express that balance using a budget multiplier ν:
.079 − .00004x_A = 1.001ν
.042 − .00001x_B = .998ν
1.001x_A + .998x_B = 1099.8.
Solving these three linear equations gives x_A≈957.385384 and x_B≈141.740712. They satisfy the assumed buy/sell directions. Direct fee recomputation gives .8739039598, and x_A+x_B+C=1100. This verifies the regional candidate; the reference solver's global gap also checks that another region or a no-trade boundary cannot improve the objective beyond tolerance.
A meaningful no-trade decision
Costs can also be the reason to do nothing. Set current risky/cash holdings to
[50,50], expected return to [-.03,0], risk aversion to zero, and both
buy/sell rates on the risky coordinate to .10. Selling one dollar avoids
only .03 dollars of expected loss but pays .10; buying is worse still.
The exact optimum is therefore the kink x=[50,50], with cost $0 and
objective −$1.50. The package’s composite gap retains the piecewise cost in
its oracle, so this ordinary no-trade optimum is certified rather than
mistakenly reported as a nonconverged midpoint-subgradient result.
A fee can change the allocator’s direction
The fee is not merely a final subtraction. In the canonical fixture, buying
asset 1 costs .001 per dollar while selling asset 2 costs .002. Those rates
shift the KKT marginal-return equations and the self-financing coefficients to
1.001 and .998; the optimizer returns about $957.39/$141.74, not the
zero-cost $1,020/$80 mix. If the rates or their horizon change, the decision
can change even when μ and Σ do not. This is why a rate without a declared
horizon is rejected.
Figure 3. The fee changes the allocation and remains visible as a ledger quantity. Open the comparison visual at full size.
Open the interactive cost ledger and make two predictions before looking at the returned rows:
- With a zero fee multiplier, the canonical fixture should reach
[1020,80,0]. Restoring the original fees should move it back to approximately[957.39,141.74,0]. - In the no-trade scenario, a negative expected return does not automatically imply a sale. Compare the expected loss avoided with the fee paid:
.03<.10, so the unchanged portfolio wins in this model.
For each result, add the final holdings and the directly computed fee. The sum must recover starting wealth. That reconciliation is more informative than a bar chart alone.
The data-to-ledger path is summarized below; the source is also available at
visuals/mermaid/transaction-cost-flow.mmd.
Running the paired implementation
From the topic directory, the tests exercise the JSON fixture in Python and the TypeScript implementation separately:
python -X utf8 -B -m unittest discover -s tests -p 'test_*.py' -v
npx --yes tsx tests/transaction-cost-aware.test.ts
The JSON adapter reads mu, covariance, lambda_risk, x0, buy_cost,
sell_cost, and optional cash_index, cost_horizon, and max_iterations.
It returns the signed trade, direct cost, cash balance, budget residual,
objective, and a composite Frank–Wolfe gap, which bounds possible improvement while retaining the cost kinks. A max-iteration result is
not labelled optimal, and invalid or non-finite input returns a structured
error rather than a numeric infinity.
What this does not know
Before I accept a fee-aware allocation, I use three quick questions:
| Question | Direct check | Why it changes the decision |
|---|---|---|
| Is the starting portfolio real on one base-currency basis? | Reconcile units, marks, corporate actions, FX, and current cash | A fictitious trade creates a fictitious fee. |
| Is the fee funded at the same rebalance horizon? | Recompute Σx+C(t)=W₀ from signed trades | A post-solve fee can make the purchase unfundable. |
| Is doing nothing allowed and honestly compared? | Evaluate the no-trade candidate at the cost kink | Positive costs can make the current portfolio optimal. |
The optimizer receives already validated dollar holdings. It does not decide whether a price series is split-adjusted, whether cash dividends have been received, whether a tax lot is eligible, or which direction an FX quote uses. Adjusted total-return prices and separately credited cash must not both be added to the same wealth path. A tax-aware model must retain acquisition basis, current mark, units, and tax cash state rather than applying this fee rate as a universal tax calculator.
I keep one habit from this lesson: inspect the fee beside the trade and cash balance. If doing nothing is feasible, compare it explicitly. A complicated allocation does not deserve preference just because it contains more trades.
Continue to Cardinality-Constrained Portfolios when the question becomes how many holdings may remain. The turnover tutorial addresses a different limit: how much weight may move.
The examples are synthetic. A named historical case remains deferred until a dated fee schedule, position record, and settlement ledger can be independently checked. The canonical contract and primary references document the implementation and source boundaries.
ReferencesPrimary sources and evidence notesExpand the source trail, evidence role, and limitations behind the engineering choices.
Expand the source trail, evidence role, and limitations behind the engineering choices.
MOSEK Portfolio Optimization Cookbook, Chapter 6
- Organization or authors: MOSEK ApS.
- Source type: Official solver documentation and portfolio-optimization cookbook.
- Publication or effective date: Online version accessed 2026-09-14.
- Version: Portfolio Optimization Cookbook 1.6.0, Chapter 6, sections 6 and 6.1.
- URL or DOI: MOSEK transaction costs
- Accessed: 2026-09-14.
- Jurisdiction: None; this is mathematical/solver documentation.
- Supports: Variable transaction costs, positive and negative trade parts, cost aggregation over an investment period, and the self-financing equality
1ᵀx + C(x−x0) = 1ᵀx0. - Limitations: The chapter’s examples do not provide a universal market fee schedule. Fixed costs, market impact, and buy-in thresholds are separate models. Its §6.4 cardinality example counts nonzero trades, not final holdings; it is not used as a holdings-cardinality claim here. The market-impact epigraph warning motivates the direct-cost audit in this package.
Boyd et al., Multi-Period Trading via Convex Optimization
- Organization or authors: Stephen Boyd and collaborators.
- Source type: University technical paper and maintained reference page.
- Publication or effective date: 2017 record; accessed 2026-09-14.
- Version:
cvx_portfolio.htmllanding page. - URL or DOI: Stanford convex portfolio page
- Accessed: 2026-09-14.
- Jurisdiction: None; this is a mathematical optimization reference.
- Supports: Convex portfolio objectives with expected return, risk, transaction costs, holding costs, and trading over portfolio positions.
- Limitations: It does not determine a user’s broker fees, taxes, settlement convention, or currency. Those are explicit data-contract inputs and are outside this one-rebalance reference.
Google discovery status
The A02 Google top-ten ledger remains queued in the bounded wave. Search-engine observations will be recorded separately with query, locale, date, observed order, page-review extent, and access blockers; no other engine will be relabeled as Google.
Full dependency-light reference implementations in both supported languages.
/**
* Self-financing linear transaction-cost portfolio optimizer.
*
* Dollar holdings use one explicit non-negative cash coordinate. The
* reduced risky holdings z=x/W0 solve
* maximize mu·z - lambdaRisk*W0*z' Sigma z - C(z-z0)
* subject to z>=0 and sum(z)+C(z-z0)<=1.
* The residual is the cash balance, so the full result satisfies
* sum(x)+C(x-x0)=sum(x0) exactly up to floating-point tolerance.
* The certificate retains the exact piecewise cost at a kink instead of using
* a midpoint subgradient as a false global gap.
*/
export const ABSOLUTE_TOLERANCE = 1e-9;
export const RELATIVE_TOLERANCE = 2e-8;
export const FEASIBILITY_TOLERANCE = 2e-10;
export const PSD_RELATIVE_TOLERANCE = 1e-10;
export const CURRENT_SUM_TOLERANCE = 1e-10;
export const RATE_TOLERANCE = 1e-12;
const KINK_SNAP_TOLERANCE = 1e-8;
export const DEFAULT_MAX_ITERATIONS = 5000;
const PROJECTION_BISECTION_ITERATIONS = 80;
const ORACLE_REGION_LIMIT = 16;
export class CostAwareValidationError extends Error {
readonly code: string;
constructor(code: string, message: string) {
super(message);
this.name = "CostAwareValidationError";
this.code = code;
}
}
export class CostAwareNumericalError extends Error {
readonly code: string;
constructor(code: string, message: string) {
super(message);
this.name = "CostAwareNumericalError";
this.code = code;
}
}
export interface CostAwareCertificate {
kind: "composite-piecewise-frank-wolfe-gap";
frankWolfeGap: number;
normalizedFrankWolfeGap: number;
objectiveScale: number;
globalOptimality: string;
}
export interface CostAwareResult {
status: "optimal" | "max-iterations";
weights: number[];
trade: number[];
buy: number[];
sell: number[];
directCost: number;
auxiliaryCost: number;
budgetResidual: number;
objective: number;
iterations: number;
cashIndex: number;
costHorizon: string;
currentWealth: number;
cashBalance: number;
certificate: CostAwareCertificate;
}
function invalid(code: string, message: string): never {
throw new CostAwareValidationError(code, message);
}
function asReal(value: unknown, name: string): number {
if (typeof value !== "number" || !Number.isFinite(value)) {
invalid("invalid-input", `${name} must be a finite real number.`);
}
return value;
}
function vector(value: unknown, name: string): number[] {
if (!Array.isArray(value) || value.length === 0) {
invalid("invalid-input", `${name} must be a non-empty numeric array.`);
}
return value.map((entry, index) => asReal(entry, `${name}[${index}]`));
}
function matrix(value: unknown, name: string, size: number): number[][] {
if (!Array.isArray(value) || value.length !== size) {
invalid("invalid-input", `${name} must have ${size} rows.`);
}
return value.map((row, rowIndex) => {
if (!Array.isArray(row) || row.length !== size) {
invalid("invalid-input", `${name} must be ${size} by ${size}.`);
}
return row.map((entry, columnIndex) => asReal(entry, `${name}[${rowIndex}][${columnIndex}]`));
});
}
function stableSum(values: readonly number[]): number {
let total = 0;
let compensation = 0;
for (const value of values) {
const candidate = total + value;
if (Math.abs(total) >= Math.abs(value)) compensation += (total - candidate) + value;
else compensation += (value - candidate) + total;
total = candidate;
}
return total + compensation;
}
function dot(left: readonly number[], right: readonly number[]): number {
return stableSum(left.map((value, index) => value * right[index]));
}
function matVec(matrixValue: readonly (readonly number[])[], value: readonly number[]): number[] {
return matrixValue.map((row) => dot(row, value));
}
function quadratic(matrixValue: readonly (readonly number[])[], value: readonly number[]): number {
return dot(value, matVec(matrixValue, value));
}
function minEigen(matrixValue: readonly (readonly number[])[]): number {
const size = matrixValue.length;
const work = matrixValue.map((row) => [...row]);
if (size === 1) return work[0][0];
const scale = Math.max(0, ...work.flat().map((entry) => Math.abs(entry)));
for (let sweep = 0; sweep < Math.max(60, 20 * size * size); sweep += 1) {
let p = 0; let q = 1; let largest = 0;
for (let row = 0; row < size; row += 1) {
for (let column = row + 1; column < size; column += 1) {
const magnitude = Math.abs(work[row][column]);
if (magnitude > largest) { largest = magnitude; p = row; q = column; }
}
}
if (largest <= 1e-14 * scale) break;
const app = work[p][p]; const aqq = work[q][q]; const apq = work[p][q];
const angle = 0.5 * Math.atan2(2 * apq, aqq - app);
const cosine = Math.cos(angle); const sine = Math.sin(angle);
for (let index = 0; index < size; index += 1) {
if (index === p || index === q) continue;
const aip = work[index][p]; const aiq = work[index][q];
work[index][p] = cosine * aip - sine * aiq; work[p][index] = work[index][p];
work[index][q] = sine * aip + cosine * aiq; work[q][index] = work[index][q];
}
work[p][p] = cosine * cosine * app - 2 * cosine * sine * apq + sine * sine * aqq;
work[q][q] = sine * sine * app + 2 * cosine * sine * apq + cosine * cosine * aqq;
work[p][q] = 0; work[q][p] = 0;
}
return Math.min(...work.map((row, index) => row[index]));
}
interface Validated {
mean: number[]; sigma: number[][]; risk: number; current: number[];
buy: number[]; sell: number[]; cash: number; horizon: string;
}
function validateInputs(
mu: unknown, covariance: unknown, lambdaRisk: unknown, x0: unknown,
buyCost: unknown, sellCost: unknown, cashIndexValue: unknown, horizonValue: unknown,
): Validated {
const mean = vector(mu, "mu");
const sigma = matrix(covariance, "covariance", mean.length);
const risk = asReal(lambdaRisk, "lambda_risk");
const current = vector(x0, "x0");
const buy = vector(buyCost, "buy_cost");
const sell = vector(sellCost, "sell_cost");
if (current.length !== mean.length || buy.length !== mean.length || sell.length !== mean.length) {
invalid("invalid-input", "mu, covariance, x0, buy_cost, and sell_cost must align.");
}
if (typeof cashIndexValue !== "number" || !Number.isInteger(cashIndexValue)) {
invalid("invalid-input", "cash_index must be an integer coordinate.");
}
const cash = cashIndexValue === -1 ? mean.length - 1 : cashIndexValue;
if (cash < 0 || cash >= mean.length) invalid("invalid-input", "cash_index is outside the holdings vector.");
if (horizonValue !== "per-rebalance") {
invalid("cost-horizon-required", "this reference scope requires cost_horizon='per-rebalance'; convert annual rates before calling.");
}
if (risk < 0) invalid("invalid-input", "lambda_risk must be nonnegative.");
if (current.some((entry) => entry < 0)) invalid("invalid-input", "x0 must be nonnegative.");
const wealth = stableSum(current);
if (!Number.isFinite(wealth)) invalid("wealth-overflow", "x0 wealth must be finite after summation.");
if (wealth <= 0) invalid("invalid-input", "x0 must have positive wealth.");
if (buy.some((entry) => entry < 0) || sell.some((entry) => entry < 0)) {
invalid("invalid-input", "buy_cost and sell_cost must be nonnegative.");
}
if (buy.some((entry) => entry >= 1) || sell.some((entry) => entry >= 1)) {
invalid("cost-rate-out-of-range", "cost rates must be below one for the reduced cash ledger.");
}
if (buy[cash] !== 0 || sell[cash] !== 0) invalid("cash-cost-unsupported", "the explicit cash coordinate must have zero variable cost.");
if (mean[cash] !== 0) invalid("cash-return-unsupported", "the reference cash coordinate must have exactly zero expected return.");
if (sigma[cash].some((entry) => entry !== 0) || sigma.some((row) => row[cash] !== 0)) {
invalid("cash-risk-unsupported", "the reference cash coordinate must have zero covariance row and column.");
}
const matrixScale = Math.max(0, ...sigma.flat().map((entry) => Math.abs(entry)));
const symmetryTolerance = 1e-12 * matrixScale;
for (let row = 0; row < mean.length; row += 1) {
for (let column = 0; column < row; column += 1) {
if (Math.abs(sigma[row][column] - sigma[column][row]) > symmetryTolerance) invalid("indefinite-covariance", "covariance must be symmetric at the matrix scale.");
}
}
if (minEigen(sigma) < -PSD_RELATIVE_TOLERANCE * matrixScale) invalid("indefinite-covariance", "covariance must be positive semidefinite.");
return { mean, sigma, risk, current, buy, sell, cash, horizon: String(horizonValue) };
}
function cost(trade: readonly number[], buy: readonly number[], sell: readonly number[]): number {
return stableSum(trade.map((entry, index) => buy[index] * Math.max(entry, 0) + sell[index] * Math.max(-entry, 0)));
}
function ledger(point: readonly number[], anchor: readonly number[], buy: readonly number[], sell: readonly number[]): number {
return stableSum(point.map((entry, index) => entry + buy[index] * Math.max(entry - anchor[index], 0) + sell[index] * Math.max(anchor[index] - entry, 0)));
}
function regionComponent(value: number, anchor: number, buyRate: number, sellRate: number, multiplier: number): number {
let upperCandidate = Math.max(anchor, value - multiplier * (1 + buyRate));
if (Math.abs(upperCandidate - anchor) <= KINK_SNAP_TOLERANCE) upperCandidate = anchor;
const candidates = [
Math.max(0, Math.min(anchor, value - multiplier * (1 - sellRate))),
upperCandidate,
];
const local = (point: number) => 0.5 * (point - value) ** 2 + multiplier * (point + buyRate * Math.max(point - anchor, 0) + sellRate * Math.max(anchor - point, 0));
return candidates.reduce((best, point) => local(point) < local(best) ? point : best, candidates[0]);
}
function project(values: readonly number[], anchor: readonly number[], buy: readonly number[], sell: readonly number[]): number[] {
const candidate = values.map((value) => Math.max(value, 0));
if (ledger(candidate, anchor, buy, sell) <= 1) return candidate;
let lower = 0; let upper = 1;
const at = (multiplier: number) => values.map((value, index) => regionComponent(value, anchor[index], buy[index], sell[index], multiplier));
while (ledger(at(upper), anchor, buy, sell) > 1) {
upper *= 2;
if (upper > 1e18) throw new CostAwareNumericalError("projection-failed", "projection multiplier exceeded finite bounds.");
}
for (let iteration = 0; iteration < PROJECTION_BISECTION_ITERATIONS; iteration += 1) {
const middle = (lower + upper) / 2;
if (ledger(at(middle), anchor, buy, sell) > 1) lower = middle; else upper = middle;
}
return at((lower + upper) / 2);
}
function costGradient(point: readonly number[], anchor: readonly number[], buy: readonly number[], sell: readonly number[]): number[] {
return point.map((value, index) => value > anchor[index] + RATE_TOLERANCE ? buy[index] : value < anchor[index] - RATE_TOLERANCE ? -sell[index] : (buy[index] - sell[index]) / 2);
}
function objective(point: readonly number[], mean: readonly number[], sigma: readonly (readonly number[])[], risk: number, wealth: number, buy: readonly number[], sell: readonly number[], anchor: readonly number[]): number {
return dot(mean, point) - risk * wealth * quadratic(sigma, point) - cost(point.map((entry, index) => entry - anchor[index]), buy, sell);
}
function solveLinearSystem(matrixValue: number[][], right: number[]): number[] | undefined {
const size = right.length;
const work = matrixValue.map((row, index) => [...row, right[index]]);
const scale = Math.max(1, ...matrixValue.flat().map((entry) => Math.abs(entry)));
for (let column = 0; column < size; column += 1) {
let pivot = column;
for (let row = column + 1; row < size; row += 1) if (Math.abs(work[row][column]) > Math.abs(work[pivot][column])) pivot = row;
if (Math.abs(work[pivot][column]) <= 1e-13 * scale) return undefined;
if (pivot !== column) [work[column], work[pivot]] = [work[pivot], work[column]];
const pivotValue = work[column][column];
for (let row = column + 1; row < size; row += 1) {
const factor = work[row][column] / pivotValue;
if (factor === 0) continue;
for (let entry = column; entry <= size; entry += 1) work[row][entry] -= factor * work[column][entry];
}
}
const solution = Array.from({ length: size }, () => 0);
for (let row = size - 1; row >= 0; row -= 1) {
let remainder = work[row][size];
for (let entry = row + 1; entry < size; entry += 1) remainder -= work[row][entry] * solution[entry];
solution[row] = remainder / work[row][row];
}
return solution;
}
function polishKkt(point: readonly number[], mean: readonly number[], sigma: readonly (readonly number[])[], risk: number, wealth: number, buy: readonly number[], sell: readonly number[], anchor: readonly number[]): number[] {
if (risk <= 0 || ledger(point, anchor, buy, sell) < 1 - 1e-8) return [...point];
const regionTolerance = 1e-8;
if (point.some((value, index) => Math.abs(value - anchor[index]) <= regionTolerance)) return [...point];
const active = point.map((value, index) => value > 1e-10 ? index : -1).filter((index) => index >= 0);
if (active.length === 0) return [...point];
const q: number[] = []; const constant: number[] = []; const linear: number[] = [];
for (let index = 0; index < point.length; index += 1) {
if (point[index] > anchor[index] + regionTolerance) { q.push(1 + buy[index]); constant.push(-buy[index] * anchor[index]); linear.push(mean[index] - buy[index]); }
else if (point[index] < anchor[index] - regionTolerance) { q.push(1 - sell[index]); constant.push(sell[index] * anchor[index]); linear.push(mean[index] + sell[index]); }
else return [...point];
}
const dimension = active.length + 1;
const system = Array.from({ length: dimension }, () => Array.from({ length: dimension }, () => 0));
const right = Array.from({ length: dimension }, () => 0);
active.forEach((index, row) => {
right[row] = linear[index];
active.forEach((other, column) => { system[row][column] = 2 * risk * wealth * sigma[index][other]; });
system[row][dimension - 1] = q[index]; system[dimension - 1][row] = q[index];
});
right[dimension - 1] = 1 - stableSum(constant);
const solution = solveLinearSystem(system, right);
if (!solution) return [...point];
const candidate = Array.from({ length: point.length }, () => 0);
active.forEach((index, row) => { candidate[index] = solution[row]; });
const multiplier = solution[dimension - 1];
for (let index = 0; index < candidate.length; index += 1) {
const value = candidate[index];
if (value < -1e-8) return [...point];
if (point[index] > anchor[index] + regionTolerance && value < anchor[index] - 1e-8) return [...point];
if (point[index] < anchor[index] - regionTolerance && value > anchor[index] + 1e-8) return [...point];
}
if (Math.abs(ledger(candidate, anchor, buy, sell) - 1) > 1e-7) return [...point];
const qGradient = candidate.map((_, index) => linear[index] - 2 * risk * wealth * dot(sigma[index], candidate) - multiplier * q[index]);
if (qGradient.some((entry, index) => !active.includes(index) && entry > 1e-7)) return [...point];
const clipped = candidate.map((value) => Math.max(0, value));
if (objective(clipped, mean, sigma, risk, wealth, buy, sell, anchor) < objective(point, mean, sigma, risk, wealth, buy, sell, anchor) - 1e-8) return [...point];
return clipped;
}
function gradient(point: readonly number[], mean: readonly number[], sigma: readonly (readonly number[])[], risk: number, wealth: number, buy: readonly number[], sell: readonly number[], anchor: readonly number[]): number[] {
const riskGradient = matVec(sigma, point);
const costs = costGradient(point, anchor, buy, sell);
return mean.map((entry, index) => entry - 2 * risk * wealth * riskGradient[index] - costs[index]);
}
function smoothGradient(point: readonly number[], mean: readonly number[], sigma: readonly (readonly number[])[], risk: number, wealth: number): number[] {
const riskGradient = matVec(sigma, point);
return mean.map((entry, index) => entry - 2 * risk * wealth * riskGradient[index]);
}
function linearOracle(gradientValue: readonly number[], anchor: readonly number[], buy: readonly number[], sell: readonly number[]): number[] {
const size = gradientValue.length;
if (size > ORACLE_REGION_LIMIT) throw new CostAwareNumericalError("certificate-too-large", `global gap certificate is limited to ${ORACLE_REGION_LIMIT} risky coordinates.`);
let bestValue = -Infinity; let bestPoint: number[] | undefined;
for (let mask = 0; mask < (1 << size); mask += 1) {
const lower: number[] = []; const rates: number[] = []; const constants: number[] = []; const capacities: Array<number | undefined> = [];
for (let index = 0; index < size; index += 1) {
if (mask & (1 << index)) { lower.push(anchor[index]); rates.push(1 + buy[index]); constants.push(-buy[index] * anchor[index]); capacities.push(undefined); }
else { lower.push(0); rates.push(1 - sell[index]); constants.push(sell[index] * anchor[index]); capacities.push(anchor[index]); }
}
let remaining = 1 - stableSum(constants) - dot(rates, lower);
if (remaining < -FEASIBILITY_TOLERANCE) continue;
const point = [...lower];
const order = [...Array(size).keys()].sort((a, b) => (gradientValue[b] / rates[b] - gradientValue[a] / rates[a]) || a - b);
for (const index of order) {
if (remaining <= FEASIBILITY_TOLERANCE || gradientValue[index] <= 0) continue;
let amount = remaining / rates[index];
if (capacities[index] !== undefined) amount = Math.min(amount, capacities[index] as number);
point[index] += amount; remaining -= rates[index] * amount;
}
const value = dot(gradientValue, point);
if (value > bestValue + 1e-15) { bestValue = value; bestPoint = point; }
}
if (!bestPoint) throw new CostAwareNumericalError("oracle-failed", "no feasible piecewise cost region was found.");
return bestPoint;
}
function compositeOracle(smooth: readonly number[], anchor: readonly number[], buy: readonly number[], sell: readonly number[]): { point: number[]; value: number } {
const size = smooth.length;
if (size > ORACLE_REGION_LIMIT) throw new CostAwareNumericalError("certificate-too-large", `global gap certificate is limited to ${ORACLE_REGION_LIMIT} risky coordinates.`);
let bestValue = -Infinity; let bestPoint: number[] | undefined;
for (let mask = 0; mask < (1 << size); mask += 1) {
const lower: number[] = []; const rates: number[] = []; const constants: number[] = []; const capacities: Array<number | undefined> = []; const costSlopes: number[] = [];
for (let index = 0; index < size; index += 1) {
if (mask & (1 << index)) { lower.push(anchor[index]); rates.push(1 + buy[index]); constants.push(-buy[index] * anchor[index]); capacities.push(undefined); costSlopes.push(buy[index]); }
else { lower.push(0); rates.push(1 - sell[index]); constants.push(sell[index] * anchor[index]); capacities.push(anchor[index]); costSlopes.push(-sell[index]); }
}
let remaining = 1 - stableSum(constants) - dot(rates, lower);
if (remaining < -FEASIBILITY_TOLERANCE) continue;
const point = [...lower];
const effective = smooth.map((entry, index) => entry - costSlopes[index]);
const order = [...Array(size).keys()].sort((a, b) => (effective[b] / rates[b] - effective[a] / rates[a]) || a - b);
for (const index of order) {
if (remaining <= FEASIBILITY_TOLERANCE || effective[index] <= 0) continue;
let amount = remaining / rates[index];
if (capacities[index] !== undefined) amount = Math.min(amount, capacities[index] as number);
point[index] += amount; remaining -= rates[index] * amount;
}
const value = dot(smooth, point) - cost(point.map((entry, index) => entry - anchor[index]), buy, sell);
if (value > bestValue + 1e-15) { bestValue = value; bestPoint = point; }
}
if (!bestPoint) throw new CostAwareNumericalError("oracle-failed", "no feasible piecewise cost region was found.");
return { point: bestPoint, value: bestValue };
}
export function optimizeTransactionCostAware(
mu: unknown, covariance: unknown, lambdaRisk: unknown, x0: unknown,
buyCost: unknown, sellCost: unknown,
options: { cashIndex?: number; costHorizon?: string; maxIterations?: number } = {},
): CostAwareResult {
const checked = validateInputs(mu, covariance, lambdaRisk, x0, buyCost, sellCost, options.cashIndex ?? -1, options.costHorizon ?? "per-rebalance");
const maxIterations = options.maxIterations ?? DEFAULT_MAX_ITERATIONS;
if (!Number.isInteger(maxIterations) || maxIterations <= 0) invalid("invalid-input", "max_iterations must be a positive integer.");
const wealth = stableSum(checked.current);
const riskyIndices = checked.current.map((_, index) => index).filter((index) => index !== checked.cash);
const mean = riskyIndices.map((index) => checked.mean[index]);
const sigma = riskyIndices.map((row) => riskyIndices.map((column) => checked.sigma[row][column]));
const buy = riskyIndices.map((index) => checked.buy[index]);
const sell = riskyIndices.map((index) => checked.sell[index]);
const anchor = riskyIndices.map((index) => checked.current[index] / wealth);
let point = [...anchor];
const maxEigen = sigma.length === 0 ? 0 : Math.max(0, ...sigma.map((row) => stableSum(row.map((entry) => Math.abs(entry)))));
const stepSize = 1 / Math.max(1, 1 + 2 * checked.risk * wealth * maxEigen);
const objectiveScale = Math.max(1e-12, ...mean.map((entry) => Math.abs(entry)), checked.risk * wealth * maxEigen, ...buy, ...sell);
let previous = objective(point, mean, sigma, checked.risk, wealth, buy, sell, anchor);
let iterations = 0;
for (let iteration = 1; iteration <= maxIterations; iteration += 1) {
iterations = iteration;
const g = gradient(point, mean, sigma, checked.risk, wealth, buy, sell, anchor);
let localStep = stepSize;
let trial = project(point.map((entry, index) => entry + localStep * g[index]), anchor, buy, sell);
let trialObjective = objective(trial, mean, sigma, checked.risk, wealth, buy, sell, anchor);
while (trialObjective < previous - 1e-14 * objectiveScale && localStep > 1e-14) {
localStep *= 0.5;
trial = project(point.map((entry, index) => entry + localStep * g[index]), anchor, buy, sell);
trialObjective = objective(trial, mean, sigma, checked.risk, wealth, buy, sell, anchor);
}
const change = Math.max(0, ...trial.map((entry, index) => Math.abs(entry - point[index])));
point = trial; previous = trialObjective;
const newSmooth = smoothGradient(point, mean, sigma, checked.risk, wealth);
const oracle = compositeOracle(newSmooth, anchor, buy, sell);
let gap = oracle.value - (dot(newSmooth, point) - cost(point.map((entry, index) => entry - anchor[index]), buy, sell));
if (gap < 0 && Math.abs(gap) <= 10 * ABSOLUTE_TOLERANCE * objectiveScale) gap = 0;
if (gap <= ABSOLUTE_TOLERANCE * objectiveScale + RELATIVE_TOLERANCE * objectiveScale || change <= 1e-14) break;
}
point = polishKkt(point, mean, sigma, checked.risk, wealth, buy, sell, anchor);
const costRate = cost(point.map((entry, index) => entry - anchor[index]), buy, sell);
const cashFraction = 1 - stableSum(point) - costRate;
if (cashFraction < -FEASIBILITY_TOLERANCE) throw new CostAwareNumericalError("infeasible-result", "cash residual is negative after optimization.");
const weights = Array.from({ length: checked.current.length }, () => 0);
riskyIndices.forEach((index, position) => { weights[index] = wealth * point[position]; });
weights[checked.cash] = wealth * Math.max(0, cashFraction);
const trade = weights.map((entry, index) => entry - checked.current[index]);
const buyParts = trade.map((entry) => Math.max(entry, 0));
const sellParts = trade.map((entry) => Math.max(-entry, 0));
const directCost = cost(trade, checked.buy, checked.sell);
const auxiliaryCost = dot(checked.buy, buyParts) + dot(checked.sell, sellParts);
const budgetResidual = stableSum(weights) + directCost - wealth;
const fullObjective = dot(checked.mean, weights) - checked.risk * quadratic(checked.sigma, weights) - directCost;
const finalSmooth = smoothGradient(point, mean, sigma, checked.risk, wealth);
const finalOracle = compositeOracle(finalSmooth, anchor, buy, sell);
const gapRate = Math.max(0, finalOracle.value - (dot(finalSmooth, point) - cost(point.map((entry, index) => entry - anchor[index]), buy, sell)));
const feasible = Math.min(...weights) >= -FEASIBILITY_TOLERANCE && Math.abs(budgetResidual) <= 10 * FEASIBILITY_TOLERANCE * Math.max(1, wealth);
const status = feasible && gapRate / objectiveScale <= RELATIVE_TOLERANCE ? "optimal" : "max-iterations";
if (![...weights, ...trade, directCost, auxiliaryCost, budgetResidual, fullObjective, gapRate].every(Number.isFinite)) throw new CostAwareNumericalError("nonfinite-result", "solver returned a non-finite diagnostic.");
return {
status, weights, trade, buy: buyParts, sell: sellParts, directCost, auxiliaryCost,
budgetResidual, objective: fullObjective, iterations, cashIndex: checked.cash,
costHorizon: checked.horizon, currentWealth: wealth, cashBalance: weights[checked.cash],
certificate: {
kind: "composite-piecewise-frank-wolfe-gap", frankWolfeGap: gapRate * wealth,
normalizedFrankWolfeGap: gapRate / objectiveScale, objectiveScale: objectiveScale * wealth,
globalOptimality: status === "optimal" ? "global concave optimum certified by piecewise-region linear oracle" : "not certified; inspect status and gap",
},
};
}
export function auditTradeParts(
tradeValue: unknown, buyCost: unknown, sellCost: unknown,
auxiliaryBuy: unknown, auxiliarySell: unknown,
): { recombinedTrade: number[]; directCost: number; auxiliaryCost: number; recombined: boolean; auxiliaryPartsMinimal: boolean } {
const trade = vector(tradeValue, "trade"); const buy = vector(buyCost, "buy_cost"); const sell = vector(sellCost, "sell_cost");
const plus = vector(auxiliaryBuy, "auxiliary_buy"); const minus = vector(auxiliarySell, "auxiliary_sell");
if (![buy, sell, plus, minus].every((entry) => entry.length === trade.length)) invalid("invalid-input", "trade and auxiliary vectors must align.");
if (buy.some((entry) => entry < 0) || sell.some((entry) => entry < 0)) invalid("negative-cost-rate", "buy_cost and sell_cost must be nonnegative.");
if (plus.some((entry) => entry < 0) || minus.some((entry) => entry < 0)) invalid("negative-auxiliary-part", "auxiliary buy and sell parts must be nonnegative.");
const recombinedTrade = plus.map((entry, index) => entry - minus[index]);
return {
recombinedTrade, directCost: cost(trade, buy, sell), auxiliaryCost: dot(buy, plus) + dot(sell, minus),
recombined: recombinedTrade.every((entry, index) => Math.abs(entry - trade[index]) <= ABSOLUTE_TOLERANCE),
auxiliaryPartsMinimal: plus.every((entry, index) => Math.abs(entry - Math.max(trade[index], 0)) <= ABSOLUTE_TOLERANCE && Math.abs(minus[index] - Math.max(-trade[index], 0)) <= ABSOLUTE_TOLERANCE),
};
}
The embedded lab now expands to its full document height, keeping the article as the only scroll surface.