Hierarchical Equal Risk Contribution: Make the Cluster Cut Auditable

When I review a portfolio allocation, I do not begin with the final percentage beside each ticker. I begin with the evidence that created those percentages. Were the price histories adjusted consistently after splits? Did a cash dividend become a fake loss because the return series ignored the distribution? Was an old imported mark even possible for that instrument on that date? Were Egyptian pounds, Japanese yen, US dollars and Saudi riyals converted on one declared basis before somebody compared the holdings?

Those questions may sound separate from Hierarchical Equal Risk Contribution (HERC), but they are the reason I care about making HERC inspectable. A hierarchy can organize a valid covariance structure. It can also organize bad data into a very convincing-looking tree. The algorithm cannot tell which one you supplied.

This tutorial builds one reproducible HERC profile and keeps every important choice visible: supplied cluster count, correlation-profile distance, deterministic single linkage, inverse-variance weights inside terminal clusters, and variance-based risk balancing down the natural dendrogram. You will be able to calculate the result by hand, reproduce it in Python or TypeScript, and explain exactly why “equal risk” here does not mean equal final risk for every asset.

The useful question HERC answers

Suppose four holdings really represent two related pairs. A flat list shows four assets; a correlation tree shows two neighborhoods. HERC asks:

If I trust this covariance matrix, this distance rule, this hierarchy, and this terminal-cluster decision, how should capital flow through the tree using a declared risk measure?

That is narrower than portfolio optimization. HERC does not estimate expected returns, maximize a Sharpe ratio, validate security identity, choose a rebalancing date, or prove that the resulting portfolio will outperform. It is a structured risk-allocation heuristic.

Thomas Raffinot introduced HERC as a method that combines hierarchical clustering with modified top-down equal-risk allocation (original SSRN record). The accessible abstract also describes downside-risk extensions. This article does not silently absorb all of those variants. Its executable profile uses variance so that every input, split and output can be checked without a solver or external dependency.

HERC is a design, not one magic formula

Six choices define a hierarchical portfolio: the distance, linkage, stopping rule, split path, intra-cluster allocation, and inter-cluster allocation. Palomar's hierarchical portfolio chapter makes those dimensions explicit. That is valuable because libraries expose different combinations.

Decision	This tutorial	Other legitimate configurations
Risk measure	Variance	Volatility, conditional value at risk, conditional drawdown at risk
Clustering input	Euclidean distance between correlation-profile columns	Direct correlation distance or other dependence measures
Linkage	Deterministic single linkage	Ward, complete, average
Terminal clusters	Caller supplies `K`	Gap statistic or another cluster-selection rule
Inside a terminal cluster	Inverse variance	Equal weight, inverse volatility, optimized risk allocation
Between branches	Follow the dendrogram and allocate inversely to summed variance proxies	Ordered bisection, equal capital, post-processing constraints

The CRAN HierPortfolios documentation uses variance, defaults to Ward linkage and applies the Gap index when no cluster count is supplied. The official skfolio HERC documentation supports many risk measures and adds an explicit post-allocation constraint procedure. Those are useful implementations, but their defaults are not universal mathematical facts.

Here, K is required. That turns the stopping decision into a model input that can be tested, compared and recorded. It also avoids pretending that an unimplemented Gap-statistic simulation chose the clusters for us.

First build the relationship profiles

Let Sigma be a covariance matrix whose rows and columns follow the same ordered asset identifiers. Convert it to correlation:

\rho_{ij}=\frac{\Sigma_{ij}} {\sqrt{\Sigma_{ii}\Sigma_{jj}}}.

Then form the pair-profile matrix

D_{ij}=\sqrt{\frac{1-\rho_{ij}}{2}}.

The clustering input is not just the direct distance between assets i and j. Each column of D is an asset's relationship to the whole universe. We compare two columns:

d(i,j)=\left\lVert D_{\cdot i}-D_{\cdot j}\right\rVert_2 =\sqrt{\sum_k(D_{ki}-D_{kj})^2}.

That second distance transformation matters. Two assets can have a similar pairwise correlation but very different relationships to everything else. The implementation runs single linkage on these profile distances and uses a lexicographic rule for exact ties, so two runtimes do not invent different trees by traversing an unordered data structure differently.

Then make the stopping decision explicit

Cut the tree into exactly K terminal clusters. A terminal cluster is not merely a colored label in a chart. It is where recursive branch allocation stops and the intra-cluster rule takes over.

For terminal cluster C, compute inverse-variance weights:

u_i^{(C)}= \frac{1/\Sigma_{ii}}{\sum_{j\in C}1/\Sigma_{jj}}, \qquad i\in C.

Its risk proxy is the variance of that terminal portfolio:

v_C={u^{(C)}}^\mathsf{T}\Sigma_Cu^{(C)}.

The word “variance” is not interchangeable with “volatility” here. Replacing v_C with its square root changes every risk ratio. A library can legitimately define another risk measure; an article cannot change it mid-calculation and still call the result reproducible.

Equal risk happens at each branch

Now follow the natural shape of the dendrogram from the root. At a branch, let the left child contain terminal clusters L and the right child contain terminal clusters R. Sum their proxy risks:

