"""Fit per-asset Tukey g-and-h return shapes to long US histories. The Monte Carlo engine draws each asset's annual log return as r = mu + sigma * (T(z) - m) / s, T(z) = (exp(g z) - 1) / g * exp(h z^2 / 2) where z is the asset's (correlated) standard normal draw, T is Tukey's g-and-h transform (Tukey 1977; Hoaglin 1985), and m, s are T's exact mean and standard deviation, so mu stays the expected log return and sigma the total volatility (``simulation_engine/tail_shapes.py``). g sets skew and h tail weight; g = h = 0 is the plain log-normal. This script estimates (g, h) for every asset class with at least ``MIN_YEARS`` of annual data by maximum likelihood, plus the "severe" alternative: the largest h inside the 95% profile-likelihood interval, i.e. the heaviest tail the record cannot statistically rule out. Asset classes with shorter histories borrow the shape of their closest long-history relative, or use the plain log-normal when there is none. It writes ``research/data/tail-shapes.json``, which the engine's constants are tested against, and the reader-facing calibration JSON and chart. Usage: python -m research.tail_calibration --out-dir frontend/public/blog-assets python tail_calibration.py --returns-dir . --out-dir . # published copy """ from __future__ import annotations import argparse import json from html import escape from math import comb, erf, sqrt from pathlib import Path from typing import Any, Callable import numpy as np DATA = Path(__file__).resolve().parent.parent / "simulation_engine/historical/data" SHAPES_JSON = Path(__file__).resolve().parent / "data/tail-shapes.json" TURBULENCE_JSON = Path(__file__).resolve().parent / "data/market-turbulence.json" TURBULENCE_WINDOW_YEARS = 30 LAST_YEAR = 2024 MIN_YEARS = 90 # long enough to contain 1929-32 and to say anything about tails CHI2_1_95 = 3.841 # Asset classes fitted directly, with the first year of their real data. FITTED = { "equity": 1871, "us_equity_small_cap": 1926, # 1871-1925 bond returns are a duration approximation from 10-year yields # with a quarter of the later volatility (1.7% vs 7.4%); pooling the two # regimes manufactures fat tails, so bonds start with measured returns. "bonds": 1926, "long_treasuries": 1928, } # h < 1/4 keeps the fourth moment finite (Martinez & Iglewicz 1984); a # simulator should not draw from a distribution with infinite kurtosis. H_MAX = 0.24 # Where a direct fit is not usable, a fitted class falls back to its parent. PARENT = {"us_equity_small_cap": "equity", "long_treasuries": "bonds"} # Shorter histories borrow a long-history relative's shape. BORROWED = { "intl_equity_developed": "equity", "emerging_markets": "equity", "reits": "equity", "tips": "bonds", "intl_bonds": "bonds", } # No long-history relative (or, for cash, returns too autocorrelated for an # i.i.d. shape to mean anything): plain log-normal. LOGNORMAL = ("cash", "gold", "commodities", "managed_futures") THRESHOLDS = (-0.10, -0.20, -0.30, -0.40) def annual_log_returns(name: str, first: int, last: int = LAST_YEAR, directory: Path = DATA) -> np.ndarray: rows = (line for line in (directory / f"{name}.csv").read_text().splitlines() if line[:1].isdigit()) series = {int(y): float(r) for y, r in (row.split(",")[:2] for row in rows)} return np.asarray(np.log1p([series[year] for year in range(first, last + 1)]), dtype=float) # --- g-and-h ------------------------------------------------------------------- def transform(z: np.ndarray, g: float, h: float) -> np.ndarray: core = z if abs(g) < 1e-9 else np.expm1(g * z) / g return np.asarray(core * np.exp(h * z * z / 2.0), dtype=float) def derivative(z: np.ndarray, g: float, h: float) -> np.ndarray: tail = np.exp(h * z * z / 2.0) if abs(g) < 1e-9: return np.asarray(tail * (1.0 + h * z * z), dtype=float) return np.asarray(tail * (np.exp(g * z) + h * z * np.expm1(g * z) / g), dtype=float) def invert(y: np.ndarray, g: float, h: float) -> np.ndarray: """Solve transform(z) = y by vectorized bisection (monotone for h >= 0).""" lo, hi = np.full_like(y, -12.0), np.full_like(y, 12.0) for _ in range(80): mid = (lo + hi) / 2.0 above = transform(mid, g, h) > y hi, lo = np.where(above, mid, hi), np.where(above, lo, mid) return (lo + hi) / 2.0 def moments(g: float, h: float) -> tuple[float, float]: """Exact mean and standard deviation of T(z), z ~ N(0, 1), for h < 1/2.""" if abs(g) < 1e-9: return 0.0, float((1.0 - 2.0 * h) ** -0.75) mean = (np.exp(g * g / (2.0 * (1.0 - h))) - 1.0) / (g * np.sqrt(1.0 - h)) second = (np.exp(2.0 * g * g / (1.0 - 2.0 * h)) - 2.0 * np.exp(g * g / (2.0 * (1.0 - 2.0 * h))) + 1.0) / ( g * g * np.sqrt(1.0 - 2.0 * h) ) return float(mean), float(np.sqrt(second - mean * mean)) def implied_moments(g: float, h: float) -> tuple[float, float]: """Exact skewness and excess kurtosis of T(z), z ~ N(0, 1), for h < 1/4. From the raw moments E[T^n] = sum_i (-1)^i C(n, i) exp(((n - i) g)^2 / (2 (1 - n h))) / (g^n sqrt(1 - n h)) (Martinez & Iglewicz 1984), or their g = 0 limit. """ def raw(n: int) -> float: if abs(g) < 1e-9: return 0.0 if n % 2 else float(np.prod(np.arange(1, n, 2)) / (1.0 - n * h) ** ((n + 1) / 2.0)) total = sum((-1) ** i * comb(n, i) * np.exp(((n - i) * g) ** 2 / (2.0 * (1.0 - n * h))) for i in range(n + 1)) return float(total / (g**n * np.sqrt(1.0 - n * h))) m1, m2, m3, m4 = raw(1), raw(2), raw(3), raw(4) var = m2 - m1 * m1 third = m3 - 3.0 * m1 * m2 + 2.0 * m1**3 fourth = m4 - 4.0 * m1 * m3 + 6.0 * m1 * m1 * m2 - 3.0 * m1**4 return float(third / var**1.5), float(fourth / var**2 - 3.0) def sample_moments(x: np.ndarray) -> tuple[float, float]: z = (x - x.mean()) / x.std() return float(np.mean(z**3)), float(np.mean(z**4) - 3.0) def bootstrap_moment_interval(x: np.ndarray, draws: int = 4000, seed: int = 0) -> dict[str, list[float]]: """95% percentile-bootstrap intervals for the sample skewness and excess kurtosis.""" rng = np.random.default_rng(seed) resampled = np.array([sample_moments(x[rng.integers(0, len(x), len(x))]) for _ in range(draws)]) lo, hi = np.percentile(resampled, [2.5, 97.5], axis=0) return {"skewness": [round(float(lo[0]), 3), round(float(hi[0]), 3)], "excess_kurtosis": [round(float(lo[1]), 3), round(float(hi[1]), 3)]} def log_likelihood(x: np.ndarray, a: float, b: float, g: float, h: float) -> float: if b <= 0 or h < 0 or h > H_MAX: return -np.inf z = invert((x - a) / b, g, h) return float(np.sum(-0.5 * z * z - 0.5 * np.log(2.0 * np.pi) - np.log(b * derivative(z, g, h)))) def nelder_mead(f: Callable[[np.ndarray], float], start: np.ndarray, step: np.ndarray, iters: int = 3000) -> tuple[np.ndarray, float]: """Minimize f; small and dependency-free, which is all four parameters need.""" simplex = [start] + [start + np.eye(len(start))[i] * step[i] for i in range(len(start))] values = [f(p) for p in simplex] for _ in range(iters): order = np.argsort(values) simplex, values = [simplex[i] for i in order], [values[i] for i in order] if abs(values[-1] - values[0]) < 1e-10: break centroid = np.mean(simplex[:-1], axis=0) reflected = centroid + (centroid - simplex[-1]) fr = f(reflected) if fr < values[0]: expanded = centroid + 2.0 * (centroid - simplex[-1]) fe = f(expanded) simplex[-1], values[-1] = (expanded, fe) if fe < fr else (reflected, fr) elif fr < values[-2]: simplex[-1], values[-1] = reflected, fr else: contracted = centroid + 0.5 * (simplex[-1] - centroid) fc = f(contracted) if fc < values[-1]: simplex[-1], values[-1] = contracted, fc else: simplex = [simplex[0] + 0.5 * (p - simplex[0]) for p in simplex] values = [f(p) for p in simplex] best = int(np.argmin(values)) return simplex[best], values[best] def fit(x: np.ndarray, h_fixed: float | None = None) -> dict[str, float]: """MLE of (a, b, g, h); with h_fixed, the profile over the other three.""" def negative(p: np.ndarray) -> float: a, log_b, g = p[0], p[1], p[2] h = h_fixed if h_fixed is not None else p[3] value = log_likelihood(x, a, float(np.exp(log_b)), g, h) return -value if np.isfinite(value) else 1e12 best: tuple[np.ndarray, float] | None = None for g0 in (-0.4, -0.1, 0.1, 0.4): start = [np.median(x), np.log(x.std()), g0] + ([] if h_fixed is not None else [0.05]) step = [0.01, 0.1, 0.1] + ([] if h_fixed is not None else [0.05]) candidate = nelder_mead(negative, np.array(start), np.array(step)) if best is None or candidate[1] < best[1]: best = candidate assert best is not None p, value = best return { "a": float(p[0]), "b": float(np.exp(p[1])), "g": float(p[2]), "h": float(h_fixed if h_fixed is not None else min(max(p[3], 0.0), H_MAX)), "log_likelihood": -float(value), } def severe_shape(x: np.ndarray, best: dict[str, float]) -> dict[str, float]: """Largest h whose profile likelihood stays inside the 95% interval.""" floor = best["log_likelihood"] - CHI2_1_95 / 2.0 inside = best for h in np.arange(round(best["h"], 2) + 0.01, H_MAX + 1e-9, 0.01): candidate = fit(x, h_fixed=float(h)) if candidate["log_likelihood"] < floor: break inside = candidate return inside def normal_log_likelihood(x: np.ndarray) -> float: var = x.var() return float(-0.5 * len(x) * (np.log(2.0 * np.pi * var) + 1.0)) def hoaglin(x: np.ndarray) -> dict[str, float]: """Quantile ("letter value") estimates of g and h (Hoaglin 1985). g is the median over tail depths of log(upper spread / lower spread) / z_p; h is the slope of log(g * spread / (e^{g z} - e^{-g z})) on z^2 / 2. Robust to the handful of extreme years that dominate a likelihood. """ median = float(np.quantile(x, 0.5)) depths = (1 / 4, 1 / 8, 1 / 16, 1 / 32) zs, spreads, gs = [], [], [] for p in depths: z = float(np.sqrt(2.0) * _erfinv(1.0 - 2.0 * p)) upper, lower = float(np.quantile(x, 1 - p)) - median, median - float(np.quantile(x, p)) zs.append(z) spreads.append(upper + lower) gs.append(np.log(upper / lower) / z) g = float(np.median(gs)) z_arr, s_arr = np.array(zs), np.array(spreads) y = np.log(g * s_arr / (np.exp(g * z_arr) - np.exp(-g * z_arr))) if abs(g) > 1e-9 else np.log(s_arr / (2 * z_arr)) slope = float(np.polyfit(z_arr**2 / 2.0, y, 1)[0]) return {"g": g, "h": max(slope, 0.0)} def _erfinv(y: float) -> float: """Inverse error function by Newton's method (no SciPy dependency).""" x = 0.0 for _ in range(60): x -= (erf(x) - y) / (2.0 / np.sqrt(np.pi) * np.exp(-x * x)) return x def candidate_shape(x: np.ndarray, best: dict[str, float], normal: float) -> tuple[str, dict[str, float] | None]: """The estimate to consider, or None when no estimate is usable. Maximum likelihood when its estimate is interior and beats the log-normal at the 5% level. If it runs to the finite-kurtosis bound with a flat profile, the quantile estimate, provided the likelihood cannot tell it from the maximum (inside the 95% interval). """ if 2.0 * (best["log_likelihood"] - normal) < CHI2_1_95: return "not distinguishable from log-normal", None if best["h"] < H_MAX - 0.005: return "maximum likelihood", best robust = hoaglin(x) if robust["h"] <= H_MAX: at_robust = fit(x, h_fixed=round(robust["h"], 3)) if best["log_likelihood"] - at_robust["log_likelihood"] <= CHI2_1_95 / 2.0: return "quantile (Hoaglin), likelihood flat up to the kurtosis bound", {"g": robust["g"], "h": robust["h"]} return "quantile estimate beyond the finite-kurtosis bound", None def choose_shape(x: np.ndarray, best: dict[str, float], normal: float, interval: dict[str, list[float]]) -> tuple[str, dict[str, float] | None]: """The candidate, if the record supports it; None falls back to the parent class. A likelihood fit only sees the years that happened, so a shape can match them well and still imply tails nothing in the record resembles. The candidate must also reproduce the sample's own skewness and excess kurtosis: its exact values must fall inside their 95% bootstrap intervals. For bonds this rejects a quantile fit implying an excess kurtosis above 100; their pooled fat tails partly reflect volatility that changed between decades, which a single i.i.d. shape cannot represent. """ method, chosen = candidate_shape(x, best, normal) if chosen is None: return method, None skew, kurt = implied_moments(chosen["g"], chosen["h"]) lo_s, hi_s = interval["skewness"] lo_k, hi_k = interval["excess_kurtosis"] if not (lo_s <= skew <= hi_s and lo_k <= kurt <= hi_k): return f"{method}, rejected: implied skewness {skew:+.2f} and excess kurtosis {kurt:+.1f} " \ "outside the sample's 95% bootstrap intervals", None return method, chosen # --- engine-equivalent tail probabilities --------------------------------------- def _phi_cdf(z: float) -> float: return 0.5 * (1.0 + erf(z / sqrt(2.0))) def tail_probability(threshold: float, g: float, h: float, mu: float, sigma: float) -> float: """P(simple return < threshold) for the engine's draw at (mu, sigma, g, h).""" m, s = moments(g, h) y = m + s * (np.log1p(threshold) - mu) / sigma return _phi_cdf(float(invert(np.array([y]), g, h)[0])) def stretch_tail_probability(threshold: float, factor: float, mu: float, sigma: float) -> float: """The same probability under the stretch this replaces (start -1, given factor).""" start = -1.0 below = _phi_cdf(start) tail_first_moment = -float(np.exp(-0.5) / np.sqrt(2.0 * np.pi)) shift = float(np.exp(-0.5) / np.sqrt(2.0 * np.pi)) + start * below + factor * (tail_first_moment - start * below) y = (np.log1p(threshold) - (mu - sigma * shift)) / sigma return _phi_cdf(start + (y - start) / factor if y < start else y) # --- run ------------------------------------------------------------------------- def turbulence(directory: Path = DATA) -> dict[str, Any]: """Volatility multipliers anchored to the calmest and roughest 30-year eras. The Market turbulence setting scales every asset's volatility by the ratio of US stock volatility in the calmest (or roughest) 30-year window on record to the volatility of the whole 1871-2024 record. Expected log returns and return shapes are unchanged, so it changes how bumpy the ride is, not long-run growth. """ x = annual_log_returns("equity", 1871, directory=directory) years = list(range(1871, LAST_YEAR + 1)) full = float(x.std(ddof=1)) windows = [ (years[i], years[i + TURBULENCE_WINDOW_YEARS - 1], float(x[i:i + TURBULENCE_WINDOW_YEARS].std(ddof=1))) for i in range(len(x) - TURBULENCE_WINDOW_YEARS + 1) ] calm = min(windows, key=lambda w: w[2]) rough = max(windows, key=lambda w: w[2]) presets = { "calmer": {"era": [calm[0], calm[1]], "stock_volatility": round(calm[2], 4), "multiplier": round(calm[2] / full, 3)}, "historical": {"era": [1871, LAST_YEAR], "stock_volatility": round(full, 4), "multiplier": 1.0}, "rougher": {"era": [rough[0], rough[1]], "stock_volatility": round(rough[2], 4), "multiplier": round(rough[2] / full, 3)}, } return {"window_years": TURBULENCE_WINDOW_YEARS, "full_record_stock_volatility": round(full, 4), "presets": presets} def calibrate(directory: Path = DATA) -> dict[str, Any]: fitted: dict[str, Any] = {} for asset, first in FITTED.items(): x = annual_log_returns(asset, first, directory=directory) z = (x - x.mean()) / x.std(ddof=1) best = fit(x) normal = normal_log_likelihood(x) interval = bootstrap_moment_interval(x) _, candidate = candidate_shape(x, best, normal) method, chosen = choose_shape(x, best, normal, interval) severe = severe_shape(x, best) if method == "maximum likelihood" else chosen or best fitted[asset] = { "method": method, "chosen": None if chosen is None else {"g": round(chosen["g"], 3), "h": round(chosen["h"], 3)}, "maximum_likelihood": {"g": round(best["g"], 3), "h": round(best["h"], 3)}, "quantile_estimate": {k: round(v, 3) for k, v in hoaglin(x).items()}, "years": [first, LAST_YEAR], "n": int(len(x)), "log_return_mean": round(float(x.mean()), 4), "log_return_volatility": round(float(x.std(ddof=1)), 4), "skewness": round(float(np.mean(z**3)), 3), "excess_kurtosis": round(float(np.mean(z**4) - 3.0), 3), "moment_interval_95": interval, "candidate_implied_moments": None if candidate is None else dict(zip( ("skewness", "excess_kurtosis"), (round(v, 3) for v in implied_moments(candidate["g"], candidate["h"])))), "severe_candidate": {"g": round(severe["g"], 3), "h": round(severe["h"], 3)}, "log_likelihood": { "normal": round(normal, 2), "g_and_h": round(best["log_likelihood"], 2), "severe": round(severe["log_likelihood"], 2) if "log_likelihood" in severe else None, }, "likelihood_ratio_vs_normal": round(2.0 * (best["log_likelihood"] - normal), 2), "aic": {"normal": round(4.0 - 2.0 * normal, 2), "g_and_h": round(8.0 - 2.0 * best["log_likelihood"], 2)}, } shapes: dict[str, Any] = {} for asset in FITTED: own = fitted[asset]["chosen"] if own is None: continue # "Severe" is about crashes, so only classes with a negative skew (a # left tail) take the heaviest tail the record allows; the rest keep # their fitted shape. severe = fitted[asset]["severe_candidate"] if own["g"] < 0 else own shapes[asset] = {"source": f"fitted 1871-2024 ({fitted[asset]['method']})" if FITTED[asset] == 1871 else f"fitted {FITTED[asset]}-{LAST_YEAR} ({fitted[asset]['method']})", "historical": own, "severe": severe} lognormal = {"g": 0.0, "h": 0.0} for asset in FITTED: if asset not in shapes and asset not in PARENT: shapes[asset] = {"source": f"log-normal ({fitted[asset]['method']})", "historical": lognormal, "severe": lognormal} for asset, parent in PARENT.items(): if asset not in shapes: shapes[asset] = {"source": f"borrowed from {parent} ({fitted[asset]['method']})", "historical": shapes[parent]["historical"], "severe": shapes[parent]["severe"]} for asset, relative in BORROWED.items(): shapes[asset] = {"source": f"borrowed from {relative}", "historical": shapes[relative]["historical"], "severe": shapes[relative]["severe"]} for asset in LOGNORMAL: shapes[asset] = {"source": "log-normal (no long-history relative)", "historical": {"g": 0.0, "h": 0.0}, "severe": {"g": 0.0, "h": 0.0}} # Tail frequencies for US stocks at the engine's default equity parameters: # the 1871-2024 mean and volatility of log returns, rounded to 0.1 point, # as research/return_calibration.py computes them for the engine. equity = annual_log_returns("equity", 1871, directory=directory) simple = np.expm1(equity) mu, sigma = round(float(equity.mean()), 3), round(float(equity.std(ddof=1)), 3) hist, sev = shapes["equity"]["historical"], shapes["equity"]["severe"] comparison = { "history": {f"{t:+.0%}": round(float(np.mean(simple < t)), 4) for t in THRESHOLDS}, "history_years": {f"{t:+.0%}": int(np.sum(simple < t)) for t in THRESHOLDS}, "none": {f"{t:+.0%}": round(tail_probability(t, 0.0, 0.0, mu, sigma), 4) for t in THRESHOLDS}, "historical": {f"{t:+.0%}": round(tail_probability(t, hist["g"], hist["h"], mu, sigma), 4) for t in THRESHOLDS}, "severe": {f"{t:+.0%}": round(tail_probability(t, sev["g"], sev["h"], mu, sigma), 4) for t in THRESHOLDS}, "previous_stretch_2_0": {f"{t:+.0%}": round(stretch_tail_probability(t, 2.0, mu, sigma), 4) for t in THRESHOLDS}, } turb = turbulence(directory) for preset in turb["presets"].values(): preset["us_stock_tail_frequencies"] = { f"{t:+.0%}": round(tail_probability(t, hist["g"], hist["h"], mu, sigma * preset["multiplier"]), 4) for t in THRESHOLDS } return { "model": "r = mu + sigma * (T(z) - m) / s, T(z) = (exp(g z) - 1) / g * exp(h z^2 / 2)", "references": [ "Tukey, J. W. (1977). Modern techniques in data analysis. NSF-sponsored regional research conference, Southeastern Massachusetts University.", "Hoaglin, D. C. (1985). Summarizing shape numerically: the g-and-h distributions. In Exploring Data Tables, Trends, and Shapes, Wiley.", "Martinez, J. and Iglewicz, B. (1984). Some properties of the Tukey g and h family of distributions. Communications in Statistics - Theory and Methods 13(3).", "Dutta, K. and Perry, J. (2006). A tale of tails: an empirical analysis of loss distribution models for estimating operational risk capital. Federal Reserve Bank of Boston Working Paper 06-13.", "Cont, R. (2001). Empirical properties of asset returns: stylized facts and statistical issues. Quantitative Finance 1(2).", ], "min_years_to_fit": MIN_YEARS, "max_h": H_MAX, "engine_equity_parameters": {"mu_log": mu, "sigma_total": sigma}, "fits": fitted, "shapes": shapes, "us_stock_tail_frequencies": comparison, "market_turbulence": turb, } SERIES = (("history", "History, 1871–2024", "#B47B16"), ("none", "Bell curve", "#9DB2CF"), ("historical", "Fitted (new default)", "#0F7B6C"), ("previous_stretch_2_0", "Previous model", "#B0532F")) def render_chart(result: dict[str, Any]) -> str: width, height = 1240, 560 left, top, bottom = 90, 150, 470 data = result["us_stock_tail_frequencies"] peak = max(v for key, _, _ in SERIES for v in data[key].values()) y_max = 0.02 * (int(peak / 0.02) + 1) group_w, bar_w = (width - left - 60) / len(THRESHOLDS), 34 def py(p: float) -> float: return bottom - (bottom - top) * p / y_max parts = [ f'', "How often a year's stock return falls below a threshold: history versus the model", 'For losses worse than 10, 20, 30 and 40 percent, the share of years in the 1871 to 2024 ' "record compared with the probability under a bell curve, under the fitted g-and-h model, and under the " "previous model.", "", '', 'Bad years: what history shows and what the model draws', 'Share of years an S&P total return fell below each threshold, ' "against the model's probability at a " f"{result['engine_equity_parameters']['mu_log']:.1%} log mean and " f"{result['engine_equity_parameters']['sigma_total']:.1%} volatility", ] tick = 0.0 while tick <= y_max + 1e-9: parts.append(f'' f'{tick:.0%}') tick += 0.02 for i, threshold in enumerate(THRESHOLDS): key = f"{threshold:+.0%}" gx = left + i * group_w + (group_w - 4 * (bar_w + 8)) / 2 for j, (series, _, colour) in enumerate(SERIES): value = data[series][key] x = gx + j * (bar_w + 8) parts.append(f'{value:.1%}') centre = gx + 2 * (bar_w + 8) - 4 parts.append(f'' f"Worse than {key.replace('-', '−')}" f'' f"{data['history_years'][key]} of 154 years") legend_x = 36 for _, label, colour in SERIES: parts.append(f'' f'{escape(label)}') legend_x += 30 + 7 * len(label) parts.append("") return "\n".join(parts) + "\n" def main() -> None: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--out-dir", required=True, type=Path) parser.add_argument("--returns-dir", type=Path, default=DATA, help="directory of year,return CSVs (default: the simulator's own series)") args = parser.parse_args() result = calibrate(args.returns_dir) args.out_dir.mkdir(parents=True, exist_ok=True) (args.out_dir / "tail-calibration.json").write_text(json.dumps(result, indent=1) + "\n") (args.out_dir / "tail-frequencies.svg").write_text(render_chart(result)) if args.returns_dir == DATA: SHAPES_JSON.parent.mkdir(parents=True, exist_ok=True) SHAPES_JSON.write_text(json.dumps(result["shapes"], indent=1, sort_keys=True) + "\n") TURBULENCE_JSON.write_text(json.dumps(result["market_turbulence"], indent=1, sort_keys=True) + "\n") print(json.dumps(result["shapes"], indent=1)) print(json.dumps(result["us_stock_tail_frequencies"], indent=1)) if __name__ == "__main__": main()