Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 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 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 | 1x 16x 16x 16x 55x 55x 55x 55x 16x 16x 16x 16x 16x 16x 16x 85x 85x 16x 16x 4x 4x 4x 4x 4x 24x 24x 24x 24x 77x 77x 24x 24x 24x 24x 6x 6x 6x 6x 6x 20x 20x 6x 6x 14x 14x 14x 14x 14x 4x 4x 4x 19x 19x 19x 4x 4x 4x 4x 7x 7x 7x 7x 7x 7x 7x 1x 1x 7x 1x 31x 31x 31x 31x 31x 31x 31x 6x 31x 11x 31x 5x 5x 4x 4x 4x 19x 19x 5x 1x 1x 1x 1x 1x 4x 4x 31x 8x 1x 1x 1x 1x 1x 7x 7x 31x 1x 1x 1x 31x 31x 1x 17x 17x 17x 17x 17x 17x 17x 17x 17x 87x 87x 1x 1x 1x 1x 87x 17x 17x 17x 17x 2x 2x 2x 8x 32x 32x 8x 2x 16x 17x 4x 4x 25x 25x 25x 25x 25x 348x 348x 348x 348x 348x 348x 348x 25x 25x 25x 25x 4x 4x 16x 16x 16x 17x 17x 17x 17x 85x 85x 17x 300x 300x 1452x 1452x 1452x 300x 300x 984x 984x 4944x 4944x 984x 984x 300x 1452x 1452x 4944x 4944x 1452x 1452x 300x 16x 16x 16x 16x 16x 16x 16x | /**
* Hierarchical forecast reconciliation strategies.
*
* Implements the four standard approaches from Hyndman et al. (2011):
*
* BU — Bottom-Up: P = [0 | I_n]
* OLS — Ordinary LS: P = (SᵀS)⁻¹ Sᵀ
* WLS — Weighted LS: P = (SᵀW⁻¹S)⁻¹ SᵀW⁻¹ (W diagonal)
* MinT — Min-Trace: P = (SᵀW⁻¹S)⁻¹ SᵀW⁻¹ (W full covariance)
*
* All strategies produce reconciled forecasts via:
* ŷ_h = S · P · ŷ_b
* where ŷ_b is the stacked base forecast vector (length m).
*/
import { Matrix, solve } from 'ml-matrix';
import type { ReconciliationStrategy, ReconcileOptions, BaseForecasts } from './types';
import type { SummingMatrixResult } from './summing-matrix';
// ---------------------------------------------------------------------------
// Helpers
// ---------------------------------------------------------------------------
/** Convert a Matrix to number[][] for serialisation / diagnostics. */
function matrixTo2D(m: Matrix): number[][] {
const rows: number[][] = [];
for (let r = 0; r < m.rows; r++) {
const row: number[] = [];
for (let c = 0; c < m.columns; c++) row.push(m.get(r, c));
rows.push(row);
}
return rows;
}
/** Build Matrix from number[][]. */
function matrixFrom2D(data: readonly (readonly number[])[]): Matrix {
const rows = data.length;
const cols = rows > 0 ? data[0].length : 0;
const m = new Matrix(rows, cols);
for (let r = 0; r < rows; r++) {
for (let c = 0; c < cols; c++) m.set(r, c, data[r][c]);
}
return m;
}
/** Build an m × m diagonal Matrix from a Float32Array of diagonal values. */
function diagFromArray(diag: Float32Array): Matrix {
const m = new Matrix(diag.length, diag.length);
for (let i = 0; i < diag.length; i++) m.set(i, i, diag[i]);
return m;
}
/** Ridge-regularised solve with automatic retry on failure. */
function solveWithRidge(A: Matrix, b: Matrix, ridge: number, maxRidge = 100): Matrix {
const p = A.rows;
const Aug = A.clone();
// Add ridge to diagonal
for (let i = 0; i < p; i++) {
Aug.set(i, i, Aug.get(i, i) + ridge);
}
try {
return solve(Aug, b);
} catch {
if (ridge >= maxRidge) {
throw new Error(
`Reconciliation failed: matrix remains singular after increasing ridge to ${ridge}. ` +
`Check for degenerate hierarchy or constant base forecasts.`,
);
}
return solveWithRidge(A, b, ridge * 10, maxRidge);
}
}
// ---------------------------------------------------------------------------
// Strategy implementations
// ---------------------------------------------------------------------------
/**
* Bottom-Up (BU): P = [0_{n×(m-n)} | I_n]
*
* Simply picks the bottom-level base forecasts — no reconciliation math
* beyond building the scaling matrix. BottomUp is the fastest strategy
* and guarantees coherence by construction.
*/
function projectionBU(S: Matrix): Matrix {
const m = S.rows;
const n = S.columns;
// P = I_n padded with zeros on the left → shape is [n × m], last n cols = I_n
const P = new Matrix(n, m);
for (let j = 0; j < n; j++) {
P.set(j, m - n + j, 1);
}
return P;
}
/**
* Ordinary Least Squares: P = (SᵀS)⁻¹ Sᵀ
*/
function projectionOLS(S: Matrix, ridge: number): Matrix {
const StS = S.transpose().mmul(S); // n × n
const St = S.transpose(); // n × m
return solveWithRidge(StS, St, ridge);
}
/**
* Weighted Least Squares: P = (SᵀW⁻¹S)⁻¹ SᵀW⁻¹
*
* W is diagonal with per-node residual variances.
*/
function projectionWLS(S: Matrix, W: Matrix, ridge: number): Matrix {
// W⁻¹ = diag(1/w_i)
const Winv = new Matrix(W.rows, W.columns);
for (let i = 0; i < W.rows; i++) {
const w = W.get(i, i);
Winv.set(i, i, w > 1e-12 ? 1 / w : 1);
}
// StWinvS = Sᵀ W⁻¹ S (n × n)
const StWinv = S.transpose().mmul(Winv);
const StWinvS = StWinv.mmul(S);
return solveWithRidge(StWinvS, StWinv, ridge);
}
/**
* Minimum Trace: P = (SᵀW⁻¹S)⁻¹ SᵀW⁻¹
*
* Same formula as WLS but W is the full residual covariance (m × m).
* Falls back to OLS when W is singular or not provided.
*/
function projectionMinT(S: Matrix, W: Matrix, ridge: number): Matrix {
// Try to invert W; fall back to OLS on failure
try {
const Winv = solve(W, Matrix.eye(W.rows));
const StWinv = S.transpose().mmul(Winv);
const StWinvS = StWinv.mmul(S);
return solveWithRidge(StWinvS, StWinv, ridge);
} catch {
return projectionOLS(S, ridge);
}
}
// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------
/**
* Compute the projection matrix P (n × m) for the given strategy.
*
* @param S The m × n summing matrix.
* @param strategy The reconciliation strategy.
* @param options Strategy-specific parameters:
* - `residualCovariance` — m × m covariance Matrix for MinT / WLS.
* - `ridge` — ridge regularisation added to diagonal before inversion (default 1e-6).
*
* @returns The n × m projection matrix P.
*/
export function computeProjectionMatrix(
S: Matrix,
strategy: ReconciliationStrategy,
options?: {
readonly residualCovariance?: Matrix;
readonly ridge?: number;
},
): Matrix {
const ridge = options?.ridge ?? 1e-6;
switch (strategy) {
case 'bu':
return projectionBU(S);
case 'ols':
return projectionOLS(S, ridge);
case 'wls': {
let W: Matrix;
if (options?.residualCovariance) {
// Use diagonal of provided covariance
const m = options.residualCovariance.rows;
W = new Matrix(m, m);
for (let i = 0; i < m; i++) {
W.set(i, i, Math.max(options.residualCovariance.get(i, i), 1e-12));
}
} else {
// No covariance → fall back to OLS with a warning
console.warn(
'[timesfm-hierarchical] WLS requires residualCovariance — falling back to OLS.',
);
return projectionOLS(S, ridge);
}
return projectionWLS(S, W, ridge);
}
case 'mint': {
if (!options?.residualCovariance) {
console.warn(
'[timesfm-hierarchical] MinT requires residualCovariance — falling back to OLS.',
);
return projectionOLS(S, ridge);
}
return projectionMinT(S, options.residualCovariance, ridge);
}
default: {
const exhaustive: never = strategy;
throw new Error(`Unknown reconciliation strategy: ${exhaustive}`);
}
}
}
// ---------------------------------------------------------------------------
// Reconcile base forecasts (model-free)
// ---------------------------------------------------------------------------
/**
* Reconcile base forecasts at all levels using the given strategy.
*
* This is a pure function that operates on pre-computed base forecasts.
* It does not require a TimesFM model — use {@link reconcileBaseForecasts}
* when you already have forecasts at all nodes (e.g., from a separate
* forecasting step).
*
* @param summing The summing matrix and node ordering (from buildSummingMatrix).
* @param baseForecasts Point forecasts per node id — every node must have an entry.
* @param options Reconciliation strategy and parameters.
*
* @returns Reconciled forecasts per node + projection matrix for diagnostics.
*
* @throws {Error} if any node id in summing is missing from baseForecasts.
*/
export function reconcileBaseForecasts(
summing: SummingMatrixResult,
baseForecasts: BaseForecasts,
options?: ReconcileOptions,
): ReconciledResult {
const { S: sMat, allNodeIds, bottomNodeIds } = summing;
const m = allNodeIds.length;
const n = bottomNodeIds.length;
// Validate every node has a base forecast
const horizon = baseForecasts[allNodeIds[0]]?.length;
for (const id of allNodeIds) {
const bf = baseForecasts[id];
if (!bf || bf.length !== horizon) {
throw new Error(
`Node "${id}" is missing from baseForecasts or has mismatched horizon length.`,
);
}
}
const strat: ReconciliationStrategy = options?.strategy ?? 'mint';
// Build S Matrix
const S = matrixFrom2D(sMat);
// Build residual covariance if provided
let residualCov: Matrix | undefined;
if (options?.residualCovariance && options.residualCovariance.length > 0) {
const cov = options.residualCovariance;
residualCov = new Matrix(cov.length, cov.length > 0 ? cov[0].length : 0);
for (let r = 0; r < cov.length; r++) {
for (let c = 0; c < cov[r].length; c++) {
residualCov.set(r, c, cov[r][c]);
}
}
}
// Auto-estimate WLS diagonal from base forecast variances when no covariance provided
let effectiveCov: Matrix | undefined = residualCov;
if ((strat === 'wls' || strat === 'mint') && !residualCov) {
// Estimate per-node variance from base forecasts for WLS
const diagVals = new Float32Array(m);
for (let i = 0; i < m; i++) {
const bf = baseForecasts[allNodeIds[i]];
let sum = 0,
sumSq = 0,
count = 0;
for (let h = 0; h < bf.length; h++) {
const v = bf[h];
if (Number.isFinite(v)) {
sum += v;
sumSq += v * v;
count++;
}
}
const mu = count > 0 ? sum / count : 0;
const variance = count > 1 ? Math.max((sumSq - count * mu * mu) / (count - 1), 1e-12) : 1e-12;
diagVals[i] = variance;
}
effectiveCov = diagFromArray(diagVals);
}
// Fallback for MinT without covariance: use the auto-estimated diagonal
const pOptions = {
residualCovariance: effectiveCov,
ridge: options?.ridge,
};
const P = computeProjectionMatrix(S, strat, pOptions);
// ── Reconcile each horizon step ─────────────────────────────────────────
const reconciled: Record<string, Float32Array> = {};
for (const id of allNodeIds) {
reconciled[id] = new Float32Array(horizon);
}
for (let h = 0; h < horizon; h++) {
// Stack base forecasts into ŷ_b: length-m vector
const yb: number[] = [];
for (let i = 0; i < m; i++) {
const bf = baseForecasts[allNodeIds[i]];
yb.push(bf[h]);
}
// ŷ_h = S · P · ŷ_b
// First: P · ŷ_b → bottom-level (length n)
const bottom: number[] = [];
for (let j = 0; j < n; j++) {
let val = 0;
for (let k = 0; k < m; k++) {
val += P.get(j, k) * yb[k];
}
bottom.push(val);
}
// Then: S · bottom → all-level (length m)
for (let i = 0; i < m; i++) {
let val = 0;
for (let j = 0; j < n; j++) {
val += sMat[i][j] * bottom[j];
}
reconciled[allNodeIds[i]][h] = val;
}
}
return {
reconciled,
projectionMatrix: matrixTo2D(P),
strategy: strat,
summingMatrix: sMat.map((row) => [...row]),
};
}
/** Return type of reconcileBaseForecasts (avoids inline type in JSDoc). */
interface ReconciledResult {
reconciled: Record<string, Float32Array>;
projectionMatrix: number[][];
strategy: ReconciliationStrategy;
summingMatrix: number[][];
}
|