v_L=\sum_{C\in L}v_C, \qquad v_R=\sum_{C\in R}v_C.

Split the parent capital inversely to those risks:

\alpha_L=\frac{v_R}{v_L+v_R}, \qquad \alpha_R=\frac{v_L}{v_L+v_R}.

The equality is local and exact:

\alpha_L v_L=\alpha_R v_R.

Multiply those factors down the tree until each selected terminal cluster is reached. Finally, multiply a cluster's capital by its internal weights.

This is the point that many short explanations lose: an unbalanced tree can equalize left and right proxy products at every split without making every terminal cluster's global proxy share equal. It certainly does not force every asset's final component-risk share to 1/N.

Rendering system map…

Work through a three-cluster portfolio

The example is synthetic by design. Every value can be redistributed and recalculated:

\Sigma=\begin{bmatrix} 0.04&0.036&0.004&0.002\\ 0.036&0.04&0.003&0.001\\ 0.004&0.003&0.01&0.008\\ 0.002&0.001&0.008&0.0225 \end{bmatrix}.

The hierarchy first groups A/B and C/D. Set K=3, leaving terminal clusters AB, C, and D.

Inside AB, the diagonal variances are equal, so the internal weights are 50% / 50%. Its variance is

v_{AB}=0.5^2(0.04)+2(0.5)(0.5)(0.036)+0.5^2(0.04)=0.038.

The singleton proxies are v_C=0.01 and v_D=0.0225. At the root, the left side is AB and the right side is C+D:

\alpha_{AB}=\frac{0.0325}{0.038+0.0325} =0.460992907801.

The right side receives 0.539007092199. Its child branch then assigns

\alpha_C=\frac{0.0225}{0.01+0.0225}=0.692307692308,

and alpha_D=0.307692307692. Multiplication down the tree produces:

Asset	Terminal cluster	Capital weight	Final component-risk share
A	A/B	23.0496%	35.6359%
B	A/B	23.0496%	34.6428%
C	C	37.3159%	19.9003%
D	D	16.5848%	9.8209%

Portfolio volatility is 0.111850951366. The component-risk shares sum to 100%, but they are not equal. That is not a failed HERC calculation. It is the expected consequence of balancing declared branch proxies instead of solving global asset-level ERC.

Three-cluster HERC tree cut, branch risks and final capital

Open the full-size allocation diagram. The red dashed line shows the model decision that changes everything below it: where the hierarchy stops.

Compare the stopping points yourself

The same covariance produces different portfolios when K changes:

K	Interpretation	Final weights A / B / C / D
1	One terminal cluster; internal inverse-variance allocation only	12.8571 / 12.8571 / 51.4286 / 22.8571%
2	Two terminal pairs	10.6881 / 10.6881 / 54.4319 / 24.1920%
3	AB remains grouped; C and D split	23.0496 / 23.0496 / 37.3159 / 16.5848%
4	Every asset is terminal; natural branches still control capital	14.4444 / 14.4444 / 49.2308 / 21.8803%

The guided playground lets you select all four states, step through the cut and branch logic, and try a singular covariance that correctly returns no portfolio. Its value is not animation for animation's sake. It makes K behave like the material model parameter it is.

How the implementation earns trust

The Python and TypeScript functions are pure: they do not mutate the input IDs or covariance, call a network provider, or choose a hidden cluster count. Each returns final weights plus the evidence needed to audit them:

correlation and distance matrices;
merge records and leaf order;
selected terminal clusters;
internal cluster weights and variance proxies;
every branch's left/right clusters, risks, allocations and balance residual;
portfolio variance, volatility and signed component-risk results; and
literal method labels for linkage, distance, stopping and risk allocation.

That output design matters more than a single clean weight vector. If a portfolio unexpectedly concentrates, I want to know whether the cause was the input covariance, the cluster cut, a terminal variance, or a branch multiplier. An audit trail turns “the model did it” into a calculation we can challenge.

The implementation normalizes covariance by its largest absolute value before factorization, then restores output risk units. It accepts only an exactly symmetric, strictly positive-definite covariance. It does not add a ridge, project to a positive-semidefinite matrix, or silently symmetrize. Those can be valid estimator policies, but they must remain visible upstream decisions.

Test more than the happy path

The shared fixture covers K=1,2,3,4, not just the portfolio used in the article. It also covers empty and duplicate identifiers, invalid cluster counts, singular and asymmetric covariance, negative variance, and nonnumeric input. Tests preserve three operational boundaries: corporate-action basis, stale or impossible prices, and quantity/FX basis.

Python independently reconstructs the canonical cluster variances and branch weights without calling implementation helpers. TypeScript starts the Python JSON-lines process and compares thirteen success/error cases directly. Both languages also test permutation identity, nonmutation, tiny positive variances, weight reconciliation, component-risk reconciliation, and branch-balance residuals.

Cross-language agreement is necessary, but it is not sufficient. Two programs can reproduce the same misconception. That is why the hand oracle begins with v_AB=0.038 and follows every split independently.

The portfolio-data trap is upstream

