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.
A covariance matrix can describe co-movement without explaining where it came from. A factor model supplies an explicit structure: assets share exposure to common drivers, while each also carries residual risk.
The benefit is not merely a smaller matrix. You can ask a concrete counterfactual: what happens to portfolio risk if a factor becomes more volatile, an exposure changes, or one asset's specific variance rises?
This topic constructs covariance from supplied loadings, factor covariance and specific variances. It does not estimate factors or claim that a chosen factor set explains the real market.
Open this figure at full size.
State the model assumptions before multiplying
Write residual asset returns as r=Bf+ε, where B is an assets-by-factors loading matrix. Let F be factor covariance and D the diagonal matrix of specific variances. Assuming factors and residuals are uncorrelated and residual cross-covariances are zero,
The MOSEK portfolio-modeling reference describes this common factor-risk structure. The diagonal residual assumption is substantive: omitted shared risks can appear as residual correlation, in which case diagonal D is incomplete.
Loadings describe exposure in the units of the chosen factors. If factors are decimal returns, dimensionless loadings are common. A factor measured in percentage points or another economic unit requires compatible loadings and covariance units. Relabeling a factor without rescaling its exposures changes the model.
One factor and two assets
Let B=[[1],[2]], F=[[0.04]], and specific variances be[0.01,0.02]. The common-risk matrix is
Add specific variances only to the diagonal to obtain
The off-diagonal 0.08 is entirely common-factor covariance under this model. Increasing the first specific variance from 0.01 to 0.03 raises Σ_11 by 0.02 and leaves all other cells unchanged. Increasing F instead changes both diagonals and the off-diagonal.
Open this figure at full size.
Download the exact worked input and expected values.
Attribute a portfolio's variance
For portfolio weights a=[0.5,0.5], aggregate factor exposure is Bᵀa=1.5. Common-factor variance is 1.5²×0.04=0.09. Specific variance is 0.5²×0.01+0.5²×0.02=0.0075. Total variance is0.0975.
Multiplying aᵀΣa from the completed matrix gives the same answer. This two-route calculation is a strong consistency check and a useful explanation for a reader: a portfolio can diversify some specific risk while retaining a large shared factor exposure.
With correlated factors, factor-attribution terms include cross-products. You cannot always assign independent nonnegative “risk slices” to each factor simply by reading its diagonal variance. The total common component BFBᵀ remains well defined; a finer attribution requires a declared allocation convention.
A model can be numerically invalid before it is economically wrong
F must be square, symmetric and positive semidefinite. The loading column count must match F's dimension. The specific-variance vector must match the asset count and contain finite nonnegative values.
These conditions imply Σ is PSD: aᵀBFBᵀa≥0 and Σa_i²D_i≥0. Positive definiteness is not automatic. Zero specific variances and insufficient factor rank can leave singular directions. If a downstream calculation needs an inverse, test that stronger requirement explicitly.
| Intervention | Matrix effect | Economic interpretation within the model |
|---|---|---|
| Raise one D_i | One diagonal rises | More asset-specific risk |
| Raise a factor variance | Multiple cells can change | More shared driver uncertainty |
| Flip one loading's sign | Some covariances can reverse | Opposite exposure to that factor |
| Add residual correlations | Requires a non-diagonal residual model | Current diagonal-D assumption is insufficient |
Open this figure at full size.
The playground is a structural counterfactual
The canonical lab uses four assets and two factors with labeled synthetic loadings. Change an exposure, a factor variance scale or a specific-variance scale. Inspect the common matrix, diagonal addition and completed covariance separately.
Step follows the matrix-construction stages and highlights the actual contribution being added, rather than pretending that an unrelated return sample estimates the supplied factors. The zero-specific-risk edge case makes rank limitations visible. An indefinite F or negative specific variance is a genuine rejected input.
Open this figure at full size.
Open the standalone guided playground. The embedded playground and runnable code are available on this page. Download the factor-model teaching input.
What would be needed to use estimated factors
A statistical or economic factor workflow must explain factor selection, exposure estimation, data availability, rebalancing and residual diagnostics. A loading estimated after a historical date cannot be used as if it were known beforehand. Factor definitions and return units must remain stable or be versioned.
The runtime intentionally requires no irrelevant returns matrix. Its inputs already specify the model. Adding unused historical returns would imply evidence that the calculation does not actually consume.
The direct transparent multiplication costs O(p²k²) for p assets and k factors; efficient matrix implementations can reduce the work. Tests check exact common and total cells, factor dimensions, symmetry, semidefiniteness, negative specific risk and portfolio arithmetic.
Use this construction when risk explanation and structured scenarios matter. Compare sample covariance to see what the model imposes beyond the observed panel. A clear factor story is useful only while its assumptions and unexplained residual behavior remain visible.
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 factor_model_covariance import calculate
data = json.loads(Path("teaching-path.json").read_text())
result = calculate(data)
print(result["latest"])
import {calculate} from './factor_model_covariance.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.
B asset-by-factor, F symmetric PSD, D nonnegative asset-specific variances. Factors/residuals uncorrelated and residual covariance diagonal by model assumption. No irrelevant returns argument.
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.
Factor-Model Covariance — calculation-flow
Factor-Model Covariance — 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: MOSEK Portfolio Optimization Cookbook — Factor models — 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.
Scope of evidence
B asset-by-factor, F symmetric PSD, D nonnegative asset-specific variances. Factors/residuals uncorrelated and residual covariance diagonal by model assumption. No irrelevant returns argument.
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-A05 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,"factor_covariance");
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-A05",title:"Factor-Model Covariance",parameters:p,...result};
}
The embedded lab now expands to its full document height, keeping the article as the only scroll surface.