Probit PD Model

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

/** Deterministic D21-F01 Probability of Default reference calculations. */

type RecordValue = Record<string, unknown>;

function numberValue(name: string, value: unknown): number {
  if (typeof value !== "number" || !Number.isFinite(value)) throw new Error(`${name} must be a finite number`);
  return value;
}

function positive(name: string, value: unknown): number {
  const parsed = numberValue(name, value);
  if (parsed <= 0) throw new Error(`${name} must be positive`);
  return parsed;
}

function probability(name: string, value: unknown, openInterval = false): number {
  const parsed = numberValue(name, value);
  const valid = openInterval ? parsed > 0 && parsed < 1 : parsed >= 0 && parsed <= 1;
  if (!valid) throw new Error(`${name} must be ${openInterval ? "strictly " : ""}between zero and one`);
  return parsed;
}

function vector(name: string, value: unknown, minimum = 1): number[] {
  if (!Array.isArray(value) || value.length < minimum) throw new Error(`${name} must contain at least ${minimum} finite numbers`);
  return value.map((item, index) => numberValue(`${name}[${index}]`, item));
}

function clean(value: number): number {
  const rounded = Math.round((value + Number.EPSILON) * 1e12) / 1e12;
  return Object.is(rounded, -0) ? 0 : rounded;
}

function sigmoid(score: number): number {
  if (score >= 0) { const tail = Math.exp(-score); return 1 / (1 + tail); }
  const head = Math.exp(score); return head / (1 + head);
}