In real portfolio work, I have seen imported prices that were not historically possible after a stock split. I have seen a cash dividend appear as a loss when the calculation looked only at the price drop. I have seen quantity and price combinations that did not represent a real position, weights change because units changed, and holdings across several countries compared before local prices and FX were put onto one declared base-currency scale.

HERC will not protect you from any of that. If a split creates a fake return, the fake return changes covariance. If a market was closed while another was open, asynchronous observations can distort correlation. If local-currency and base-currency returns are mixed, the hierarchy answers a question no one intended to ask.

The production boundary should therefore be explicit:

identify the security and verify point-in-time price plausibility;
reconcile splits, cash distributions and other corporate actions;
separate quantities, local marks, FX rates and base-currency market values;
align sessions and return timestamps;
document missing-value, outlier and covariance-estimation policies; then
call HERC with the validated covariance and recorded cluster-count decision.

This is additional value the weight formula alone cannot provide: the tree is the end of an evidence pipeline, not the beginning.

HERC versus nearby methods

Method	Equality or allocation target	Uses a hierarchy?	Key distinction
Inverse-Volatility Weighting	Reciprocal standalone volatility scores	No	Ignores correlations during weight construction.
Equal Risk Contribution	Equal final asset component-risk shares	No	Solves a global asset-level condition.
Risk Budgeting	User-supplied final asset risk shares	No	Targets need not be equal.
Hierarchical Risk Parity	Inverse-variance branch allocation along an ordered bisection	Yes	No explicit terminal-cluster cut in the frozen HRP profile.
This HERC profile	Equal left/right proxy products at every natural dendrogram split	Yes	Explicit terminal clusters determine where recursion stops.

With K=2 on this balanced four-asset tree, HERC happens to match the related HRP worked example. That equality is a property of this tree and stopping point—not an identity between the methods. K=3 exposes the difference.

Where HERC can still fail

Single linkage can create a chaining tree. Ward or another linkage may be more suitable, but it is a different declared configuration.
An arbitrary K can be unstable. A Gap statistic can support selection, but its reference simulation and randomness must be reproducible.
Covariance estimates can change across windows, estimators and return bases.
Variance treats upside and downside symmetrically; downside-risk HERC is a separate contract.
Long-only positive weights do not enforce issuer, sector, country, liquidity, turnover, tax or transaction-cost constraints.
A stable calculation is not an empirical performance result.

What you should be able to explain now

HERC first discovers a relationship tree, then makes a stopping decision, builds a risk proxy for every terminal cluster, and balances those proxies as capital moves down natural branches. The final weights are products of all those decisions. They are not global ERC weights and they are not evidence that the input history was valid.

Before accepting an allocation, ask four questions: What return evidence built the covariance? What exact distance and linkage built the tree? Why was this cluster count chosen? Which risk measure and terminal allocation produced each branch proxy? If the result cannot answer those questions, it is not ready for a portfolio review.

The complete executable contract, full-precision fixture, code and references are in the canonical topic package. The next family moves from risk-only allocation to Black–Litterman, where expected-return views and their uncertainty enter the portfolio decision.

Workflows

Rendered from the canonical Mermaid sources linked by this article.

HERC computation and audit flow

Purpose: show where the terminal-cluster decision sits between tree discovery and capital allocation.

Rendering system map…

Takeaway: HERC does not ask for global asset-level ERC. It balances declared cluster risk proxies locally as it walks the selected tree.

References6 primary sources and evidence notes

Expand the source trail, evidence role, and limitations behind the engineering choices.

Accessed on 2026-09-15. The executable package uses a declared fixed-cluster, single-linkage, variance profile. Sources describing Ward linkage, automatic cluster selection, downside risk, or constrained libraries are evidence for variants, not hidden defaults.

R01 — The Hierarchical Equal Risk Contribution Portfolio

Organization or authors: Thomas Raffinot.
Source type: original working-paper record and accessible abstract.
Publication date: August 23, 2018; SSRN posting September 20, 2018.
Version: SSRN 3237540.
URL or DOI: SSRN record, DOI 10.2139/ssrn.3237540.
Accessed: 2026-09-15.
Jurisdiction: method research; not jurisdiction-specific.
Supports: HERC authorship; combination of hierarchical clustering with a modified top-down equal-risk allocation; extension to downside-risk measures.
Limitations: the full PDF was access-restricted from the available retrieval path. The package does not attribute paper-body formulas or empirical results that were not independently inspectable. The abstract's performance language is not repeated as a package claim.

R02 — Hierarchical Clustering-Based Portfolios, section 12.3.3

Organization or authors: Daniel P. Palomar, Portfolio Optimization.
Source type: open technical textbook chapter.
Publication/effective date: continuously published edition inspected 2026-09-15.
Version: online chapter 12.3.
URL: HERC section and hierarchical portfolio comparison.
Accessed: 2026-09-15.
Jurisdiction: mathematical education; not jurisdiction-specific.
Supports: the six explicit design choices of a hierarchical allocation; HERC's early-stopping, dendrogram-following, intra-cluster and inter-cluster distinctions; linkage sensitivity; warning against conclusions from one anecdotal backtest.
Limitations: secondary synthesis, used to explain and compare—not to claim that this package reproduces every Raffinot configuration.

