Mean-CVaR Optimization: Make the Tail Auditable

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

export interface MeanCvarResult {
    branch?: string;
    status: string;
    variant: string;
    method: string;
    assetIds: string[] | null;
    weights: number[] | null;
    expectedReturn: number | null;
    cvar: number | null;
    thresholdZ: number | null;
    tailSlacks: number[] | null;
    scenarioLosses: number[] | null;
    objective: number | null;
    budgetResidual: number | null;
    lowerBoundResidual: number | null;
    probabilityResidual: number | null;
    iterations: number;
    maxIterations: number;
    solutionClass: string;
    diagnostics: Record<string, unknown>;
    warnings: string[];
}
const DEFAULT_MAX_ITERATIONS = 10000;
const TOL_PROBABILITY = 1e-12;
const TOL_FEASIBILITY = 1e-10;
const TOL_OBJECTIVE = 1e-10;
const MAX_VERTEX_VARIABLES = 18;
const MAX_VERTEX_CANDIDATES = 200000;
function finiteReal(value: unknown): value is number {
    return typeof value === "number" && Number.isFinite(value);
}
function dot(left: number[], right: number[]): number {
    return left.reduce((sum, value, index) => sum + value * right[index], 0);
}
function solveLinear(matrix: number[][], vector: number[]): number[] | null {
    const size = vector.length;
    const augmented = matrix.map((row, index) => [...row, vector[index]]);
    for (let column = 0; column < size; column += 1) {
        let pivot = column;
        for (let row = column + 1; row < size; row += 1)
            if (Math.abs(augmented[row][column]) > Math.abs(augmented[pivot][column]))
                pivot = row;
        if (Math.abs(augmented[pivot][column]) <= 1e-13)
            return null;
        [augmented[column], augmented[pivot]] = [augmented[pivot], augmented[column]];
        const divisor = augmented[column][column];
        for (let j = column; j <= size; j += 1)
            augmented[column][j] /= divisor;
        for (let row = 0; row < size; row += 1) {
            if (row === column)
                continue;
            const factor = augmented[row][column];
            if (factor === 0)
                continue;
            for (let j = column; j <= size; j += 1)
                augmented[row][j] -= factor * augmented[column][j];
        }
    }
    return augmented.map((row) => row[size]);
}
function combinations(size: number, choose: number): number[][] {
    const output: number[][] = [];
    const current: number[] = [];
    const visit = (start: number): void => {
        if (current.length === choose) {
            output.push(current.slice());
            return;
        }
        for (let index = start; index <= size - (choose - current.length); index += 1) {
            current.push(index);
            visit(index + 1);
            current.pop();
        }
    };
    visit(0);
    return output;
}
export function weightedVarCvar(losses: number[], probabilities: number[], beta: number): {
    threshold: number;
    cvar: number;
} {
    const ordered = losses.map((loss, index) => ({ loss, probability: probabilities[index] })).sort((a, b) => a.loss - b.loss);
    let cumulative = 0;
    let threshold = ordered[ordered.length - 1].loss;
    for (const item of ordered) {
        cumulative += item.probability;
        if (cumulative >= beta) {
            threshold = item.loss;
            break;
        }
    }
    const cvar = threshold + losses.reduce((sum, loss, index) => sum + probabilities[index] * Math.max(0, loss - threshold), 0) / (1 - beta);
    return { threshold, cvar };
}
function failure(status: string, message: string, details: Partial<MeanCvarResult> = {}, diagnostics: Record<string, unknown> = {}): MeanCvarResult {
    return { status, variant: "mean-cvar-positive-gamma-long-only", method: "lp-vertex-enumeration", assetIds: details.assetIds ?? null, weights: null, expectedReturn: null, cvar: null, thresholdZ: null, tailSlacks: null, scenarioLosses: null, objective: null, budgetResidual: null, lowerBoundResidual: null, probabilityResidual: details.probabilityResidual ?? null, iterations: details.iterations ?? 0, maxIterations: details.maxIterations ?? DEFAULT_MAX_ITERATIONS, solutionClass: "unknown", diagnostics: { message, ...diagnostics }, warnings: [] };
}
function meanCvarNormalized(assetIdsInput: unknown, returnsInput: unknown, probabilitiesInput: unknown, betaInput: unknown, gammaInput: unknown, options: {
    maxIterations?: unknown;
} = {}): MeanCvarResult {
    if (!Array.isArray(assetIdsInput))
        return failure("invalid_input", "assetIds must be an ordered array.");
    const assetIds = assetIdsInput.slice();
    if (!assetIds.length || assetIds.some((value) => typeof value !== "string" || value.trim() === ""))
        return failure("invalid_input", "assetIds must contain non-empty strings.");
    if (new Set(assetIds).size !== assetIds.length)
        return failure("invalid_input", "assetIds must be unique.");
    const size = assetIds.length;
    if (!Array.isArray(returnsInput) || returnsInput.length === 0)
        return failure("invalid_input", "returns must be a non-empty T by N matrix.", { assetIds });
    if (returnsInput.some((row) => !Array.isArray(row) || row.length !== size))
        return failure("invalid_input", "each return row must align to assetIds.", { assetIds });
    if ((returnsInput as unknown[][]).some((row) => row.some((value) => !finiteReal(value))))
        return failure("invalid_input", "returns must contain finite real numbers.", { assetIds });
    const returns = (returnsInput as unknown[][]).map((row) => row.slice() as number[]);
    const observations = returns.length;
    if (!Array.isArray(probabilitiesInput) || probabilitiesInput.length !== observations)
        return failure("invalid_input", "probabilities must align to return rows.", { assetIds });
    if ((probabilitiesInput as unknown[]).some((value) => !finiteReal(value) || value < 0))
        return failure("invalid_input", "probabilities must be finite and non-negative.", { assetIds });
    const probabilities = (probabilitiesInput as unknown[]).slice() as number[];
    const probabilityResidual = Math.abs(probabilities.reduce((sum, value) => sum + value, 0) - 1);
    if (probabilityResidual > TOL_PROBABILITY)
        return failure("invalid_input", "probabilities must sum to one; no silent normalization is applied.", { assetIds, probabilityResidual });
    if (!finiteReal(betaInput) || betaInput <= 0 || betaInput >= 1)
        return failure("invalid_input", "beta must be a confidence level strictly between zero and one.", { assetIds });
    const beta = betaInput;
    if (!finiteReal(gammaInput) || gammaInput <= 0)
        return failure("invalid_input", "gamma must be strictly positive in this bounded LP variant.", { assetIds });
    const gamma = gammaInput;
    const maxIterations = options.maxIterations ?? DEFAULT_MAX_ITERATIONS;
    if (!Number.isInteger(maxIterations) || (maxIterations as number) < 0 || (maxIterations as number) > DEFAULT_MAX_ITERATIONS)
        return failure("invalid_input", "maxIterations must be an integer from 0 through 10,000.", { assetIds });
    const iterationBudget = maxIterations as number;
    if (iterationBudget === 0)
        return failure("numerical_issue", "maxIterations=0 cannot produce an LP certificate.", { assetIds, maxIterations: iterationBudget });
    const variables = size + 1 + observations;
    if (variables > MAX_VERTEX_VARIABLES)
        return failure("numerical_issue", "the deterministic vertex core is bounded at 18 variables.", { assetIds, maxIterations: iterationBudget }, { variables, maxVertexVariables: MAX_VERTEX_VARIABLES });
    const meanReturns = Array.from({ length: size }, (_value, index) => returns.reduce((sum, row, scenario) => sum + probabilities[scenario] * row[index], 0));
    const inequalities: Array<{
        a: number[];
        b: number;
    }> = [];
    for (let index = 0; index < size; index += 1) {
        const lower = Array(variables).fill(0);
        lower[index] = -1;
        inequalities.push({ a: lower, b: 0 });
        const upper = Array(variables).fill(0);
        upper[index] = 1;
        inequalities.push({ a: upper, b: 1 });
    }
    for (let scenario = 0; scenario < observations; scenario += 1) {
        const lower = Array(variables).fill(0);
        lower[size + 1 + scenario] = -1;
        inequalities.push({ a: lower, b: 0 });
        const tail = Array(variables).fill(0);
        for (let index = 0; index < size; index += 1)
            tail[index] = -returns[scenario][index];
        tail[size] = -1;
        tail[size + 1 + scenario] = -1;
        inequalities.push({ a: tail, b: 0 });
    }
    const equality = Array(variables).fill(0);
    for (let index = 0; index < size; index += 1)
        equality[index] = 1;
    const objective = Array(variables).fill(0);
    meanReturns.forEach((value, index) => { objective[index] = value / Math.max(1, gamma); });
    objective[size] = -gamma / Math.max(1, gamma);
    probabilities.forEach((value, scenario) => { objective[size + 1 + scenario] = -(gamma / Math.max(1, gamma)) * value / (1 - beta); });
    const activeCount = variables - 1;
    const combinationCount = activeCount <= inequalities.length ? chooseCount(inequalities.length, activeCount) : 0;
    if (combinationCount > MAX_VERTEX_CANDIDATES)
        return failure("numerical_issue", "the LP vertex count exceeds the deterministic candidate limit.", { assetIds, maxIterations: iterationBudget }, { verticesRequired: combinationCount, maxVertexCandidates: MAX_VERTEX_CANDIDATES });
    const budget = Math.min(iterationBudget, combinationCount);
    const activeSets = combinations(inequalities.length, activeCount);
    let bestValue = -Infinity;
    let bestSolution: number[] | null = null;
    let checked = 0;
    for (const active of activeSets) {
        if (checked >= budget)
            break;
        checked += 1;
        const candidate = solveLinear([equality, ...active.map((index) => inequalities[index].a)], [1, ...active.map((index) => inequalities[index].b)]);
        if (candidate === null || candidate.some(x => !Number.isFinite(x)) || inequalities.some((constraint) => dot(constraint.a, candidate) > constraint.b + TOL_FEASIBILITY))
            continue;
        const value = dot(objective, candidate);
        if (bestSolution === null || value > bestValue + TOL_OBJECTIVE * Math.max(Math.abs(value), Math.abs(bestValue), 1e-300)) {
            bestValue = value;
            bestSolution = candidate;
        }
    }
    if (checked < combinationCount)
        return failure("numerical_issue", "maxIterations exhausted before all LP vertices were checked.", { assetIds, maxIterations: iterationBudget, iterations: checked }, { verticesChecked: checked, verticesRequired: combinationCount });
    if (bestSolution === null)
        return failure("numerical_issue", "no feasible LP vertex was found.", { assetIds, maxIterations: iterationBudget, iterations: checked });
    const weights = bestSolution.slice(0, size);
    const thresholdFromLp = bestSolution[size];
    const scenarioLosses = returns.map(row => -dot(row, weights));
    const expectedReturn = dot(meanReturns, weights);
    const { threshold, cvar } = weightedVarCvar(scenarioLosses, probabilities, beta);
    const tailSlacks = scenarioLosses.map(loss => Math.max(0, loss - threshold));
    const objectiveValue = expectedReturn - gamma * cvar;
    const lpObjective = bestValue * Math.max(1, gamma);
    const maxTailResidual = Math.max(0, ...scenarioLosses.map((loss, t) => loss - thresholdFromLp - bestSolution![size + 1 + t]));
    const lowerResidual = Math.max(0, -Math.min(...weights));
    const budgetResidual = Math.abs(weights.reduce((a, b) => a + b, 0) - 1);
    const primalResidual = Math.max(maxTailResidual, lowerResidual, budgetResidual);
    const objectiveScale = Math.max(Math.abs(lpObjective), Math.abs(objectiveValue), Math.abs(expectedReturn), Math.abs(gamma * cvar), 1e-15);
    const objectiveGap = Math.abs(lpObjective - objectiveValue);
    if (![lpObjective, objectiveValue, cvar, objectiveGap].every(Number.isFinite) || primalResidual > TOL_FEASIBILITY || objectiveGap > TOL_OBJECTIVE * objectiveScale)
        return failure('numerical_issue', 'the selected LP vertex failed the independent primal/objective audit.', { assetIds, maxIterations: iterationBudget, iterations: checked }, { verticesChecked: checked, verticesRequired: combinationCount, primalResidual, objectiveGap, objectiveScale, lpObjective, directObjective: objectiveValue });
    const tiedWeights: number[][] = [];
    for (const active of activeSets) {
        const candidate = solveLinear([equality, ...active.map(i => inequalities[i].a)], [1, ...active.map(i => inequalities[i].b)]);
        if (candidate === null || candidate.some(x => !Number.isFinite(x)) || inequalities.some(c => dot(c.a, candidate) > c.b + TOL_FEASIBILITY))
            continue;
        const ws = candidate.slice(0, size);
        const losses = returns.map(row => -dot(row, ws));
        const value = dot(meanReturns, ws) - gamma * weightedVarCvar(losses, probabilities, beta).cvar;
        if (Math.abs(value - objectiveValue) <= TOL_OBJECTIVE * Math.max(Math.abs(value), Math.abs(objectiveValue), 1e-300) && !tiedWeights.some(prior => Math.max(...ws.map((w, j) => Math.abs(w - prior[j]))) <= 1e-9))
            tiedWeights.push(ws);
    }
    return { status: 'optimal', variant: 'mean-cvar-long-only', method: 'lp-vertex-enumeration', branch: 'positive_gamma_lp', assetIds, weights, expectedReturn, cvar, thresholdZ: threshold, tailSlacks, scenarioLosses, objective: objectiveValue, budgetResidual, lowerBoundResidual: lowerResidual, probabilityResidual, iterations: checked, maxIterations: iterationBudget, solutionClass: tiedWeights.length > 1 ? 'non_unique' : 'unique', diagnostics: { beta, gamma, meanReturns, tailMass: 1 - beta, lpObjective, objectiveGap, objectiveScale, directObjective: objectiveValue, primalResidual, certificate: 'lp-vertex-enumeration', verticesChecked: checked, verticesRequired: combinationCount, thresholdFromLp, quantileFromLosses: threshold, distinctOptimalWeights: tiedWeights.length }, warnings: [] };
}
export const meanCvarOptimization = meanCvar;
export function meanCvar(assetIdsInput: unknown, returnsInput: unknown, probabilitiesInput: unknown, betaInput: unknown, gammaInput: unknown, options: {
    maxIterations?: unknown;
} = {}): MeanCvarResult {
    let scale = 1;
    let rows = returnsInput;
    if (Array.isArray(rows) && rows.length && rows.every(row => Array.isArray(row) && row.length && row.every(x => typeof x === 'number' && Number.isFinite(x)))) {
        scale = rows.reduce((m: number, row: number[]) => row.reduce((n, x) => Math.max(n, Math.abs(x)), m), 0) || 1;
        rows = rows.map((row: number[]) => row.map(x => x / scale));
    }
    const result = meanCvarNormalized(assetIdsInput, rows, probabilitiesInput, betaInput, gammaInput, options);
    const restore = (value: unknown): unknown => Array.isArray(value) ? value.map(restore) : typeof value === 'number' ? value * scale : value;
    const record = result as unknown as Record<string, unknown>;
    for (const key of ["expectedReturn", "cvar", "thresholdZ", "tailSlacks", "scenarioLosses", "objective"])
        if (record[key] !== null && record[key] !== undefined)
            record[key] = restore(record[key]);
    for (const key of ["meanReturns", "lpObjective", "objectiveGap", "objectiveScale", "directObjective", "thresholdFromLp", "quantileFromLosses"])
        if (key in result.diagnostics)
            result.diagnostics[key] = restore(result.diagnostics[key]);
    result.diagnostics.returnScale = scale;
    const allFinite = (value: unknown): boolean => typeof value === 'number' ? Number.isFinite(value) : Array.isArray(value) ? value.every(allFinite) : value !== null && typeof value === 'object' ? Object.values(value).every(allFinite) : true;
    if (result.status === 'optimal' && !allFinite(result))
        return failure('numerical_issue', 'output arithmetic exceeded finite return units.', { assetIds: assetIdsInput as string[] });
    return result;
}
function chooseCount(n: number, k: number): number { let v = 1; for (let i = 1; i <= k; i++)
    v = v * (n - k + i) / i; return Math.round(v); }