// Abramowitz-Stegun 7.1.26; adequate for the package's 1e-7 cross-language tolerance.
function erf(value: number): number {
  const sign = value < 0 ? -1 : 1;
  const x = Math.abs(value);
  const t = 1 / (1 + 0.3275911 * x);
  const polynomial = ((((1.061405429 * t - 1.453152027) * t + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t;
  const y = 1 - polynomial * Math.exp(-x * x);
  return sign * y;
}

function normalCdf(value: number): number { return 0.5 * (1 + erf(value / Math.sqrt(2))); }

function linearScore(intercept: unknown, coefficients: unknown, features: unknown): [number, number[]] {
  const constant = numberValue("intercept", intercept);
  const beta = vector("coefficients", coefficients);
  const x = vector("features", features);
  if (beta.length !== x.length) throw new Error("coefficients and features must have equal length");
  const contributions = beta.map((item, index) => clean(item * x[index]));
  return [clean(constant + contributions.reduce((sum, item) => sum + item, 0)), contributions];
}

export function logisticPdModel(intercept: unknown, coefficients: unknown, features: unknown, alert_threshold: unknown): RecordValue {
  const [score, contributions] = linearScore(intercept, coefficients, features);
  const threshold = probability("alert_threshold", alert_threshold, true);
  const pd = sigmoid(score);
  return {
    linear_score: clean(score), feature_contributions: contributions,
    probability_of_default: clean(pd), survival_probability: clean(1 - pd),
    odds_of_default: clean(pd / (1 - pd)), alert_threshold: clean(threshold),
    state: pd >= threshold ? "at-or-above-alert" : "below-alert",
    reason: "supplied-logit-score-transformed",
  };
}

export function probitPdModel(intercept: unknown, coefficients: unknown, features: unknown, alert_threshold: unknown): RecordValue {
  const [score, contributions] = linearScore(intercept, coefficients, features);
  const threshold = probability("alert_threshold", alert_threshold, true);
  const pd = normalCdf(score);
  return {
    latent_score: clean(score), feature_contributions: contributions,
    probability_of_default: clean(pd), survival_probability: clean(1 - pd),
    alert_threshold: clean(threshold), state: pd >= threshold ? "at-or-above-alert" : "below-alert",
    reason: "supplied-probit-score-transformed",
  };
}

export function throughTheCyclePd(annual_obligors: unknown, annual_defaults: unknown, current_year_index: unknown, minimum_years: unknown): RecordValue {
  const obligors = vector("annual_obligors", annual_obligors, 2);
  const defaults = vector("annual_defaults", annual_defaults, 2);
  if (obligors.length !== defaults.length) throw new Error("annual_obligors and annual_defaults must have equal length");
  const minimum = positive("minimum_years", minimum_years);
  if (!Number.isInteger(minimum) || obligors.length < minimum) throw new Error("minimum_years must be an integer no greater than the history length");
  const index = numberValue("current_year_index", current_year_index);
  if (!Number.isInteger(index) || index < 0 || index >= obligors.length) throw new Error("current_year_index must select an observed year");
  const rates = obligors.map((population, year) => {
    const count = defaults[year];
    if (population <= 0) throw new Error(`annual_obligors[${year}] must be positive`);
    if (!Number.isInteger(count) || count < 0 || count > population) throw new Error(`annual_defaults[${year}] must be an integer from zero to obligors`);
    return count / population;
  });
  const ttc = rates.reduce((sum, item) => sum + item, 0) / rates.length;
  const pooled = defaults.reduce((sum, item) => sum + item, 0) / obligors.reduce((sum, item) => sum + item, 0);
  const current = rates[index];
  return {
    annual_default_rates: rates.map(clean), through_the_cycle_pd: clean(ttc), pooled_default_rate: clean(pooled),
    current_observed_default_rate: clean(current), cycle_gap: clean(current - ttc), observation_years: rates.length,
    state: current > ttc ? "current-above-long-run" : current < ttc ? "current-below-long-run" : "current-equals-long-run",
    reason: "simple-average-of-annual-one-year-rates",
  };
}

export function pointInTimePd(through_the_cycle_pd_value: unknown, borrower_log_odds_shift: unknown, macro_factor_z: unknown, macro_sensitivity: unknown, alert_threshold: unknown): RecordValue {
  const ttc = probability("through_the_cycle_pd", through_the_cycle_pd_value, true);
  const borrowerShift = numberValue("borrower_log_odds_shift", borrower_log_odds_shift);
  const macro = numberValue("macro_factor_z", macro_factor_z);
  const sensitivity = numberValue("macro_sensitivity", macro_sensitivity);
  if (sensitivity < 0) throw new Error("macro_sensitivity must be nonnegative under this package convention");
  const threshold = probability("alert_threshold", alert_threshold, true);
  const baselineLogOdds = Math.log(ttc / (1 - ttc));
  const macroShift = sensitivity * macro;
  const pd = sigmoid(baselineLogOdds + borrowerShift + macroShift);
  return {
    through_the_cycle_pd: clean(ttc), baseline_log_odds: clean(baselineLogOdds), borrower_log_odds_shift: clean(borrowerShift),
    macro_log_odds_shift: clean(macroShift), point_in_time_pd: clean(pd), cycle_uplift: clean(pd - ttc),
    alert_threshold: clean(threshold), state: pd >= threshold ? "at-or-above-alert" : "below-alert",
    reason: "declared-log-odds-overlay-not-universal-ifrs-or-regulatory-formula",
  };
}

function mertonTerms(asset: number, sigma: number, debt: number, rate: number, horizon: number): [number, number] {
  const scale = sigma * Math.sqrt(horizon);
  const d1 = (Math.log(asset / debt) + (rate + 0.5 * sigma * sigma) * horizon) / scale;
  return [d1, d1 - scale];
}

export function mertonDistanceToDefault(equity_value: unknown, equity_volatility: unknown, debt_face_value: unknown, risk_free_rate: unknown, asset_drift: unknown, horizon_years: unknown, tolerance: unknown, max_iterations: unknown): RecordValue {
  const equity = positive("equity_value", equity_value), sigmaEquity = positive("equity_volatility", equity_volatility);
  const debt = positive("debt_face_value", debt_face_value), rate = numberValue("risk_free_rate", risk_free_rate);
  const drift = numberValue("asset_drift", asset_drift), horizon = positive("horizon_years", horizon_years);
  const tol = positive("tolerance", tolerance), limit = positive("max_iterations", max_iterations);
  if (!Number.isInteger(limit) || limit > 10000) throw new Error("max_iterations must be an integer no greater than 10000");
  let asset = equity + debt * Math.exp(-rate * horizon), sigmaAsset = Math.min(3, Math.max(1e-6, sigmaEquity * equity / asset));
  let residual = Number.POSITIVE_INFINITY, converged = false, iterations = 0;
  for (iterations = 1; iterations <= limit; iterations += 1) {
    const [d1, d2] = mertonTerms(asset, sigmaAsset, debt, rate, horizon), n1 = normalCdf(d1), n2 = normalCdf(d2);
    if (n1 <= 1e-14) break;
    const nextAsset = (equity + debt * Math.exp(-rate * horizon) * n2) / n1;
    const nextSigma = sigmaEquity * equity / (nextAsset * n1);
    residual = Math.max(Math.abs(nextAsset - asset) / Math.max(asset, 1), Math.abs(nextSigma - sigmaAsset));
    asset = nextAsset; sigmaAsset = nextSigma;
    if (residual <= tol) { converged = true; break; }
  }
  if (!converged) return { asset_value: null, asset_volatility: null, distance_to_default: null, physical_default_probability: null, risk_neutral_default_probability: null, iterations: Math.min(iterations, limit), residual: Number.isFinite(residual) ? clean(residual) : null, state: "not-converged", reason: "merton-equity-system-did-not-converge" };
  const [, d2] = mertonTerms(asset, sigmaAsset, debt, rate, horizon);
  const dd = (Math.log(asset / debt) + (drift - 0.5 * sigmaAsset * sigmaAsset) * horizon) / (sigmaAsset * Math.sqrt(horizon));
  return { asset_value: clean(asset), asset_volatility: clean(sigmaAsset), distance_to_default: clean(dd), physical_default_probability: clean(normalCdf(-dd)), risk_neutral_default_probability: clean(normalCdf(-d2)), iterations, residual: clean(residual), state: "converged", reason: "merton-equity-system-solved" };
}

export function chsDistressProbability(nimtaavg: unknown, tlmta: unknown, exretavg: unknown, sigma: unknown, rsize: unknown, cashmta: unknown, market_to_book: unknown, log_price: unknown): RecordValue {
  const values: Record<string, number> = { nimtaavg: numberValue("nimtaavg", nimtaavg), tlmta: numberValue("tlmta", tlmta), exretavg: numberValue("exretavg", exretavg), sigma: numberValue("sigma", sigma), rsize: numberValue("rsize", rsize), cashmta: numberValue("cashmta", cashmta), market_to_book: numberValue("market_to_book", market_to_book), log_price: numberValue("log_price", log_price) };
  if (values.sigma < 0 || values.tlmta < 0 || values.cashmta < 0) throw new Error("sigma, tlmta, and cashmta must be nonnegative");
  const beta: Record<string, number> = { nimtaavg: -20.264, tlmta: 1.416, exretavg: -7.129, sigma: 1.411, rsize: -0.045, cashmta: -2.132, market_to_book: 0.075, log_price: -0.058 };
  const contributions: Record<string, number> = {};
  Object.keys(values).forEach(key => { contributions[key] = clean(values[key] * beta[key]); });
  const score = -9.164 + Object.values(contributions).reduce((sum, item) => sum + item, 0);
  return { published_intercept: -9.164, score_contributions: contributions, failure_log_odds: clean(score), distress_probability: clean(sigmoid(score)), state: "published-score-replication", reason: "chs-table-4-twelve-month-lag-coefficients" };
}

export function bharathShumwayNaiveDistanceToDefault(equity_value: unknown, debt_face_value: unknown, equity_volatility: unknown, prior_year_equity_return: unknown, horizon_years: unknown): RecordValue {
  const equity = positive("equity_value", equity_value), debt = positive("debt_face_value", debt_face_value);
  const sigmaEquity = positive("equity_volatility", equity_volatility), drift = numberValue("prior_year_equity_return", prior_year_equity_return), horizon = positive("horizon_years", horizon_years);
  const firm = equity + debt, sigmaDebt = 0.05 + 0.25 * sigmaEquity;
  const sigmaAsset = (equity / firm) * sigmaEquity + (debt / firm) * sigmaDebt;
  const dd = (Math.log(firm / debt) + (drift - 0.5 * sigmaAsset * sigmaAsset) * horizon) / (sigmaAsset * Math.sqrt(horizon));
  return { naive_firm_value: clean(firm), naive_debt_volatility: clean(sigmaDebt), naive_asset_volatility: clean(sigmaAsset), naive_distance_to_default: clean(dd), naive_default_probability: clean(normalCdf(-dd)), state: "calculated", reason: "bharath-shumway-naive-approximation" };
}

export function calculate(topicId: string, inputs: RecordValue): RecordValue {
  if (!inputs || typeof inputs !== "object" || Array.isArray(inputs)) throw new Error("inputs must be an object");
  switch (topicId) {
    case "D21-F01-A01": return logisticPdModel(inputs.intercept, inputs.coefficients, inputs.features, inputs.alert_threshold);
    case "D21-F01-A02": return probitPdModel(inputs.intercept, inputs.coefficients, inputs.features, inputs.alert_threshold);
    case "D21-F01-A03": return throughTheCyclePd(inputs.annual_obligors, inputs.annual_defaults, inputs.current_year_index, inputs.minimum_years);
    case "D21-F01-A04": return pointInTimePd(inputs.through_the_cycle_pd_value, inputs.borrower_log_odds_shift, inputs.macro_factor_z, inputs.macro_sensitivity, inputs.alert_threshold);
    case "D21-F01-A05": return mertonDistanceToDefault(inputs.equity_value, inputs.equity_volatility, inputs.debt_face_value, inputs.risk_free_rate, inputs.asset_drift, inputs.horizon_years, inputs.tolerance, inputs.max_iterations);
    case "D21-F01-A06": return chsDistressProbability(inputs.nimtaavg, inputs.tlmta, inputs.exretavg, inputs.sigma, inputs.rsize, inputs.cashmta, inputs.market_to_book, inputs.log_price);
    case "D21-F01-A07": return bharathShumwayNaiveDistanceToDefault(inputs.equity_value, inputs.debt_face_value, inputs.equity_volatility, inputs.prior_year_equity_return, inputs.horizon_years);
    default: throw new Error(`unsupported topic_id: ${topicId}`);
  }
}