R03 — HierPortfolios HERC_Portfolio documentation and source

Organization or authors: Carlos Trucios and Moon Jun Kwon; CRAN package HierPortfolios 1.0.2.
Source type: official maintained package documentation and GPL-2 source.
Publication/effective date: package date 2025-09-12; CRAN publication 2025-09-12 15:00:22 UTC.
Version: 1.0.2.
URL: CRAN manual, HERC documentation, maintained source.
Accessed: 2026-09-15.
Jurisdiction: software implementation; not jurisdiction-specific.
Supports: covariance input, covariance-to-correlation-to-profile-distance construction, selectable linkage, optional explicit cluster count, variance risk, inverse-variance terminal weights, dendrogram-recursive allocation.
Limitations: an implementation following its named references; defaults to Ward and Gap selection. This package instead declares single linkage and a required cluster count for deterministic dependency-free teaching.

R04 — skfolio HierarchicalEqualRiskContribution

Organization or authors: skfolio maintainers, documentation and BSD-3-Clause source.
Source type: official maintained library documentation/source.
Publication/effective date: current documentation inspected 2026-09-15.
Version: main documentation at access time.
URL: HERC estimator documentation, source.
Accessed: 2026-09-15.
Jurisdiction: software implementation; not jurisdiction-specific.
Supports: terminal-cluster inverse-risk allocation, top-down recursion that follows the dendrogram, distinction from HRP bisection, risk-measure variants, and the fact that post-allocation constraints are library choices.
Limitations: library defaults, cluster estimator, constraints, and broad risk-measure support are not universal HERC rules and are excluded here.

R05 — Estimating the Number of Clusters in a Data Set via the Gap Statistic

Organization or authors: Robert Tibshirani, Guenther Walther, Trevor Hastie.
Source type: original peer-reviewed paper.
Publication date: 2001.
Version: Journal of the Royal Statistical Society, Series B, 63(2), 411–423.
URL or DOI: author-hosted paper, DOI 10.1111/1467-9868.00293.
Accessed: 2026-09-15.
Jurisdiction: statistical method; not jurisdiction-specific.
Supports: the Gap statistic as an explicit cluster-count-selection method.
Limitations: this package requires clusterCount; it explains Gap selection as a production extension and does not implement or claim its result.

R06 — Building Diversified Portfolios that Outperform Out of Sample

Organization or authors: Marcos López de Prado.
Source type: original HRP paper.
Publication date: 2016.
Version: Journal of Portfolio Management 42(4), 59–69.
URL or DOI: accessible paper mirror, DOI 10.3905/jpm.2016.42.4.059.
Accessed: 2026-09-15.
Jurisdiction: method research; not jurisdiction-specific.
Supports: HRP comparison—single-linkage hierarchy, quasi-diagonalization and recursive bisection—so readers can see what HERC changes.
Limitations: paper title is not a performance promise from this package.

Publication and data boundary

All numeric inputs in this topic are synthetic and redistributable. No provider payload, customer holding, employer model, or live security value is included. The user's prior split, cash-dividend, impossible-price, quantity, currency and FX experiences motivate upstream validation checks but are not represented as dated historical facts.

Implementation

Full dependency-light reference implementations in both supported languages.

hierarchical-equal-risk-contribution.ts

/**
 * Pure D14-F02-A05 fixed-count variance HERC allocator.
 *
 * Frozen profile: deterministic single linkage on the paper's distance
 * profiles, an explicit terminal cluster count, inverse-variance weights
 * within each terminal cluster, and a dendrogram-following top-down variance
 * split. Other HERC variants remain explicit alternatives.
 */

export const WEIGHT_SUM_ABSOLUTE_TOLERANCE = 1e-12;
export const BRANCH_BALANCE_TOLERANCE = 1e-12;
export const PAIR_DISTANCE_TOLERANCE = 1e-12;

export type HierarchicalEqualRiskContributionErrorCode =
  | "invalid-input"
  | "invalid-covariance"
  | "numerical-range-invalid";

export interface HercBranchSplit {
  nodeAssets: number[];
  leftClusters: number[];
  rightClusters: number[];
  leftRisk: number;
  rightRisk: number;
  leftAllocation: number;
  rightAllocation: number;
  balanceResidual: number;
}

export interface HierarchicalEqualRiskContributionResult {
  assetIds: string[];
  weights: number[];
  componentRiskContributions: number[];
  componentRiskShares: number[];
  portfolioVariance: number;
  portfolioVolatility: number;
  sumWeights: number;
  quasiDiagonalOrder: number[];
  correlationMatrix: number[][];
  distanceMatrix: number[][];
  merges: Array<{ left: number[]; right: number[]; distance: number }>;
  clusters: number[][];
  clusterWeights: number[];
  withinClusterWeights: number[][];
  clusterRiskProxies: number[];
  clusterRiskProxyShares: number[];
  branchSplits: HercBranchSplit[];
  maxBranchBalanceResidual: number;
  method: "hierarchical-equal-risk-contribution";
  linkage: "single";
  distance: "euclidean-distance-between-sqrt((1-correlation)/2)-profiles";
  clusterCount: number;
  clusterRiskAllocation: "dendrogram-recursive-inverse-variance-risk";
  withinClusterAllocation: "inverse-variance";
  treeVariant: "fixed-count-single-linkage";
  covarianceScale: number;
  status: "ok";
}

