By Islam BarakaFor analysts and developers who know aligned return vectors, matrix multiplication, and covariance versus correlation. The goal is to reproduce the mechanism, inspect its failure states, and decide what the output can legitimately tell you—not to fit or endorse a trading strategy.
Two implementations can both say “OAS” and return different shrinkage intensities without either having a simple coding typo. The original finite-dimensional formula retains terms that a common software implementation deliberately omits.
This tutorial fixes the original finite-p OAS coefficient with ML empirical covariance. You will calculate it, identify the isotropic boundary, and learn how to make a meaningful comparison with Ledoit–Wolf or a library's simplified OAS.
Open this figure at full size.
The same target, a different intensity derivation
Let C be the centered covariance with denominator n, p the number of assets, T=tr(C), and Q=tr(C²). For a symmetric matrix, Q is the sum of squared entries. The target is μI with μ=T/p, and C_hat=(1−ρ)C+ρμI.
The coefficient used here is
The original OAS paper develops the shrinkage construction under a Gaussian covariance-estimation setting. “Oracle approximating” does not mean the runtime observes the true population covariance. It replaces unknown quantities using an approximation derived under that framework.
A hand-checkable two-asset case
Take n=4 and C=[[1,0.5],[0.5,1]]. Then p=2, T=2 and Q=2.5. Because 2/p=1, the Q term in the numerator vanishes. The numerator is 4; the denominator is(4+1−1)×(2.5−4/2)=4×0.5=2.
The raw ratio is 2, clipped to ρ=1. Since μ=1, the estimate is the identity matrix. The four centered synthetic observations producing this C are included in the executable fixture; the test does not simply feed a desired covariance into a different interface.
This endpoint is simple enough to verify by hand. The independent numerical suite also checks nontrivial partial-shrinkage cases using the original formula evaluated through NumPy matrix operations.
Open this figure at full size.
Download the exact worked input and expected values.
Why the library comparison needs a convention note
The scikit-learn OAS documentation explicitly states that its implementation omits the original 2/p terms as a simplification. Those terms can matter when p is small.
Before comparing numbers, record the formula, covariance denominator, centering rule and sample count. An n−1 covariance passed into an n-based derivation changes the output scale even when the intensity happens to remain unchanged. A different formula can change the intensity too.
A reproducible test should therefore compare this implementation with the original finite-p equation, not declare failure merely because a simplified OAS API differs. Conversely, the convention note is not an excuse to accept arbitrary discrepancies: the exact frozen expression is independently evaluated.
The zero-denominator boundary has meaning
Q−T²/p measures dispersion of covariance eigenvalues around their average. If C=μI, this term is zero: the empirical matrix already matches the target. There is nothing to change.
For a positive isotropic matrix, the package reportsρ=1 and returns the same matrix. For an all-zero matrix it reportsρ=0 and returns zeros. These explicit conventions avoid 0/0 and make diagnostics deterministic. A tiny denominator is treated relative to matrix scale, not with an arbitrary absolute threshold that behaves differently in percentages and decimals.
| Situation | Expected interpretation |
|---|---|
| Large empirical eigenvalue spread | Target is meaningfully different |
| Isotropic positive C | Target equals input; intensity is not economically identifiable |
| All-zero C | No observed variation; no inverse is promised |
| Strongly non-Gaussian data | Evaluate the derivation's applicability, not just code validity |
Open this figure at full size.
What the playground should explain
Step through the 64-row synthetic panel. Inspect C, μ, T, Q, numerator, denominator and the clipped coefficient. The matrix view shows the empirical estimate beside the shrunk result, so a change inρ has a visible numerical consequence.
Use the sample-prefix control to compare a shorter and longer endpoint sample. Both the matrix and coefficient can change. Try the zero or collinear edge case and inspect the explicit boundary state rather than expecting every covariance estimate to be invertible.
Open this figure at full size.
Open the standalone guided playground. The embedded playground and runnable code are available on this page. Download the 64-observation teaching input.
OAS is an estimator, not a performance certificate
OAS and Ledoit–Wolf use the same scaled-identity target here but different intensity constructions. Neither is automatically best for every asset universe, sample size, tail distribution or downstream objective. The Gaussian derivation is a reason to investigate assumptions, not to relabel every return panel Gaussian.
The reference accepts a complete finite aligned matrix, estimates endpoint means and uses O(np²) direct computation. Missing-data treatment, robust covariance, factor structure and time-varying covariance are separate choices. Keep asset units and ordering consistent.
For a real application, compare held-out covariance or portfolio-risk behavior under an explicit chronological protocol. Do not select the estimator using the final evaluation period and then present that same result as independent evidence.
The practical value of this topic is precision about conventions. You should be able to explain a disagreement down to a trace term or denominator, not stop at “both libraries implement OAS.”
Reproduce and inspect the calculation
The Python and TypeScript tabs contain standalone implementations, not imports into an unseen runtime. Both expose calculate(input_data). Feed the worked JSON's input object into that entry point. For the longer experiment, use the teaching-path JSON directly.
import json
from pathlib import Path
from oracle_approximating_shrinkage import calculate
data = json.loads(Path("teaching-path.json").read_text())
result = calculate(data)
print(result["latest"])
import {calculate} from './oracle_approximating_shrinkage.ts';
const result = calculate(inputData); // inputData is the downloaded JSON object
console.log(result.latest);
Place the downloaded input beside your script and the standalone source on its import path. The Python reference uses the standard library; the TypeScript reference has no external runtime dependency. Shared tests include independent numeric anchors, valid boundaries, rejected inputs and cross-language output comparisons. They establish arithmetic, not forecasting performance.
Evidence and scope
This article uses authored synthetic calculations and primary technical references, reviewed 2026-09-10. Historical market examples are deferred until identity, adjustment basis, chronology and redistribution rights can be verified. No personal trading history or search-ranking superiority is asserted.
Original finite-p OAS coefficient with ML covariance C; retains 2/p terms omitted by sklearn simplification. Isotropic positive covariance rho=1; zero matrix rho=0.
- scikit-learn — OAS API and implementation notes
- Chen et al. — Shrinkage Algorithms for MMSE Covariance Estimation
Continue the investigation
- Sample Covariance: compare its assumptions and information boundary before comparing the numbers.
- Graphical-Lasso Covariance: compare its assumptions and information boundary before comparing the numbers.
Rendered from the canonical Mermaid sources linked by this article.
Oracle Approximating Shrinkage — calculation-flow
Oracle Approximating Shrinkage — decision-boundary
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.
Reviewed 2026-09-10. Primary technical documentation and papers; synthetic arithmetic is author-derived. This is a targeted source review, not a verified review of Google's top ten results and not a claim of ranking superiority.
- S1: scikit-learn — OAS API and implementation notes — accessed 2026-09-10. Rolling official documentation snapshot; exact package versions used for numerical comparisons are recorded in the repair numeric-evidence.json. Supports the definition and declared convention, not investment performance. Jurisdiction: not applicable to this mathematical reference.
- S2: Chen et al. — Shrinkage Algorithms for MMSE Covariance Estimation — accessed 2026-09-10. Original research publication; not a current market observation. Supports the definition and declared convention, not investment performance. Jurisdiction: not applicable to this mathematical reference.
Scope of evidence
Original finite-p OAS coefficient with ML covariance C; retains 2/p terms omitted by sklearn simplification. Isotropic positive covariance rho=1; zero matrix rho=0.
Historical case: deferred. No public provider dataset, historical performance claim, or personal trading anecdote is used. Synthetic examples demonstrate arithmetic, not market efficacy. Sources are not copied as article prose.
Accessed: 2026-09-10.
Supports: estimator definition and the explicitly declared variants.
Limitations: technical documentation does not verify a real market feed, author experience, forecast efficacy or search-result superiority. Original-paper access limitations are recorded in the repair report.
Full dependency-light reference implementations in both supported languages.
/** Standalone D10-F04-A04 reference. Generated from validated D10 v2 source. */
export class ContractError extends Error {}
type RecordValue=Record<string, any>;
type Matrix=number[][];
const sum=(x:number[]):number=>x.reduce((a,b)=>a+b,0);
function requireValue(ok: unknown, code: string, message: string): asserts ok {
if (!ok) throw new ContractError(`${code}: ${message}`);
}
function finite(x: unknown, name: string): number {
requireValue(typeof x === 'number' && Number.isFinite(x), 'NUMBER', `${name} must be a finite number`);
return x;
}
function integer(x: unknown, name: string, minimum = 0, maximum = 10000): number {
const v = finite(x, name);
requireValue(Number.isInteger(v) && v >= minimum && v <= maximum, 'INTEGER', `${name} must be an integer in [${minimum}, ${maximum}]`);
return v;
}
function param(p: RecordValue, key: string, fallback: number): number {
return finite(Object.hasOwn(p, key) ? p[key] : fallback, key);
}
function option(p: RecordValue, key: string, fallback: any): any {
return Object.hasOwn(p, key) ? p[key] : fallback;
}
function positive(p: RecordValue, key: string, fallback: number): number {
const v = param(p, key, fallback);
requireValue(v > 0, 'RANGE', `${key} must be positive`);
return v;
}
function vector(value: unknown, name: string, minimum = 1): number[] {
requireValue(Array.isArray(value) && value.length >= minimum, 'SHAPE', `${name} needs ${minimum} or more values`);
return value.map((v, i) => finite(v, `${name}[${i}]`));
}
function matrix(value: unknown, name: string, minRows = 1, minCols = 1): Matrix {
requireValue(Array.isArray(value) && value.length >= minRows, 'SHAPE', `${name}: too few rows`);
const rows = value.map(row => vector(row, name, minCols));
requireValue(rows.every(row => row.length === rows[0].length), 'SHAPE', `${name} must be rectangular`);
return rows;
}
function eye(n: number, value = 1): Matrix {
return Array.from({length: n}, (_, i) => Array.from({length: n}, (_, j) => i === j ? value : 0));
}
function psd(a: Matrix, name: string, strict = false): number | null {
const n = a.length;
requireValue(a.every(row => row.length === n), 'SHAPE', `${name} must be square`);
const tolerance = 1e-12 * Math.max(...a.flat().map(Math.abs), 1e-300);
requireValue(a.every((row, i) => row.every((v, j) => Math.abs(v - a[j][i]) <= tolerance)), 'PSD', `${name} must be symmetric`);
const lower = eye(n), pivots: number[] = [];
for (let j = 0; j < n; j++) {
const pivot = a[j][j] - sum(pivots.map((v, k) => lower[j][k] ** 2 * v));
requireValue(strict ? pivot > tolerance : pivot >= -tolerance, 'PSD', `${name} must be positive semidefinite (strict when requested)`);
pivots.push(pivot > tolerance ? pivot : 0);
for (let i = j + 1; i < n; i++) {
const residual = a[i][j] - sum(pivots.slice(0, j).map((v, k) => lower[i][k] * lower[j][k] * v));
requireValue(pivots[j] !== 0 || Math.abs(residual) <= tolerance, 'PSD', `${name} has invalid zero pivot`);
lower[i][j] = pivots[j] ? residual / pivots[j] : 0;
}
}
return strict ? sum(pivots.map(Math.log)) : null;
}
function inverse(a: Matrix): Matrix {
const n = a.length, identity = eye(n), work = a.map((row, i) => [...row, ...identity[i]]);
for (let j = 0; j < n; j++) {
let k = j;
for (let i = j + 1; i < n; i++) if (Math.abs(work[i][j]) > Math.abs(work[k][j])) k = i;
[work[j], work[k]] = [work[k], work[j]];
requireValue(work[j][j] !== 0, 'NUMERIC', 'singular inverse');
const pivot = work[j][j];
work[j] = work[j].map(v => v / pivot);
for (let i = 0; i < n; i++) if (i !== j) {
const factor = work[i][j];
work[i] = work[i].map((v, l) => v - factor * work[j][l]);
}
}
return work.map(row => row.slice(n));
}
function glassoCovariance(S: Matrix, alpha: number, maxIter = 200, tol = 1e-8): [Matrix, RecordValue] {
const n = S.length;
requireValue(n >= 2 && n <= 12 && alpha > 0 && Math.min(...S.map((row, i) => row[i])) > 0, 'RANGE', 'glasso needs 2..12 nonconstant assets and positive alpha');
const scale = Math.max(...S.map((row, i) => row[i])), C = S.map(row => row.map(v => v / scale)), penalty = alpha / scale;
// Dual-feasible SPD initialization; a diagonal start is not generally feasible.
const maximumOff = Math.max(...C.flatMap((row, i) => row.filter((_, j) => i !== j).map(Math.abs)));
const blend = maximumOff ? Math.min(.1, penalty / maximumOff) : .1;
const W = C.map((row, i) => row.map((v, j) => i === j ? v : (1 - blend) * v)), history: RecordValue[] = [];
for (let iteration = 0; iteration < maxIter; iteration++) {
for (let j = 0; j < n; j++) {
const indices = Array.from({length: n}, (_, i) => i).filter(i => i !== j), beta = indices.map(() => 0);
for (let inner = 0; inner < 2000; inner++) {
let change = 0;
indices.forEach((i, k) => {
const partial = C[i][j] - sum(indices.map((q, l) => l === k ? 0 : W[i][q] * beta[l]));
const updated = Math.sign(partial) * Math.max(Math.abs(partial) - penalty, 0) / W[i][i];
change = Math.max(change, Math.abs(updated - beta[k])); beta[k] = updated;
});
if (change < tol * 1e-5) break;
}
const values = indices.map(i => sum(indices.map((q, l) => W[i][q] * beta[l])));
indices.forEach((i, k) => { W[i][j] = W[j][i] = values[k]; });
}
const logdet = psd(W, 'glasso covariance', true)!, precision = inverse(W);
let residual = 0, trace = 0, l1 = 0;
for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
const gradient = C[i][j] - W[i][j];
const trial = precision[i][j] - gradient;
const prox = Math.sign(trial) * Math.max(Math.abs(trial) - penalty, 0);
const error = i === j ? Math.abs(gradient) : Math.abs(precision[i][j] - prox);
residual = Math.max(residual, error); trace += C[i][j] * precision[j][i];
if (i !== j) l1 += penalty * Math.abs(precision[i][j]);
}
const gap = trace - n + l1;
history.push({iteration: iteration + 1, kkt_residual: residual, dual_gap: gap, objective: trace + logdet + l1 + n * Math.log(scale)});
if (residual <= tol && Math.abs(gap) <= tol * Math.max(1, n)) return [W.map(row => row.map(v => v * scale)),
{iterations: iteration + 1, converged: true, alpha, precision: precision.map(row => row.map(v => v / scale)), kkt_residual: residual, dual_gap: gap, history}];
}
throw new ContractError('CONVERGENCE: graphical lasso did not satisfy KKT and dual-gap tolerances; iteration budget exhausted');
}
function covariance(data: RecordValue, p: RecordValue, kind: string): RecordValue {
if (kind === 'factor_covariance') {
const B = matrix(data.loadings, 'loadings', 2), F = matrix(data.factor_covariance, 'factor_covariance'), D = vector(data.specific_variances, 'specific_variances');
const assets = B.length, factors = B[0].length;
requireValue(F.length === factors && D.length === assets, 'SHAPE', 'factor dimensions do not match');
psd(F, 'factor_covariance'); requireValue(Math.min(...D) >= 0, 'RANGE', 'specific variances must be nonnegative');
const common = B.map(row => B.map(other => sum(row.map((v, k) => sum(other.map((u, l) => v * F[k][l] * u))))));
const result = common.map((row, i) => row.map((v, j) => v + (i === j ? D[i] : 0)));
return {matrix: result, latest: result, ready: true, ready_at: 0, diagnostics: {common, specific_variances: D, assets, factors, causal: true}};
}
const X = matrix(data.returns, 'returns', 2, 2), n = X.length, dim = X[0].length;
const means = X[0].map((_, j) => sum(X.map(row => row[j])) / n), Z = X.map(row => row.map((v, j) => v - means[j]));
const C = X[0].map((_, i) => X[0].map((_, j) => sum(Z.map(z => z[i] * z[j])) / n));
const diagnostics: RecordValue = {observations: n, assets: dim, means, causal: true, ml_covariance: C};
let result: Matrix;
if (kind === 'sample_covariance') { result = C.map(row => row.map(v => v * n / (n - 1))); diagnostics.denominator = n - 1; }
else if (kind === 'ewma_covariance') {
const decay = param(p, 'decay', .94);
requireValue(decay > 0 && decay < 1, 'RANGE', 'decay must be in (0,1)');
result = Object.hasOwn(data, 'initial_covariance') ? matrix(data.initial_covariance, 'initial_covariance') : eye(dim, 0);
requireValue(result.length === dim, 'SHAPE', 'initial covariance dimensions mismatch'); psd(result, 'initial_covariance');
const states: Matrix[] = [];
for (const row of X) { result = result.map((old, i) => old.map((v, j) => decay * v + (1 - decay) * row[i] * row[j])); states.push(result); }
Object.assign(diagnostics, {states, mean_convention: 'supplied zero-mean residuals', weight_mass: 1 - decay ** n, seed_weight: decay ** n});
} else if (kind === 'graphical_lasso') {
const solver = glassoCovariance(C, positive(p, 'alpha', .00002), integer(option(p, 'max_iter', 200), 'max_iter', 1, 2000), positive(p, 'tolerance', 1e-8));
result = solver[0]; diagnostics.graphical_lasso = solver[1];
} else {
const trace = sum(C.map((row, i) => row[i])), mu = trace / dim, tr2 = sum(C.flat().map(v => v * v));
const delta = sum(C.flatMap((row, i) => row.map((v, j) => (v - (i === j ? mu : 0)) ** 2)));
let rho: number;
if (kind === 'ledoit_wolf') {
const noise = sum(Z.map(z => sum(C.flatMap((row, i) => row.map((v, j) => (z[i] * z[j] - v) ** 2))))) / n ** 2;
rho = delta > 0 ? Math.min(1, Math.max(0, noise / delta)) : 0;
Object.assign(diagnostics, {noise_estimate: noise, target_distance: delta});
} else {
const numerator = (1 - 2 / dim) * tr2 + trace * trace, denominator = (n + 1 - 2 / dim) * (tr2 - trace * trace / dim);
rho = tr2 === 0 ? 0 : denominator > 1e-14 * tr2 ? Math.min(1, Math.max(0, numerator / denominator)) : 1;
Object.assign(diagnostics, {numerator, denominator, variant: 'original finite-p OAS'});
}
result = C.map((row, i) => row.map((v, j) => (1 - rho) * v + (i === j ? rho * mu : 0)));
Object.assign(diagnostics, {shrinkage: rho, target: mu});
}
return {matrix: result, latest: result, ready: true, ready_at: 0, diagnostics};
}
export function calculate(data: RecordValue): RecordValue {
requireValue(data && typeof data === 'object' && !Array.isArray(data),'SHAPE','input must be an object');
const p=Object.hasOwn(data,'parameters')?data.parameters:{};
requireValue(p && typeof p === 'object' && !Array.isArray(p),'SHAPE','parameters must be an object');
const result=covariance(data,p,"oas");
function check(v:any):void {
if(typeof v==='number')requireValue(Number.isFinite(v),'NUMERIC','nonfinite computed output');
else if(Array.isArray(v))v.forEach(check);
else if(v && typeof v==='object')Object.values(v).forEach(check);
}
check(result);
return {topic_id:"D10-F04-A04",title:"Oracle Approximating Shrinkage",parameters:p,...result};
}
The embedded lab now expands to its full document height, keeping the article as the only scroll surface.