export class HierarchicalEqualRiskContributionValidationError extends Error {
  readonly code: HierarchicalEqualRiskContributionErrorCode;
  constructor(code: HierarchicalEqualRiskContributionErrorCode, message: string) {
    super(message);
    this.name = "HierarchicalEqualRiskContributionValidationError";
    this.code = code;
  }
}

function invalid(code: HierarchicalEqualRiskContributionErrorCode, message: string): never {
  throw new HierarchicalEqualRiskContributionValidationError(code, message);
}

function sum(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(row: readonly number[], vector: readonly number[]): number {
  const products: number[] = [];
  for (let index = 0; index < vector.length; index += 1) products.push(row[index] * vector[index]);
  return sum(products);
}

function validateIds(assetIds: readonly unknown[]): string[] {
  if (!Array.isArray(assetIds)) invalid("invalid-input", "assetIds must be an ordered array.");
  if (assetIds.length === 0) invalid("invalid-input", "the asset universe must not be empty.");
  const ids: string[] = [];
  for (let index = 0; index < assetIds.length; index += 1) {
    const value = assetIds[index];
    if (typeof value !== "string" || value.trim() === "") invalid("invalid-input", `assetIds[${index}] must be a non-empty string.`);
    ids.push(value);
  }
  if (new Set(ids).size !== ids.length) invalid("invalid-input", "assetIds must be unique.");
  return ids;
}

function validateClusterCount(value: unknown, n: number): number {
  if (typeof value !== "number" || !Number.isInteger(value)) invalid("invalid-input", "clusterCount must be an integer.");
  if (value < 1 || value > n) invalid("invalid-input", "clusterCount must be between one and the number of assets.");
  return value;
}

function validateCovariance(covariance: readonly unknown[], n: number): { matrix: number[][]; scale: number; normalized: number[][] } {
  if (!Array.isArray(covariance)) invalid("invalid-input", "covariance must be a square numeric matrix.");
  if (covariance.length !== n) invalid("invalid-covariance", "covariance must have one row per asset.");
  const matrix: number[][] = [];
  for (let i = 0; i < covariance.length; i += 1) {
    const rawRow = covariance[i];
    if (!Array.isArray(rawRow)) invalid("invalid-covariance", `covariance[${i}] must be a numeric row.`);
    if (rawRow.length !== n) invalid("invalid-covariance", "covariance must be square.");
    const row: number[] = [];
    for (let j = 0; j < rawRow.length; j += 1) {
      const value = rawRow[j];
      if (typeof value !== "number") invalid("invalid-covariance", `covariance[${i}][${j}] must be numeric data.`);
      if (!Number.isFinite(value)) invalid("numerical-range-invalid", "covariance must contain only finite values.");
      row.push(value);
    }
    matrix.push(row);
  }
  for (let i = 0; i < n; i += 1) {
    for (let j = i + 1; j < n; j += 1) if (matrix[i][j] !== matrix[j][i]) invalid("invalid-covariance", "covariance must be exactly symmetric; no repair is performed.");
    if (!(matrix[i][i] > 0)) invalid("invalid-covariance", "covariance diagonal must be strictly positive.");
  }
  let scale = 0;
  for (const row of matrix) for (const value of row) if (Math.abs(value) > scale) scale = Math.abs(value);
  if (!Number.isFinite(scale) || scale <= 0) invalid("numerical-range-invalid", "covariance scale is zero or non-finite.");
  const normalized = matrix.map((row) => row.map((value) => value / scale));
  const lower = Array.from({ length: n }, () => Array<number>(n).fill(0));
  for (let i = 0; i < n; i += 1) {
    for (let j = 0; j <= i; j += 1) {
      let remainder = normalized[i][j];
      if (j > 0) remainder -= sum(Array.from({ length: j }, (_, k) => lower[i][k] * lower[j][k]));
      if (i === j) {
        if (!Number.isFinite(remainder) || remainder <= 0) invalid("invalid-covariance", "HERC requires a strictly positive-definite covariance.");
        lower[i][j] = Math.sqrt(remainder);
      } else {
        if (!(lower[j][j] > 0)) invalid("invalid-covariance", "covariance Cholesky factor is degenerate.");
        lower[i][j] = remainder / lower[j][j];
      }
    }
  }
  return { matrix, scale, normalized };
}

function correlationDistance(normalized: readonly number[][]): { correlation: number[][]; distance: number[][] } {
  const n = normalized.length;
  const correlation = Array.from({ length: n }, () => Array<number>(n).fill(0));
  const correlationDistance = Array.from({ length: n }, () => Array<number>(n).fill(0));
  for (let i = 0; i < n; i += 1) {
    correlation[i][i] = 1;
    for (let j = i + 1; j < n; j += 1) {
      // Preserve representable correlations when the diagonal product would
      // underflow; sqrt(diag_i) * sqrt(diag_j) is the stable denominator.
      const denominator = Math.sqrt(normalized[i][i]) * Math.sqrt(normalized[j][j]);
      if (!Number.isFinite(denominator) || denominator <= 0) invalid("numerical-range-invalid", "correlation denominator is outside the finite positive range.");
      let rho = normalized[i][j] / denominator;
      if (!Number.isFinite(rho) || rho < -1 - PAIR_DISTANCE_TOLERANCE || rho > 1 + PAIR_DISTANCE_TOLERANCE) invalid("invalid-covariance", "correlation is outside [-1, 1].");
      rho = Math.min(1, Math.max(-1, rho));
      const d = Math.sqrt((1 - rho) / 2);
      correlation[i][j] = rho; correlation[j][i] = rho; correlationDistance[i][j] = d; correlationDistance[j][i] = d;
    }
  }
  const distance = Array.from({ length: n }, () => Array<number>(n).fill(0));
  for (let i = 0; i < n; i += 1) {
    for (let j = i + 1; j < n; j += 1) {
      const value = Math.sqrt(sum(Array.from({ length: n }, (_, k) => (correlationDistance[k][i] - correlationDistance[k][j]) ** 2)));
      if (!Number.isFinite(value)) invalid("numerical-range-invalid", "profile distance is not finite.");
      distance[i][j] = value; distance[j][i] = value;
    }
  }
  return { correlation, distance };
}

type Node = { leaf: number } | { merge: [Node, Node] };
type Cluster = number[];
function compareClusters(left: readonly number[], right: readonly number[]): number {
  const length = Math.min(left.length, right.length);
  for (let i = 0; i < length; i += 1) if (left[i] !== right[i]) return left[i] - right[i];
  return left.length - right.length;
}

function comparePairs(
  left: readonly number[],
  right: readonly number[],
  otherLeft: readonly number[],
  otherRight: readonly number[],
): number {
  const leftComparison = compareClusters(left, otherLeft);
  return leftComparison !== 0 ? leftComparison : compareClusters(right, otherRight);
}
function singleLinkage(distance: readonly number[][]): { order: number[]; merges: Array<{ left: number[]; right: number[]; distance: number }>; root: Node } {
  let clusters: Cluster[] = distance.map((_, i) => [i]);
  let nodes: Node[] = distance.map((_, i) => ({ leaf: i }));
  const merges: Array<{ left: number[]; right: number[]; distance: number }> = [];
  while (clusters.length > 1) {
    let bestDistance = Number.POSITIVE_INFINITY; let bestLeft = 0; let bestRight = 1;
    for (let left = 0; left < clusters.length; left += 1) {
      for (let right = left + 1; right < clusters.length; right += 1) {
        let d = Number.POSITIVE_INFINITY;
        for (const i of clusters[left]) for (const j of clusters[right]) if (distance[i][j] < d) d = distance[i][j];
        const currentLeft = clusters[bestLeft]; const currentRight = clusters[bestRight];
        if (d < bestDistance || (d === bestDistance && comparePairs(clusters[left], clusters[right], currentLeft, currentRight) < 0)) {
          bestDistance = d; bestLeft = left; bestRight = right;
        }
      }
    }
    let leftCluster = clusters[bestLeft]; let rightCluster = clusters[bestRight];
    if (compareClusters(leftCluster, rightCluster) > 0) {
      const swap = leftCluster; leftCluster = rightCluster; rightCluster = swap;
      const indexSwap = bestLeft; bestLeft = bestRight; bestRight = indexSwap;
    }
    const merged = [...leftCluster, ...rightCluster].sort((a, b) => a - b);
    merges.push({ left: [...leftCluster], right: [...rightCluster], distance: bestDistance });
    const newClusters: Cluster[] = []; const newNodes: Node[] = [];
    for (let index = 0; index < clusters.length; index += 1) if (index !== bestLeft && index !== bestRight) { newClusters.push(clusters[index]); newNodes.push(nodes[index]); }
    newClusters.push(merged); newNodes.push({ merge: [nodes[bestLeft], nodes[bestRight]] });
    const pairs = newClusters.map((cluster, index) => ({ cluster, node: newNodes[index] }));
    pairs.sort((a, b) => compareClusters(a.cluster, b.cluster));
    clusters = pairs.map((pair) => pair.cluster); nodes = pairs.map((pair) => pair.node);
  }
  function leaves(node: Node): number[] {
    return "leaf" in node ? [node.leaf] : [...leaves(node.merge[0]), ...leaves(node.merge[1])];
  }
  return { order: leaves(nodes[0]), merges, root: nodes[0] };
}

function cutClusters(merges: readonly { left: number[]; right: number[]; distance: number }[], n: number, count: number): number[][] {
  let clusters: Cluster[] = Array.from({ length: n }, (_, i) => [i]);
  for (const merge of merges.slice(0, n - count)) {
    const leftIndex = clusters.findIndex((cluster) => compareClusters(cluster, merge.left) === 0);
    const rightIndex = clusters.findIndex((cluster) => compareClusters(cluster, merge.right) === 0);
    if (leftIndex < 0 || rightIndex < 0 || leftIndex === rightIndex) invalid("invalid-covariance", "hierarchy cut does not match merge records.");
    const left = clusters[leftIndex]; const right = clusters[rightIndex];
    clusters = clusters.filter((_, index) => index !== leftIndex && index !== rightIndex);
    clusters.push([...left, ...right].sort((a, b) => a - b));
    clusters.sort(compareClusters);
  }
  return clusters.map((cluster) => [...cluster]);
}

function nodeLeaves(node: Node): number[] {
  return "leaf" in node ? [node.leaf] : [...nodeLeaves(node.merge[0]), ...nodeLeaves(node.merge[1])];
}

function inverseVarianceWeights(covariance: readonly number[][]): { weights: number[]; variance: number } {
  const variances = covariance.map((row, index) => row[index]);
  if (variances.some((value) => !Number.isFinite(value) || value <= 0)) {
    invalid("numerical-range-invalid", "terminal-cluster variance is outside the finite positive range.");
  }
  const minimum = Math.min(...variances);
  const scores = variances.map((value) => minimum / value);
  const total = sum(scores);
  if (!Number.isFinite(total) || total <= 0 || scores.some((value) => !Number.isFinite(value) || value <= 0)) {
    invalid("numerical-range-invalid", "terminal-cluster inverse-variance normalization is invalid.");
  }
  const weights = scores.map((value) => value / total);
  const product = covariance.map((row) => dot(row, weights));
  const variance = dot(weights, product);
  if (!Number.isFinite(variance) || variance <= 0) {
    invalid("numerical-range-invalid", "terminal-cluster variance is not positive and finite.");
  }
  return { weights, variance };
}

function hierarchicalAllocation(
  covariance: readonly number[][],
  clusters: readonly number[][],
  root: Node,
): {
  withinClusterWeights: number[][];
  clusterRiskProxies: number[];
  clusterWeights: number[];
  clusterRiskProxyShares: number[];
  branchSplits: HercBranchSplit[];
  maxBranchBalanceResidual: number;
} {
  const withinClusterWeights: number[][] = [];
  const clusterRiskProxies: number[] = [];
  for (const cluster of clusters) {
    const local = cluster.map((i) => cluster.map((j) => covariance[i][j]));
    const solved = inverseVarianceWeights(local);
    withinClusterWeights.push(solved.weights);
    clusterRiskProxies.push(solved.variance);
  }

  const clusterWeights = Array<number>(clusters.length).fill(1);
  const clusterSets = clusters.map((cluster) => new Set(cluster));
  const branchSplits: HercBranchSplit[] = [];
  const subset = (left: Set<number>, right: Set<number>): boolean => [...left].every((value) => right.has(value));
  const equalSet = (left: Set<number>, right: Set<number>): boolean => left.size === right.size && subset(left, right);

  const recurse = (node: Node): void => {
    const nodeAssets = new Set(nodeLeaves(node));
    if (clusterSets.some((cluster) => equalSet(nodeAssets, cluster))) return;
    if ("leaf" in node) invalid("invalid-covariance", "dendrogram stopped before a terminal cluster was reached.");
    const leftAssets = new Set(nodeLeaves(node.merge[0]));
    const rightAssets = new Set(nodeLeaves(node.merge[1]));
    const leftClusters = clusterSets.map((cluster, index) => subset(cluster, leftAssets) ? index : -1).filter((index) => index >= 0);
    const rightClusters = clusterSets.map((cluster, index) => subset(cluster, rightAssets) ? index : -1).filter((index) => index >= 0);
    if (leftClusters.length === 0 || rightClusters.length === 0) invalid("invalid-covariance", "dendrogram cut does not define two terminal-cluster sides.");
    const leftRisk = sum(leftClusters.map((index) => clusterRiskProxies[index]));
    const rightRisk = sum(rightClusters.map((index) => clusterRiskProxies[index]));
    const totalRisk = leftRisk + rightRisk;
    if (!Number.isFinite(totalRisk) || totalRisk <= 0) invalid("numerical-range-invalid", "dendrogram branch risk is not positive and finite.");
    const leftAllocation = rightRisk / totalRisk;
    const rightAllocation = leftRisk / totalRisk;
    for (const index of leftClusters) clusterWeights[index] *= leftAllocation;
    for (const index of rightClusters) clusterWeights[index] *= rightAllocation;
    const leftBalancedRisk = leftAllocation * leftRisk;
    const rightBalancedRisk = rightAllocation * rightRisk;
    const balanceResidual = Math.abs(leftBalancedRisk - rightBalancedRisk) / (leftBalancedRisk + rightBalancedRisk);
    branchSplits.push({nodeAssets:[...nodeAssets].sort((a,b)=>a-b),leftClusters,rightClusters,leftRisk,rightRisk,leftAllocation,rightAllocation,balanceResidual});
    recurse(node.merge[0]);
    recurse(node.merge[1]);
  };

  recurse(root);
  const total = sum(clusterWeights);
  if (!Number.isFinite(total) || total <= 0 || Math.abs(total - 1) > WEIGHT_SUM_ABSOLUTE_TOLERANCE) {
    invalid("numerical-range-invalid", "dendrogram cluster weights do not satisfy the sum tolerance.");
  }
  const proxyContributions = clusterWeights.map((weight, index) => weight * clusterRiskProxies[index]);
  const proxyTotal = sum(proxyContributions);
  if (!Number.isFinite(proxyTotal) || proxyTotal <= 0) invalid("numerical-range-invalid", "dendrogram cluster risk total is not representable.");
  const clusterRiskProxyShares = proxyContributions.map((value) => value / proxyTotal);
  const maxBranchBalanceResidual = branchSplits.length === 0 ? 0 : Math.max(...branchSplits.map((split) => split.balanceResidual));
  if (maxBranchBalanceResidual > BRANCH_BALANCE_TOLERANCE) invalid("numerical-range-invalid", "dendrogram branch risk balance exceeded tolerance.");
  return { withinClusterWeights, clusterRiskProxies, clusterWeights, clusterRiskProxyShares, branchSplits, maxBranchBalanceResidual };
}

export function hierarchicalEqualRiskContributionWeights(
  assetIds: readonly unknown[],
  covariance: readonly unknown[],
  clusterCount: unknown,
): HierarchicalEqualRiskContributionResult {
  const ids = validateIds(assetIds);
  const count = validateClusterCount(clusterCount, ids.length);
  const validated = validateCovariance(covariance, ids.length);
  const { correlation, distance } = correlationDistance(validated.normalized);
  const hierarchy = singleLinkage(distance);
  const clusters = cutClusters(hierarchy.merges, ids.length, count);
  const allocation = hierarchicalAllocation(validated.normalized, clusters, hierarchy.root);
  const weights = Array<number>(ids.length).fill(0);
  clusters.forEach((cluster, clusterIndex) => cluster.forEach((assetIndex, localIndex) => { weights[assetIndex] = allocation.clusterWeights[clusterIndex] * allocation.withinClusterWeights[clusterIndex][localIndex]; }));
  const sumWeights = sum(weights);
  if (!Number.isFinite(sumWeights) || Math.abs(sumWeights - 1) > WEIGHT_SUM_ABSOLUTE_TOLERANCE) invalid("numerical-range-invalid", "HERC weights do not satisfy the sum tolerance.");
  const product = validated.normalized.map((row) => dot(row, weights));
  const variance = dot(weights, product);
  if (!Number.isFinite(variance) || variance <= 0) invalid("numerical-range-invalid", "HERC portfolio variance is not representable.");
  const scaleRoot = Math.sqrt(validated.scale); const portfolioVariance = variance * validated.scale; const portfolioVolatility = Math.sqrt(variance) * scaleRoot;
  if (![scaleRoot, portfolioVariance, portfolioVolatility].every(Number.isFinite) || portfolioVariance <= 0 || portfolioVolatility <= 0) invalid("numerical-range-invalid", "restoring HERC risk scale lost positive risk.");
  const normalizedVolatility = Math.sqrt(variance);
  // Divide before multiplying by a tiny weight so a finite contribution is
  // not lost to an avoidable product underflow.
  const componentRiskContributions = weights.map((weight, index) => weight * (product[index] / normalizedVolatility) * scaleRoot);
  const componentRiskShares = weights.map((weight, index) => weight * (product[index] / variance));
  if (componentRiskContributions.some((value) => !Number.isFinite(value)) || componentRiskShares.some((value) => !Number.isFinite(value))) invalid("numerical-range-invalid", "HERC component risk is non-finite.");
  const restoredClusterRiskProxies = allocation.clusterRiskProxies.map((value) => value * validated.scale);
  const restoredBranchSplits = allocation.branchSplits.map((split) => ({
    ...split,
    leftRisk: split.leftRisk * validated.scale,
    rightRisk: split.rightRisk * validated.scale,
  }));
  return {
    assetIds: ids, weights, componentRiskContributions, componentRiskShares,
    portfolioVariance, portfolioVolatility, sumWeights,
    quasiDiagonalOrder: hierarchy.order, correlationMatrix: correlation, distanceMatrix: distance, merges: hierarchy.merges,
    clusters, clusterWeights: allocation.clusterWeights, withinClusterWeights: allocation.withinClusterWeights,
    clusterRiskProxies: restoredClusterRiskProxies, clusterRiskProxyShares: allocation.clusterRiskProxyShares,
    branchSplits: restoredBranchSplits, maxBranchBalanceResidual: allocation.maxBranchBalanceResidual,
    method: "hierarchical-equal-risk-contribution", linkage: "single",
    distance: "euclidean-distance-between-sqrt((1-correlation)/2)-profiles", clusterCount: count,
    clusterRiskAllocation: "dendrogram-recursive-inverse-variance-risk",
    withinClusterAllocation: "inverse-variance",
    treeVariant: "fixed-count-single-linkage", covarianceScale: validated.scale, status: "ok",
  };
}

Playground

The embedded lab now expands to its full document height, keeping the article as the only scroll surface.

Full-height labplaygroundOpen full screen

Written by

Islam Baraka

Fintech engineer building market-data and financial systems, and the author of every article, glossary record, and reference implementation on The Fintech Builder.