Skip to content

quanttoolbox.sustainable_finance

sustainable_finance.risk

Python alternatives

Keep — sector-decomposed quadratic-form risk building blocks (quadratic_form, quadratic_form_risk) and portfolio/sector modified-duration & DTS aggregation (bond_portfolio_metrics), generic to sector-based portfolio risk models. No general-purpose equivalent found; see Library alternatives for the full reasoning.

quanttoolbox.sustainable_finance.risk

Generic quadratic-form risk building blocks used across sector-based portfolio risk models -- e.g. bond-portfolio modified-duration/DTS risk in quanttoolbox.bond.pricing, which these functions were factored out of.

Ported from HSF toolbox hsf/{quadratic_form,quadratic_form_risk, bond_portfolio_metrics}.m.

Translation notes:

  • quadratic_form_risk.m builds a quadratic-form penalty for deviating a sector-level risk measure (e.g. modified duration) from per-sector targets: for each sector j it contributes an outer-product term Q_j = (s_j * risk)(s_j * risk)' (s_j the sector-j 0/1 indicator), summed across sectors into a single (Q, R, c) triple. Kept as a general-purpose helper rather than folded into bond/pricing.py, since the original file lives in hsf/ (sustainable-finance-general), not bond/, and nothing here is bond-specific.

BondPortfolioMetrics(md_portfolio, dts_portfolio, md_by_sector, dts_by_sector, unique_sector) dataclass

Portfolio- and sector-level modified duration (MD) and duration-times-spread (DTS).

QuadraticFormResult(qf, q, r, c, n, n_sector, unique_sector, q_j, r_j, c_j) dataclass

Sector-decomposed quadratic-risk form: qf = 0.5 w'Qw - w'R + c, with q_j/r_j/c_j holding each sector's own contribution to Q/R/c (summed into q/r/c).

bond_portfolio_metrics(sector, md, dts, w)

Portfolio- and sector-level modified duration and DTS, weighted by portfolio weights w.

Original: hsf/bond_portfolio_metrics.m

Source code in src/quanttoolbox/sustainable_finance/risk.py 
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
def bond_portfolio_metrics(
    sector: np.ndarray,
    md: np.ndarray,
    dts: np.ndarray,
    w: np.ndarray,
) -> BondPortfolioMetrics:
    """Portfolio- and sector-level modified duration and DTS, weighted by
    portfolio weights w.

    Original: hsf/bond_portfolio_metrics.m
    """
    sector = np.asarray(sector)
    md = np.asarray(md, dtype=float).flatten()
    dts = np.asarray(dts, dtype=float).flatten()
    w = np.asarray(w, dtype=float).flatten()

    unique_sector = np.unique(sector)
    n_sector = unique_sector.shape[0]

    md_portfolio = float(np.sum(w * md))
    dts_portfolio = float(np.sum(w * dts))

    md_by_sector = np.zeros(n_sector)
    dts_by_sector = np.zeros(n_sector)
    for j, s in enumerate(unique_sector):
        s_j = (sector == s).astype(float)
        md_by_sector[j] = np.sum(s_j * w * md)
        dts_by_sector[j] = np.sum(s_j * w * dts)

    return BondPortfolioMetrics(
        md_portfolio=md_portfolio,
        dts_portfolio=dts_portfolio,
        md_by_sector=md_by_sector,
        dts_by_sector=dts_by_sector,
        unique_sector=unique_sector,
    )

quadratic_form(x, q, r, c)

Evaluate the quadratic form qf(x) = 0.5 x'Qx - x'R + c.

Original: hsf/quadratic_form.m

Source code in src/quanttoolbox/sustainable_finance/risk.py 
27
28
29
30
31
32
33
34
35
def quadratic_form(x: np.ndarray, q: np.ndarray, r: np.ndarray, c: float) -> float:
    """Evaluate the quadratic form qf(x) = 0.5 x'Qx - x'R + c.

    Original: hsf/quadratic_form.m
    """
    x = np.asarray(x, dtype=float)
    q = np.asarray(q, dtype=float)
    r = np.asarray(r, dtype=float)
    return float(0.5 * x @ q @ x - x @ r + c)

quadratic_form_risk(sector, risk, risk_star, w=None)

Build a quadratic-form penalty for a sector-level risk measure (e.g. modified duration or DTS) deviating from per-sector targets risk_star, and (if w is given) evaluate it at portfolio weights w.

Original: hsf/quadratic_form_risk.m

Source code in src/quanttoolbox/sustainable_finance/risk.py 
 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
def quadratic_form_risk(
    sector: np.ndarray,
    risk: np.ndarray,
    risk_star: np.ndarray,
    w: np.ndarray | None = None,
) -> QuadraticFormResult:
    """Build a quadratic-form penalty for a sector-level risk measure (e.g.
    modified duration or DTS) deviating from per-sector targets
    `risk_star`, and (if `w` is given) evaluate it at portfolio weights w.

    Original: hsf/quadratic_form_risk.m
    """
    sector = np.asarray(sector)
    risk = np.asarray(risk, dtype=float).flatten()
    risk_star = np.asarray(risk_star, dtype=float).flatten()
    n = risk.shape[0]

    unique_sector = np.unique(sector)
    n_sector = unique_sector.shape[0]

    q = np.zeros((n, n))
    r = np.zeros(n)
    c = 0.0
    q_j = np.zeros((n, n, n_sector))
    r_j = np.zeros((n, n_sector))
    c_j = np.zeros(n_sector)

    for j, s in enumerate(unique_sector):
        s_j = (sector == s).astype(float)
        s_j_risk = s_j * risk

        q_j[:, :, j] = np.outer(s_j_risk, s_j_risk)
        r_j[:, j] = s_j_risk * risk_star[j]
        c_j[j] = 0.5 * risk_star[j] ** 2

        q += q_j[:, :, j]
        r += r_j[:, j]
        c += c_j[j]

    qf = quadratic_form(w, q, r, c) if w is not None else float("nan")

    return QuadraticFormResult(
        qf=qf,
        q=q,
        r=r,
        c=c,
        n=n,
        n_sector=n_sector,
        unique_sector=unique_sector,
        q_j=q_j,
        r_j=r_j,
        c_j=c_j,
    )

sustainable_finance.carbon

Python alternatives

Keep — closed-form cumulative-emissions ("carbon budget") integrals under linear, compound-decline, piecewise-linear, and GDP-adjusted trajectories. No general-purpose equivalent found; the piecewise closed-form sum and each other formula are cross-checked against scipy.integrate.quad in this module's tests. Two near-duplicate originals (carbon_budget_linear.m/carbon_budget_linear_trend.m) were merged, and one (carbon_budget_linear_reduction.m) was recognized as a strict special case of carbon_budget_Reduction.m and not ported separately — see the module docstring.

quanttoolbox.sustainable_finance.carbon

Cumulative carbon-budget calculations: the total emissions integral of CE(s) ds over [t0, t] under various assumed emissions trajectories CE(s) (constant decline rates, GDP-adjusted compound decline, and piecewise-linear historical/target paths).

Ported from HSF toolbox hsf/{carbon_budget_linear, carbon_budget_linear_trend,carbon_budget_linear_reduction, carbon_budget_piecewise,carbon_budget_compound_reduction, carbon_budget_Reduction}.m.

Translation notes:

  • carbon_budget_linear.m and carbon_budget_linear_trend.m are the exact same computation (beta0 * (t - t0) + 0.5 * beta1 * (t^2 - t0^2), the closed-form integral of beta0 + beta1 * s) -- the original files differ only in which of their two outputs (closed form vs. a scipy.integrate.quad-equivalent numerical cross-check) comes first. Merged into one function, carbon_budget_linear.
  • carbon_budget_linear_reduction.m is a strict special case of carbon_budget_Reduction.m's default ("linear rate") method, under the substitution r = reduction * ce_t0 (verified algebraically and in tests) -- not ported as a separate function; use carbon_budget(t0, t, ce_t0, r=reduction * ce_t0, method=1) instead.
  • Every original file returns both a closed-form value and a integral()-computed numerical cross-check of the same quantity (a self-verification pattern, not two different pieces of information). Only the closed-form value is exposed here; the numerical cross-check is instead used in this module's own test suite (via scipy.integrate.quad), which is a strictly better place for a redundant-by-construction correctness check than a return value every caller pays for.
  • MATLAB's integral maps to scipy.integrate.quad where a numerical integral genuinely differs from a closed form (carbon_budget_piecewise's piecewise-linear-interpolated integral, verified in tests against the closed-form sum).

carbon_budget(t0, t, ce_t0, r, method=1)

Cumulative emissions over [t0, t] starting from CE(t0), declining (or growing) at rate r under one of three conventions:

  • method=1 (default, "linear rate"): CE(s) = CE(t0) - r * (s - t0)
  • method=2 ("compound rate"): CE(s) = CE(t0) * (1 - r)^(s - t0)
  • method=3 ("growth rate"): CE(s) = CE(t0) * exp(-r * (s - t0))

Original: hsf/carbon_budget_Reduction.m

Source code in src/quanttoolbox/sustainable_finance/carbon.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def carbon_budget(t0: float, t: float, ce_t0: float, r: float, method: int = 1) -> float:
    """Cumulative emissions over [t0, t] starting from CE(t0), declining (or
    growing) at rate `r` under one of three conventions:

    - method=1 (default, "linear rate"): CE(s) = CE(t0) - r * (s - t0)
    - method=2 ("compound rate"): CE(s) = CE(t0) * (1 - r)^(s - t0)
    - method=3 ("growth rate"): CE(s) = CE(t0) * exp(-r * (s - t0))

    Original: hsf/carbon_budget_Reduction.m
    """
    dt = t - t0
    if method == 3:
        return (1.0 - np.exp(-r * dt)) / r * ce_t0
    if method == 2:
        return ((1.0 - r) ** dt - 1.0) / np.log(1.0 - r) * ce_t0
    return dt * ce_t0 - 0.5 * dt**2 * r

carbon_budget_compound_reduction(t0, t, delta_r, r_minus, ce_t0, g_y=None, discrete=False)

Cumulative emissions over [t0, t], for CE(s) declining at compound rate delta_r from an initial level CE(t0) * (1 - r_minus), optionally compounded against a GDP growth rate g_y (CE(s) then grows/shrinks at the net rate of decarbonization vs. growth).

discrete=False (default, matches the original) treats CE(s) as continuous and integrates the closed-form geometric series over [t0, t]. discrete=True instead sums CE(s) at each integer year s = t0, t0+1, ..., t (t0 and t must be integers) -- added for callers that budget year-by-year rather than continuously (e.g. HSF-Notebooks chapter 11g); the two conventions can differ by several percent on the same inputs, so this is a real modeling choice, not a numerical-precision detail.

Original: hsf/carbon_budget_compound_reduction.m (discrete=True is not part of the original, which is continuous-only)

Source code in src/quanttoolbox/sustainable_finance/carbon.py 
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
def carbon_budget_compound_reduction(
    t0: float,
    t: float,
    delta_r: float,
    r_minus: float,
    ce_t0: float,
    g_y: float | None = None,
    discrete: bool = False,
) -> float:
    """Cumulative emissions over [t0, t], for CE(s) declining at compound
    rate `delta_r` from an initial level ``CE(t0) * (1 - r_minus)``,
    optionally compounded against a GDP growth rate `g_y` (CE(s) then
    grows/shrinks at the net rate of decarbonization vs. growth).

    `discrete=False` (default, matches the original) treats `CE(s)` as
    continuous and integrates the closed-form geometric series over
    ``[t0, t]``. `discrete=True` instead sums `CE(s)` at each integer year
    ``s = t0, t0+1, ..., t`` (``t0`` and `t` must be integers) -- added for
    callers that budget year-by-year rather than continuously (e.g.
    HSF-Notebooks chapter 11g); the two conventions can differ by several
    percent on the same inputs, so this is a real modeling choice, not a
    numerical-precision detail.

    Original: hsf/carbon_budget_compound_reduction.m (`discrete=True` is
    not part of the original, which is continuous-only)
    """
    if discrete:
        years = np.arange(t0, t + 1)
        decay = (1.0 - delta_r) * (1.0 if g_y is None else 1.0 + g_y)
        return float(np.sum(ce_t0 * (1.0 - r_minus) * decay ** (years - t0)))

    dt = t - t0
    if g_y is None:
        return ((1.0 - delta_r) ** dt - 1.0) / np.log(1.0 - delta_r) * (1.0 - r_minus) * ce_t0
    return (
        ((1.0 + g_y) ** dt * (1.0 - delta_r) ** dt - 1.0)
        / (np.log(1.0 + g_y) + np.log(1.0 - delta_r))
        * (1.0 - r_minus)
        * ce_t0
    )

carbon_budget_linear(t0, t, beta0, beta1)

Cumulative emissions over [t0, t] under a linear emissions trend CE(s) = beta0 + beta1 * s.

Original: hsf/carbon_budget_linear.m (identical to hsf/carbon_budget_linear_trend.m -- see module docstring)

Source code in src/quanttoolbox/sustainable_finance/carbon.py 
43
44
45
46
47
48
49
50
def carbon_budget_linear(t0: float, t: float, beta0: float, beta1: float) -> float:
    """Cumulative emissions over [t0, t] under a linear emissions trend
    CE(s) = beta0 + beta1 * s.

    Original: hsf/carbon_budget_linear.m (identical to
    hsf/carbon_budget_linear_trend.m -- see module docstring)
    """
    return beta0 * (t - t0) + 0.5 * beta1 * (t**2 - t0**2)

carbon_budget_piecewise(t0, t, t_k, ce_k)

Cumulative emissions over [t0, t], for a piecewise-linear-interpolated emissions trajectory given at knots (t_k, ce_k).

Original: hsf/carbon_budget_piecewise.m (this is CB1/CB3, the two closed-form sums the original computes -- verified algebraically identical, see this module's tests)

Source code in src/quanttoolbox/sustainable_finance/carbon.py 
 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
def carbon_budget_piecewise(t0: float, t: float, t_k: np.ndarray, ce_k: np.ndarray) -> float:
    """Cumulative emissions over [t0, t], for a piecewise-linear-interpolated
    emissions trajectory given at knots (`t_k`, `ce_k`).

    Original: hsf/carbon_budget_piecewise.m (this is `CB1`/`CB3`, the two
    closed-form sums the original computes -- verified algebraically
    identical, see this module's tests)
    """
    t_k = np.asarray(t_k, dtype=float)
    ce_k = np.asarray(ce_k, dtype=float)

    ce_t0 = float(np.interp(t0, t_k, ce_k))
    ce_t = float(np.interp(t, t_k, ce_k))

    mask = (t_k >= t0) & (t_k <= t)
    t_k = t_k[mask]
    ce_k = ce_k[mask]

    if t_k.shape[0] == 0 or t_k[0] != t0:
        t_k = np.concatenate(([t0], t_k))
        ce_k = np.concatenate(([ce_t0], ce_k))
    if t_k[-1] != t:
        t_k = np.concatenate((t_k, [t]))
        ce_k = np.concatenate((ce_k, [ce_t]))

    t_k1, t_k2 = t_k[:-1], t_k[1:]
    ce_k1, ce_k2 = ce_k[:-1], ce_k[1:]
    dt_k = t_k2 - t_k1

    beta0 = (t_k2 / dt_k) * ce_k1 - (t_k1 / dt_k) * ce_k2
    beta1 = (ce_k2 - ce_k1) / dt_k

    segments = beta0 * (t_k2 - t_k1) + 0.5 * beta1 * (t_k2**2 - t_k1**2)
    return float(np.sum(segments))

sustainable_finance.esg

Python alternatives

Keep — Roncalli's implied-ESG-beta minimum-variance tilt (esg_beta_star, esg_minimum_variance) and Pedersen-Fitzgibbons-Pomorski (2021)'s ESG-efficient-frontier portfolio (pedersen_portfolio). No general-purpose equivalent found. hsf/cdp_filter.m (a CDP-dataset-specific data-loading/filtering script) was not ported — see the module docstring.

quanttoolbox.sustainable_finance.esg

ESG-tilted portfolio construction: the "implied ESG beta" a minimum-variance investor effectively targets (Roncalli's ESG beta-star model), the resulting minimum-variance-plus-ESG-tilt portfolio, Pedersen, Fitzgibbons & Pomorski (2021)'s ESG-efficient-frontier portfolio, ESG score-bucket transition/turnover matrices, and a Monte Carlo power-law-weighted-portfolio-score simulation.

Ported from HSF toolbox hsf/{compute_esg_beta_star, compute_esg_minimum_variance,compute_pedersen_portfolio}.m. esg_transition_matrix and mc_weighted_score have no standalone .m source in hfs-archive -- see their own docstrings -- and are instead promoted from HSF-Notebooks chapter 16a.

Translation notes:

  • hsf/cdp_filter.m (CDP -- Carbon Disclosure Project -- data loading and regional/sectoral filtering, plus trend estimation) is not ported: it loads a specific data/chap9_cdp3.mat dataset this package does not ship, and its region/sector categories are hardcoded to that one textbook chapter's data -- a data-loading script tied to one dataset, not a general-purpose library function (same category as the untranslated tools/ display helpers noted in docs/migration_map.md).
  • compute_esg_beta_star.m's local e = ones(n,1) is unrelated to compute_esg_minimum_variance.m's e parameter (a universe-selection mask) despite the shared MATLAB variable name -- kept as an internal local (ones_n) here to avoid the naming collision.
  • MATLAB's logical(e) (boolean-mask submatrix selection) is Sigma[mask][:, mask] in numpy, where mask = e.astype(bool).

EsgMinimumVarianceResult(x, sigma, sigma_x, beta_star, beta_esg_star, sigma_tilde_matrix, x_tilde) dataclass

Minimum-variance-plus-ESG-tilt portfolio weights x, the implied covariance matrix sigma, its portfolio volatility sigma_x, the beta-star pair the tilt was built from, and (sigma_tilde_matrix, x_tilde) the plain (untilted) minimum-variance portfolio restricted to the selected universe, for comparison.

EsgTransitionMatrixResult(p, cdf, to, to_system) dataclass

ESG score-bucket transition probabilities p (p[j, k] = the probability of moving from bucket j to bucket k on remeasurement), the equilibrium bucket probabilities cdf, the per-bucket turnover to = 1 - diag(p), and the aggregate (probability-weighted) turnover to_system.

PedersenPortfolioResult(w, w_r, sigma_bar, s_bar, lambda1, lambda2, pi_w, sigma_w, s_w, sr_w, c_x_y) dataclass

Pedersen-Fitzgibbons-Pomorski ESG-efficient portfolios: w (shape (n_assets, n_scenarios)) holds one portfolio per (sigma_bar, s_bar) scenario, with w_r its residual (cash/risk-free) weight and the rest the per-scenario risk/return/ESG-score/Sharpe-ratio diagnostics.

esg_beta_star(beta, sigma_m, beta_esg, sigma_esg, sigma_tilde)

The "implied" market beta and ESG beta (beta_star, beta_esg_star) at which a minimum-variance investor's residual risk is fully diversified, given per-asset market/ESG factor loadings (beta, beta_esg), factor volatilities (sigma_m, sigma_esg), and idiosyncratic volatilities (sigma_tilde).

Original: hsf/compute_esg_beta_star.m

Source code in src/quanttoolbox/sustainable_finance/esg.py 
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
def esg_beta_star(
    beta: np.ndarray,
    sigma_m: float,
    beta_esg: np.ndarray,
    sigma_esg: float,
    sigma_tilde: np.ndarray,
) -> tuple[float, float]:
    """The "implied" market beta and ESG beta (``beta_star``,
    ``beta_esg_star``) at which a minimum-variance investor's residual risk
    is fully diversified, given per-asset market/ESG factor loadings
    (`beta`, `beta_esg`), factor volatilities (`sigma_m`, `sigma_esg`), and
    idiosyncratic volatilities (`sigma_tilde`).

    Original: hsf/compute_esg_beta_star.m
    """
    beta = np.asarray(beta, dtype=float)
    beta_esg = np.asarray(beta_esg, dtype=float)
    sigma_tilde = np.asarray(sigma_tilde, dtype=float)
    ones_n = np.ones_like(beta)

    sigma_m_sqr = sigma_m**2
    sigma_esg_sqr = sigma_esg**2
    sigma_m_esg_sqr = (sigma_m * sigma_esg) ** 2

    beta_tilde = beta / sigma_tilde**2
    varphi_m = beta @ beta_tilde

    beta_esg_tilde = beta_esg / sigma_tilde**2
    varphi_esg = beta_esg @ beta_esg_tilde

    varphi_m_esg = beta @ beta_esg_tilde

    omega0 = (
        1.0
        + sigma_m_sqr * varphi_m
        + sigma_esg_sqr * varphi_esg
        + sigma_m_esg_sqr * (varphi_m * varphi_esg - varphi_m_esg**2)
    )
    omega1 = varphi_esg * (beta_tilde @ ones_n) - varphi_m_esg * (beta_esg_tilde @ ones_n)
    omega1 = sigma_m_sqr * ((beta_tilde @ ones_n) + sigma_esg_sqr * omega1)
    omega2 = varphi_m * (beta_esg_tilde @ ones_n) - varphi_m_esg * (beta_tilde @ ones_n)
    omega2 = sigma_esg_sqr * ((beta_esg_tilde @ ones_n) + sigma_m_sqr * omega2)

    beta_star = omega0 / omega1
    beta_esg_star = omega0 / omega2
    return float(beta_star), float(beta_esg_star)

esg_minimum_variance(beta, sigma_m, beta_esg, sigma_esg, sigma_tilde, e=None)

The minimum-variance portfolio tilted toward the ESG beta-star of esg_beta_star, restricted to an optional universe-selection mask e (default: the full universe).

Original: hsf/compute_esg_minimum_variance.m

Source code in src/quanttoolbox/sustainable_finance/esg.py 
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
def esg_minimum_variance(
    beta: np.ndarray,
    sigma_m: float,
    beta_esg: np.ndarray,
    sigma_esg: float,
    sigma_tilde: np.ndarray,
    e: np.ndarray | None = None,
) -> EsgMinimumVarianceResult:
    """The minimum-variance portfolio tilted toward the ESG beta-star of
    `esg_beta_star`, restricted to an optional universe-selection mask `e`
    (default: the full universe).

    Original: hsf/compute_esg_minimum_variance.m
    """
    beta = np.asarray(beta, dtype=float)
    beta_esg = np.asarray(beta_esg, dtype=float)
    sigma_tilde = np.asarray(sigma_tilde, dtype=float)
    n = beta.shape[0]

    d = np.diag(sigma_tilde**2)
    sigma_m_sqr = sigma_m**2
    sigma_esg_sqr = sigma_esg**2
    sigma = np.outer(beta, beta) * sigma_m_sqr + np.outer(beta_esg, beta_esg) * sigma_esg_sqr + d

    if e is None:
        e = np.ones(n)
    e = np.asarray(e, dtype=float)

    beta_masked = e * beta
    beta_esg_masked = e * beta_esg
    beta_star, beta_esg_star = esg_beta_star(
        beta_masked, sigma_m, beta_esg_masked, sigma_esg, sigma_tilde
    )

    x = 1.0 / sigma_tilde**2 * (1.0 - beta_masked / beta_star - beta_esg_masked / beta_esg_star)
    x = e * x
    x = x / np.sum(x)
    sigma_x = float(np.sqrt(x @ sigma @ x))

    mask = e.astype(bool)
    sigma_tilde_matrix = sigma[np.ix_(mask, mask)]
    n_tilde = sigma_tilde_matrix.shape[0]
    e_tilde = np.ones(n_tilde)
    inv_sigma_tilde = np.linalg.inv(sigma_tilde_matrix)
    x_tilde = (inv_sigma_tilde @ e_tilde) / (e_tilde @ inv_sigma_tilde @ e_tilde)

    return EsgMinimumVarianceResult(
        x=x,
        sigma=sigma,
        sigma_x=sigma_x,
        beta_star=beta_star,
        beta_esg_star=beta_esg_star,
        sigma_tilde_matrix=sigma_tilde_matrix,
        x_tilde=x_tilde,
    )

esg_transition_matrix(s, sigma)

The bucket-to-bucket transition probability matrix and turnover of an ESG score-bucket scheme with cut points s (K = 1 + len(s) buckets, delimited by -inf = s_0 < s_1 < ... < s_{K-1} < s_K = +inf), under a "true" ESG factor g1 ~ N(0, 1) and a noisy remeasurement g2 = g1 + sigma*eps (eps ~ N(0, 1) independent of g1), so (g1, g2) is bivariate normal with Cov(g1, g2) = 1 and Var(g2) = 1 + sigma**2.

p[j, k] = P(g2 in bucket k | g1 in bucket j) is computed via the two-sided bivariate-normal rectangle probability (_rectangle_prob). to = 1 - diag(p) is the per-bucket turnover (probability of not staying in the same bucket); to_system = sum(cdf * to) is the overall, probability-weighted turnover.

Not ported from a standalone MATLAB HSF toolbox library function -- no hsf/esg_transition_matrix.m file exists in hfs-archive. The identical algorithm does appear as a local (embedded) function inside HSF/16. Exercise Solution/chap16_chap2_exercise4.m, of which HSF-Notebooks chapter 16a is itself a direct Python port; promoted here from that notebook.

Source code in src/quanttoolbox/sustainable_finance/esg.py 
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
354
355
356
def esg_transition_matrix(s: np.ndarray, sigma: float) -> EsgTransitionMatrixResult:
    """The bucket-to-bucket transition probability matrix and turnover of
    an ESG score-bucket scheme with cut points `s` (`K = 1 + len(s)`
    buckets, delimited by `-inf = s_0 < s_1 < ... < s_{K-1} < s_K = +inf`),
    under a "true" ESG factor `g1 ~ N(0, 1)` and a noisy remeasurement
    `g2 = g1 + sigma*eps` (`eps ~ N(0, 1)` independent of `g1`), so
    `(g1, g2)` is bivariate normal with `Cov(g1, g2) = 1` and
    `Var(g2) = 1 + sigma**2`.

    `p[j, k] = P(g2 in bucket k | g1 in bucket j)` is computed via the
    two-sided bivariate-normal rectangle probability (`_rectangle_prob`).
    `to = 1 - diag(p)` is the per-bucket turnover (probability of *not*
    staying in the same bucket); `to_system = sum(cdf * to)` is the
    overall, probability-weighted turnover.

    Not ported from a standalone MATLAB HSF toolbox library function --
    no `hsf/esg_transition_matrix.m` file exists in `hfs-archive`. The
    identical algorithm does appear as a local (embedded) function inside
    `HSF/16. Exercise Solution/chap16_chap2_exercise4.m`, of which
    HSF-Notebooks chapter 16a is itself a direct Python port; promoted
    here from that notebook.
    """
    s = np.asarray(s, dtype=float)
    k = 1 + s.shape[0]
    s_ext = np.sort(np.concatenate([[-np.inf], s, [np.inf]]))

    p = np.zeros((k, k))
    cdf = np.zeros(k)
    mvn = multivariate_normal(
        mean=[0.0, 0.0], cov=[[1.0, 1.0], [1.0, 1.0 + sigma**2]], allow_singular=True
    )

    for j in range(k):
        cdf1 = norm.cdf(s_ext[j + 1]) - norm.cdf(s_ext[j])
        cdf[j] = cdf1
        for kk in range(k):
            cdf2 = _rectangle_prob(mvn, [s_ext[j], s_ext[kk]], [s_ext[j + 1], s_ext[kk + 1]])
            p[j, kk] = cdf2 / cdf1

    to = 1.0 - np.diag(p)
    to_system = float(np.sum(cdf * to))

    return EsgTransitionMatrixResult(p=p, cdf=cdf, to=to, to_system=to_system)

mc_weighted_score(alpha_vec, nS, n, rng, g1=None, g2=None)

Monte Carlo simulation of S(x) = sum(x*s), the score of a power-law-weighted portfolio, for nS draws at once (vectorized over the draw axis), for each power-law exponent in alpha_vec.

Portfolio weights x are built by drawing n iid uniforms u1, raising them to the power 1/alpha (a power-law reweighting of an otherwise-uniform allocation), and normalizing to sum to 1.

  • g1=None (default): u1 and the per-asset scores s are drawn fresh and independent of one another (s uncorrelated with the weights x).
  • g1 given: u1 = Phi(g1) (so the weights are tilted from the standard-normal factor g1) and s = g2 -- letting g1/g2 be correlated standard normals induces a correlation between the weights and the scores.

Not ported from a standalone MATLAB HSF toolbox library function -- no .m file in hfs-archive implements this as a reusable function. HSF/16. Exercise Solution/chap16_chap2_exercise4.m performs the same Monte Carlo computation inline, in per-scenario loops (Questions 2.e/2.g/2.h); HSF-Notebooks chapter 16a factored that inline computation into this function, and it is promoted here from that notebook.

Source code in src/quanttoolbox/sustainable_finance/esg.py 
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
def mc_weighted_score(
    alpha_vec: np.ndarray,
    nS: int,
    n: int,
    rng: np.random.Generator,
    g1: np.ndarray | None = None,
    g2: np.ndarray | None = None,
) -> np.ndarray:
    """Monte Carlo simulation of `S(x) = sum(x*s)`, the score of a
    power-law-weighted portfolio, for `nS` draws at once (vectorized over
    the draw axis), for each power-law exponent in `alpha_vec`.

    Portfolio weights `x` are built by drawing `n` iid uniforms `u1`,
    raising them to the power `1/alpha` (a power-law reweighting of an
    otherwise-uniform allocation), and normalizing to sum to 1.

    - `g1=None` (default): `u1` and the per-asset scores `s` are drawn
      fresh and independent of one another (`s` uncorrelated with the
      weights `x`).
    - `g1` given: `u1 = Phi(g1)` (so the weights are tilted from the
      standard-normal factor `g1`) and `s = g2` -- letting `g1`/`g2` be
      correlated standard normals induces a correlation between the
      weights and the scores.

    Not ported from a standalone MATLAB HSF toolbox library function --
    no `.m` file in `hfs-archive` implements this as a reusable function.
    `HSF/16. Exercise Solution/chap16_chap2_exercise4.m` performs the same
    Monte Carlo computation inline, in per-scenario loops (Questions
    2.e/2.g/2.h); HSF-Notebooks chapter 16a factored that inline
    computation into this function, and it is promoted here from that
    notebook.
    """
    if g1 is None:
        u1 = rng.random((n, nS))
        s = rng.standard_normal((n, nS))
    else:
        u1 = norm.cdf(g1)
        s = g2

    s_x = np.empty((nS, len(alpha_vec)))
    for j, a in enumerate(alpha_vec):
        v = u1 ** (1.0 / a)
        x = v / v.sum(axis=0, keepdims=True)
        s_x[:, j] = (x * s).sum(axis=0)
    return s_x

pedersen_portfolio(mu, r, sigma, s, sigma_bar, s_bar)

Pedersen, Fitzgibbons & Pomorski (2021)'s ESG-efficient-frontier portfolio: mean-variance-optimal subject to a target volatility sigma_bar and a target portfolio ESG score s_bar, swept over one portfolio per (sigma_bar, s_bar) pair (broadcast to a common length).

Original: hsf/compute_pedersen_portfolio.m

Source code in src/quanttoolbox/sustainable_finance/esg.py 
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
def pedersen_portfolio(
    mu: np.ndarray,
    r: float,
    sigma: np.ndarray,
    s: np.ndarray,
    sigma_bar: np.ndarray | float,
    s_bar: np.ndarray | float,
) -> PedersenPortfolioResult:
    """Pedersen, Fitzgibbons & Pomorski (2021)'s ESG-efficient-frontier
    portfolio: mean-variance-optimal subject to a target volatility
    `sigma_bar` and a target portfolio ESG score `s_bar`, swept over one
    portfolio per (`sigma_bar`, `s_bar`) pair (broadcast to a common
    length).

    Original: hsf/compute_pedersen_portfolio.m
    """
    mu = np.asarray(mu, dtype=float)
    s = np.asarray(s, dtype=float)
    sigma = np.asarray(sigma, dtype=float)
    n = mu.shape[0]

    pi_ = mu - r
    inv_sigma = np.linalg.inv(sigma)
    ones_n = np.ones(n)

    c_1_pi = ones_n @ inv_sigma @ pi_
    c_1_s = ones_n @ inv_sigma @ s
    c_s_pi = s @ inv_sigma @ pi_
    c_s_s = s @ inv_sigma @ s
    c_1_1 = ones_n @ inv_sigma @ ones_n
    c_pi_pi = pi_ @ inv_sigma @ pi_
    c_x_y = np.array([c_1_pi, c_s_pi, c_s_s, c_1_s, c_1_1, c_pi_pi])

    sigma_bar_arr = np.atleast_1d(np.asarray(sigma_bar, dtype=float))
    s_bar_arr = np.atleast_1d(np.asarray(s_bar, dtype=float))
    n_iters = max(sigma_bar_arr.shape[0], s_bar_arr.shape[0])
    sigma_bar_arr = np.broadcast_to(sigma_bar_arr, (n_iters,))
    s_bar_arr = np.broadcast_to(s_bar_arr, (n_iters,))

    all_lambda1 = np.zeros(n_iters)
    all_lambda2 = np.zeros(n_iters)
    all_w = np.zeros((n, n_iters))
    all_pi_w = np.zeros(n_iters)
    all_sigma_w = np.zeros(n_iters)
    all_s_w = np.zeros(n_iters)
    all_sr_w = np.zeros((n_iters, 2))

    for it in range(n_iters):
        sb, tb = sigma_bar_arr[it], s_bar_arr[it]
        denom = c_s_s - 2 * c_1_s * tb + c_1_1 * tb**2
        lambda2 = (c_1_pi * tb - c_s_pi) / denom
        aux = c_pi_pi - (c_1_pi * tb - c_s_pi) ** 2 / denom
        lambda1 = -1.0 / (2 * sb) * np.sqrt(aux)
        w = (-1.0 / (2 * lambda1)) * inv_sigma @ (pi_ + lambda2 * (s - tb))

        all_lambda1[it] = lambda1
        all_lambda2[it] = lambda2
        all_w[:, it] = w
        all_pi_w[it] = w @ pi_
        all_sigma_w[it] = np.sqrt(w @ sigma @ w)
        all_s_w[it] = (w @ s) / np.sum(w)
        all_sr_w[it, 0] = (w @ pi_) / sb
        all_sr_w[it, 1] = np.sqrt(aux)

    w_r = 1.0 - np.sum(all_w, axis=0)

    return PedersenPortfolioResult(
        w=all_w,
        w_r=w_r,
        sigma_bar=sigma_bar_arr,
        s_bar=s_bar_arr,
        lambda1=all_lambda1,
        lambda2=all_lambda2,
        pi_w=all_pi_w,
        sigma_w=all_sigma_w,
        s_w=all_s_w,
        sr_w=all_sr_w,
        c_x_y=c_x_y,
    )

sustainable_finance.climate

Python alternatives

Keep — the DICE (Dynamic Integrated Climate-Economy) model's carbon-cycle/temperature state-transition matrices and forward simulation. No general-purpose equivalent found; physical constants are Nordhaus's DICE-2016 calibration, hardcoded as in the original.

quanttoolbox.sustainable_finance.climate

The DICE (Dynamic Integrated Climate-Economy) model's carbon-cycle and temperature submodules: the state-transition matrices linking industrial emissions to atmospheric carbon concentration and global temperature, and a full forward simulation of both given exogenous GDP and mitigation-rate paths.

Ported from HSF toolbox hsf/{dice_temperature_matrix, dice_temperature_simulation}.m.

Translation notes:

  • All physical constants (carbon-cycle transfer coefficients, radiative forcing parameters, initial 2015 conditions) are Nordhaus's DICE-2016 calibration, hardcoded in the original -- kept as-is (not exposed as parameters), matching the original's scope.
  • dice_temperature_simulation.m accepts a parameters argument that is never read anywhere in the function body (a vestigial/dead parameter in the original) -- dropped here rather than carried forward as a do-nothing kwarg.
  • Y_fn/mu_fn (GDP and mitigation-rate paths, as functions of time) are plain Python callables, matching the original's function-handle arguments exactly.

DiceTemperatureMatrices(phi_cc, b_cc, xi_t, b_t, xi_t_5, b_t_5, xi1, xi2, xi3, xi4, c_at, c_lo, lambda_, beta) dataclass

The DICE carbon-cycle (phi_cc, b_cc) and temperature (xi_t, b_t) state-transition matrices for one time step of length delta_t, plus the intermediate physical constants they were built from.

dice_temperature_matrix(delta_t, scale=None, method=1)

Build the DICE model's carbon-cycle and temperature state-transition matrices for a time step of length delta_t.

scale converts delta_t into seconds (default: delta_t is in years, scale = 365.25 * 24 * 3600). method=2 instead derives xi_t from the standard 5-year-step calibration matrix via a matrix power (xi_t_5 ** (delta_t_years / 5)) rather than reconstructing it from the underlying continuous-time physical constants.

Original: hsf/dice_temperature_matrix.m

Source code in src/quanttoolbox/sustainable_finance/climate.py 
 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
def dice_temperature_matrix(
    delta_t: float, scale: float | None = None, method: int = 1
) -> DiceTemperatureMatrices:
    """Build the DICE model's carbon-cycle and temperature state-transition
    matrices for a time step of length `delta_t`.

    `scale` converts `delta_t` into seconds (default: `delta_t` is in
    years, `scale = 365.25 * 24 * 3600`). `method=2` instead derives `xi_t`
    from the standard 5-year-step calibration matrix via a matrix power
    (``xi_t_5 ** (delta_t_years / 5)``) rather than reconstructing it from
    the underlying continuous-time physical constants.

    Original: hsf/dice_temperature_matrix.m
    """
    if scale is None:
        scale = 365.25 * 24 * 3600

    phi1 = 0.2727
    phi_cc = np.array(
        [
            [0.9120, 0.0383, 0.0],
            [0.0880, 0.9592, 0.0003],
            [0.0, 0.0025, 0.9997],
        ]
    )
    b_cc = np.array([phi1, 0.0, 0.0])

    xi_t_5 = np.array([[86.30, 0.8624], [2.50, 97.50]]) / 100.0
    b_t_5 = np.array([0.098, 0.0])

    delta_5_seconds = 5.0 * scale
    delta_t_seconds = delta_t * scale

    xi1, xi2, xi3, xi4 = 0.098, 3.8 / 2.9, 0.088, 0.025
    c_at = delta_5_seconds / xi1
    lambda_ = xi2
    beta = xi3
    c_lo = delta_5_seconds * beta / xi4

    xi1_prime = 1.0 - (lambda_ + beta) * delta_t_seconds / c_at
    xi2_prime = beta * delta_t_seconds / c_at
    xi3_prime = beta * delta_t_seconds / c_lo
    xi4_prime = 1.0 - beta * delta_t_seconds / c_lo
    xi_t = np.array([[xi1_prime, xi2_prime], [xi3_prime, xi4_prime]])
    b_t = np.array([delta_t_seconds / c_at, 0.0])

    if method == 2:
        xi_t = fractional_matrix_power(xi_t_5, 1.0 / 5.0).real

    return DiceTemperatureMatrices(
        phi_cc=phi_cc,
        b_cc=b_cc,
        xi_t=xi_t,
        b_t=b_t,
        xi_t_5=xi_t_5,
        b_t_5=b_t_5,
        xi1=xi1,
        xi2=xi2,
        xi3=xi3,
        xi4=xi4,
        c_at=c_at,
        c_lo=c_lo,
        lambda_=lambda_,
        beta=beta,
    )

dice_temperature_simulation(t0, t_end, delta_t, y_fn, mu_fn, numeric=False)

Forward-simulate the DICE model's industrial emissions, carbon-cycle concentrations, radiative forcing, and temperature, from t0 to t_end in steps of delta_t, given a GDP path y_fn(t) and a mitigation-rate path mu_fn(t).

numeric=True uses the 5-year-calibration temperature matrices (xi_t_5, b_t_5) directly instead of the continuous-time reconstruction from dice_temperature_matrix.

Returns an array of shape (n_iters + 1, 10) with columns [t, CE_t, sigma_t, CC_AT_t, CC_UP_t, CC_LO_t, F_EX_t, F_RAD_t, T_AT_t, T_LO_t].

Original: hsf/dice_temperature_simulation.m (the unused parameters argument is dropped -- see module docstring)

Source code in src/quanttoolbox/sustainable_finance/climate.py 
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
def dice_temperature_simulation(
    t0: float,
    t_end: float,
    delta_t: float,
    y_fn: Callable[[float], float],
    mu_fn: Callable[[float], float],
    numeric: bool = False,
) -> np.ndarray:
    """Forward-simulate the DICE model's industrial emissions, carbon-cycle
    concentrations, radiative forcing, and temperature, from `t0` to
    `t_end` in steps of `delta_t`, given a GDP path `y_fn(t)` and a
    mitigation-rate path `mu_fn(t)`.

    `numeric=True` uses the 5-year-calibration temperature matrices
    (`xi_t_5`, `b_t_5`) directly instead of the continuous-time
    reconstruction from `dice_temperature_matrix`.

    Returns an array of shape ``(n_iters + 1, 10)`` with columns
    ``[t, CE_t, sigma_t, CC_AT_t, CC_UP_t, CC_LO_t, F_EX_t, F_RAD_t,
    T_AT_t, T_LO_t]``.

    Original: hsf/dice_temperature_simulation.m (the unused `parameters`
    argument is dropped -- see module docstring)
    """
    n_iters = int(round((t_end - t0) / delta_t))

    ce_land_0 = 3.3
    delta_land = 0.20
    sigma_0 = 0.5491
    g_sigma_0 = 0.01
    delta_sigma = 0.001

    cc_0 = np.array([830.4, 1527.0, 10010.0])

    eta = 3.8
    cc_at_1750 = 588.0
    f_ex_0 = 0.25
    f_ex_2100 = 0.70
    delta_f_ex = 5.0 * (f_ex_2100 - f_ex_0) / 90.0
    f_rad_0 = (eta / np.log(2)) * np.log(cc_0[0] / cc_at_1750) + f_ex_0

    t_0_vec = np.array([0.8, 0.0068])

    matrices = dice_temperature_matrix(delta_t, scale=1.0)
    phi_cc, b_cc = matrices.phi_cc, matrices.b_cc
    if numeric:
        xi_t = np.array([[86.30, 0.8624], [2.50, 97.50]]) / 100.0
        b_t = np.array([0.098, 0.0])
    else:
        xi_t, b_t = matrices.xi_t, matrices.b_t

    y_0 = y_fn(t0)
    mu_0 = mu_fn(t0)
    ce_industry_0 = (1 - mu_0) * sigma_0 * y_0
    ce_0 = ce_industry_0 + ce_land_0

    t = t0
    sigma_t = sigma_0
    ce_land_t = ce_land_0
    ce_t = ce_0
    g_sigma_t = g_sigma_0
    cc_t = cc_0.copy()
    f_ex_t = f_ex_0
    f_rad_t = f_rad_0
    t_t = t_0_vec.copy()

    results = np.zeros((n_iters + 1, 10))
    results[0, :] = np.concatenate(([t, ce_t, sigma_t], cc_t, [f_ex_t, f_rad_t], t_t))

    for it in range(1, n_iters + 1):
        t = t + delta_t

        y_t = y_fn(t)
        mu_t = mu_fn(t)

        ce_industry_t = (1 - mu_t) * sigma_t * y_t
        ce_land_t = ce_land_t * (1 - delta_land)
        ce_t = ce_industry_t + ce_land_t

        g_sigma_t = 1.0 / (1 + delta_sigma) * g_sigma_t
        sigma_t = (1 + g_sigma_t) * sigma_t

        cc_t = phi_cc @ cc_t + b_cc * ce_t
        cc_at_t = cc_t[0]

        if t <= 2100:
            f_ex_t = f_ex_t + delta_f_ex
        f_rad_t = (eta / np.log(2)) * np.log(cc_at_t / cc_at_1750) + f_ex_t

        t_t = xi_t @ t_t + b_t * f_rad_t
        results[it, :] = np.concatenate(([t, ce_t, sigma_t], cc_t, [f_ex_t, f_rad_t], t_t))

    return results

sustainable_finance.ecology

Python alternatives

Keep — species-area/endemics-area relationships, species-abundance-distribution histogram binning (including Preston's log2 "octave" classes), and Hurlbert's rarefaction estimator. No general-purpose equivalent found (these are ecology-specific biodiversity measures, not general statistics).

quanttoolbox.sustainable_finance.ecology

Species-area/endemics-area relationships, species-abundance distributions, Hurlbert's rarefaction estimator, and a colonization- extinction equilibrium model's net species-richness growth rate -- classic biodiversity measures from ecology, used in this toolbox's biodiversity-risk chapters.

Ported from HSF toolbox hsf/{species_area_relationship, endemics_area_relationship,species_abundance_distribution,hurlbert}.m. fn_delta_t has no standalone .m source in hfs-archive -- see its own docstring -- and is instead promoted from HSF-Notebooks chapter 16c.

Translation notes:

  • species_area_relationship/endemics_area_relationship accept a as a scalar or array of target sub-areas, and can work either directly from per-species individual counts n_i, or from a pre-binned species-abundance histogram (s_j = number of species with exactly j individuals) -- typically the output of species_abundance_distribution. Pass s_j=None for the first mode.
  • species_abundance_distribution.m's three-way branch on its second argument (brk) is preserved as three explicit modes rather than MATLAB's type-based dispatch: breaks=None (unique-count histogram), breaks="octave" (Preston's log2 "octave" abundance classes), or breaks=<array> (custom breakpoints). The "octave" and <array> modes' bucket-assignment rule was corrected to match numpy's own np.histogram half-open-bin convention exactly (floor(log2(n)) for octaves; np.histogram(n_i, bins=[0, *breaks]) for custom breakpoints) -- the original bucket-assignment rule here used a different, non-numpy-native inclusive convention, which HSF-Notebooks chapter 05c's independent reimplementation of this same function was found to disagree with on real (Whittaker/BCI) survey data before this fix; see species_abundance_distribution's own docstring.
  • hurlbert.m's method=2 (log-gamma) branch is numerically more stable than method=1's direct scipy.special.comb ratio for large sample sizes (avoids overflow in the individual binomial coefficients); both compute the same quantity and are verified to agree in this module's tests.

endemics_area_relationship(n_i, s_j, area_total, area)

Expected number of species confined entirely within a sub-area area out of a total surveyed area area_total ("endemics"), given either per-species individual counts n_i (s_j=None), or a species-abundance histogram s_j.

Original: hsf/endemics_area_relationship.m

Source code in src/quanttoolbox/sustainable_finance/ecology.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def endemics_area_relationship(
    n_i: np.ndarray, s_j: np.ndarray | None, area_total: float, area: np.ndarray | float
) -> np.ndarray:
    """Expected number of species confined entirely within a sub-area
    `area` out of a total surveyed area `area_total` ("endemics"), given
    either per-species individual counts `n_i` (`s_j=None`), or a
    species-abundance histogram `s_j`.

    Original: hsf/endemics_area_relationship.m
    """
    area = np.asarray(area, dtype=float)
    ratio = area / area_total

    if s_j is None:
        n_i = np.asarray(n_i, dtype=float)
        return np.sum(ratio[..., None] ** n_i, axis=-1)

    s_j = np.asarray(s_j, dtype=float)
    j = np.arange(1, s_j.shape[0] + 1)
    return np.sum(s_j * ratio[..., None] ** j, axis=-1)

fn_delta_t(S_t, lambda_0, beta1, mu_s, beta2, beta3, beta4, S_s)

Net species-richness growth rate delta(S) = lambda(S) - mu(S) of a colonization-extinction equilibrium model (the theory of island biogeography, TIB), where mu(S) = mu^long(S) + mu^short(S):

  • lambda(S): colonization rate, =lambda_0 at S=0, decaying to 0 at the saturation richness S=S_s.
  • mu^long(S): long-term extinction rate, =0 at S=0, rising to mu_s at S=S_s.
  • mu^short(S) = beta3*lambda(S)*exp(-beta4*S): a short-term extinction component proportional to the colonization rate (recently-arrived species are more extinction-prone).

S_t is clipped to S_t = max(0, S_t) before evaluation, so the ODE/root-finder driving dS/dt = delta(S) never evaluates the model at a negative richness -- this clip protects only this one function evaluation, not the ODE trajectory as a whole (an unclipped S(t) can still run negative between evaluations if delta(S) stays negative all the way down).

Returns (delta_t, lambda_t, mu_long_t, mu_short_t, mu_t). Used for equilibrium root-finding (delta_t == 0), ODE trajectory integration (dS/dt = delta_t), and comparative-statics sweeps over the model parameters.

Not ported from a standalone MATLAB HSF toolbox library function -- no hsf/fn_delta_t.m file exists in hfs-archive. The identical model (fn_lambda_t, fn_mu_long_t, fn_mu_short_t, fn_delta_t) does appear as local (embedded) functions inside HSF/16. Exercise Solution/chap16_chap5_exercise2.m, of which HSF-Notebooks chapter 16c is itself a direct Python port; promoted here from that notebook.

Source code in src/quanttoolbox/sustainable_finance/ecology.py 
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
def fn_delta_t(
    S_t: np.ndarray | float,
    lambda_0: float,
    beta1: float,
    mu_s: float,
    beta2: float,
    beta3: float,
    beta4: float,
    S_s: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Net species-richness growth rate `delta(S) = lambda(S) - mu(S)` of a
    colonization-extinction equilibrium model (the theory of island
    biogeography, TIB), where `mu(S) = mu^long(S) + mu^short(S)`:

    - `lambda(S)`: colonization rate, `=lambda_0` at `S=0`, decaying to
      `0` at the saturation richness `S=S_s`.
    - `mu^long(S)`: long-term extinction rate, `=0` at `S=0`, rising to
      `mu_s` at `S=S_s`.
    - `mu^short(S) = beta3*lambda(S)*exp(-beta4*S)`: a short-term
      extinction component proportional to the colonization rate
      (recently-arrived species are more extinction-prone).

    `S_t` is clipped to `S_t = max(0, S_t)` before evaluation, so the
    ODE/root-finder driving `dS/dt = delta(S)` never evaluates the model
    at a negative richness -- this clip protects only this one function
    evaluation, not the ODE trajectory as a whole (an unclipped `S(t)` can
    still run negative between evaluations if `delta(S)` stays negative
    all the way down).

    Returns `(delta_t, lambda_t, mu_long_t, mu_short_t, mu_t)`. Used for
    equilibrium root-finding (`delta_t == 0`), ODE trajectory integration
    (`dS/dt = delta_t`), and comparative-statics sweeps over the model
    parameters.

    Not ported from a standalone MATLAB HSF toolbox library function -- no
    `hsf/fn_delta_t.m` file exists in `hfs-archive`. The identical model
    (`fn_lambda_t`, `fn_mu_long_t`, `fn_mu_short_t`, `fn_delta_t`) does
    appear as local (embedded) functions inside `HSF/16. Exercise
    Solution/chap16_chap5_exercise2.m`, of which HSF-Notebooks chapter 16c
    is itself a direct Python port; promoted here from that notebook.
    """
    S_t = np.maximum(0.0, np.asarray(S_t, dtype=float))
    lambda_t = _fn_lambda_t(S_t, lambda_0, beta1, S_s)
    mu_long_t = _fn_mu_long_t(S_t, mu_s, beta2, S_s)
    mu_short_t = _fn_mu_short_t(S_t, lambda_t, beta3, beta4, S_s)
    mu_t = mu_long_t + mu_short_t
    delta_t = lambda_t - mu_t
    return delta_t, lambda_t, mu_long_t, mu_short_t, mu_t

hurlbert(n_i, m, method=1)

Hurlbert's rarefaction estimator: expected number of species present in a random sample of m individuals drawn (without replacement) from a community with per-species counts n_i.

method=2 uses a log-gamma formulation (via scipy.special.gammaln) for numerical stability at large sample sizes; method=1 (default) computes the binomial-coefficient ratio directly.

Original: hsf/hurlbert.m

Source code in src/quanttoolbox/sustainable_finance/ecology.py 
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
def hurlbert(n_i: np.ndarray, m: int, method: int = 1) -> float:
    """Hurlbert's rarefaction estimator: expected number of species present
    in a random sample of `m` individuals drawn (without replacement) from
    a community with per-species counts `n_i`.

    `method=2` uses a log-gamma formulation (via `scipy.special.gammaln`)
    for numerical stability at large sample sizes; `method=1` (default)
    computes the binomial-coefficient ratio directly.

    Original: hsf/hurlbert.m
    """
    n_i = np.asarray(n_i, dtype=float)
    n = np.sum(n_i)

    sac = 0.0
    for n_s in n_i:
        if (n - n_s) >= m:
            if method == 2:
                log_q = (
                    gammaln(n - n_s + 1)
                    + gammaln(n - m + 1)
                    - (gammaln(n + 1) + gammaln(n - n_s - m + 1))
                )
                q = np.exp(log_q)
            else:
                q = comb(n - n_s, m) / comb(n, m)
            sac += 1.0 - q
        else:
            sac += 1.0

    return float(sac)

species_abundance_distribution(n_i, breaks=None)

Bin per-species individual counts n_i into a species-abundance histogram: s[k] species fall into abundance class j[k].

  • breaks=None (default): one class per distinct count value in n_i.
  • breaks="octave": Preston's log2 "octave" classes, class k (k=0,1,2,...) covering 2**k <= n < 2**(k+1) (so k=0 is n=1; k=1 is n in {2,3}; k=2 is n in {4,...,7}; etc.), assigned via k = floor(log2(n)).
  • breaks=<array>: custom upper breakpoints; classes are the numpy.histogram bins of n_i over [0, *breaks] (all classes half-open [lower, upper) except the last, which is closed at both ends); j[k] is the midpoint of class k's range.

Returns (j, s, breaks) -- breaks is None for the default mode.

Original: hsf/species_abundance_distribution.m

Note: the "octave" and <array> modes' bucket-assignment rule below was corrected to this numpy-native np.histogram/floor(log2(n)) convention from an earlier, non-numpy-native inclusive rule (each class taken as (prev_break, this_break], closed at the top and open at the bottom, with the first class starting from n=1/n=0 rather than a numpy-style half-open [0, this_break)). On real survey data with non-integer abundance counts (Whittaker's Siskiyou Mountains tree survey), the two conventions disagree: HSF-Notebooks chapter 05c's independent reimplementation of this function -- built and used without reference to this toolbox version -- was cross-checked against it and found to classify some individual counts into different buckets. This toolbox's binning was corrected to match the notebook's (and numpy's own) convention exactly; see tests/sustainable_finance/test_ecology.py for the real-data cross-check.

Source code in src/quanttoolbox/sustainable_finance/ecology.py 
 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
def species_abundance_distribution(
    n_i: np.ndarray, breaks: np.ndarray | str | None = None
) -> tuple[np.ndarray, np.ndarray, np.ndarray | None]:
    """Bin per-species individual counts `n_i` into a species-abundance
    histogram: `s[k]` species fall into abundance class `j[k]`.

    - `breaks=None` (default): one class per distinct count value in `n_i`.
    - `breaks="octave"`: Preston's log2 "octave" classes, class `k`
      (`k=0,1,2,...`) covering `2**k <= n < 2**(k+1)` (so `k=0` is `n=1`;
      `k=1` is `n` in `{2,3}`; `k=2` is `n` in `{4,...,7}`; etc.), assigned
      via `k = floor(log2(n))`.
    - `breaks=<array>`: custom upper breakpoints; classes are the
      `numpy.histogram` bins of `n_i` over `[0, *breaks]` (all classes
      half-open `[lower, upper)` except the last, which is closed at both
      ends); `j[k]` is the midpoint of class `k`'s range.

    Returns `(j, s, breaks)` -- `breaks` is `None` for the default mode.

    Original: hsf/species_abundance_distribution.m

    Note: the `"octave"` and `<array>` modes' bucket-assignment rule below
    was corrected to this numpy-native `np.histogram`/`floor(log2(n))`
    convention from an earlier, non-numpy-native inclusive rule (each
    class taken as `(prev_break, this_break]`, closed at the top and open
    at the bottom, with the first class starting from `n=1`/`n=0` rather
    than a numpy-style half-open `[0, this_break)`). On real survey data
    with non-integer abundance counts (Whittaker's Siskiyou Mountains tree
    survey), the two conventions disagree: HSF-Notebooks chapter 05c's
    independent reimplementation of this function -- built and used
    without reference to this toolbox version -- was cross-checked against
    it and found to classify some individual counts into different
    buckets. This toolbox's binning was corrected to match the notebook's
    (and numpy's own) convention exactly; see
    `tests/sustainable_finance/test_ecology.py` for the real-data
    cross-check.
    """
    n_i = np.asarray(n_i, dtype=float)

    if breaks is None:
        j = np.unique(n_i)
        s = np.array([np.sum(n_i == jj) for jj in j], dtype=float)
        return j, s, None

    if isinstance(breaks, str):
        if breaks.lower() != "octave":
            raise ValueError("breaks must be None, 'octave', or an array of breakpoints")
        k = np.floor(np.log2(n_i)).astype(int)
        k_max = int(k.max())
        octaves = np.arange(0, k_max + 1)
        s = np.array([np.sum(k == kk) for kk in octaves], dtype=float)
        octave_breaks = 2.0 ** (octaves + 1)
        return octaves.astype(float), s, octave_breaks

    breaks = np.asarray(breaks, dtype=float)
    edges = np.concatenate([[0.0], breaks])
    s, _ = np.histogram(n_i, bins=edges)
    j = 0.5 * (edges[:-1] + edges[1:])
    return j, s.astype(float), breaks

species_area_relationship(n_i, s_j, area_total, area)

Expected number of species found in a sub-area area out of a total surveyed area area_total, given either per-species individual counts n_i (s_j=None), or a species-abundance histogram s_j (number of species with exactly j = 1, 2, ... individuals).

Original: hsf/species_area_relationship.m

Source code in src/quanttoolbox/sustainable_finance/ecology.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def species_area_relationship(
    n_i: np.ndarray, s_j: np.ndarray | None, area_total: float, area: np.ndarray | float
) -> np.ndarray:
    """Expected number of species found in a sub-area `area` out of a total
    surveyed area `area_total`, given either per-species individual counts
    `n_i` (`s_j=None`), or a species-abundance histogram `s_j` (number of
    species with exactly `j` = 1, 2, ... individuals).

    Original: hsf/species_area_relationship.m
    """
    area = np.asarray(area, dtype=float)
    ratio = 1.0 - area / area_total

    if s_j is None:
        n_i = np.asarray(n_i, dtype=float)
        total_species = n_i.shape[0]
        term = np.sum(ratio[..., None] ** n_i, axis=-1)
        return total_species - term

    s_j = np.asarray(s_j, dtype=float)
    total_species = np.sum(s_j)
    j = np.arange(1, s_j.shape[0] + 1)
    term = np.sum(s_j * ratio[..., None] ** j, axis=-1)
    return total_species - term

sustainable_finance.entropy

Python alternatives

Keep — Shannon entropy/mutual-information decomposition (shannon_entropy, shannon_entropy_markov_chain) and the Israel-Rosenthal-Wei (2001) Markov-generator-matrix repair technique (estimate_markov_generator), used here for rating-migration-generator estimation. No general-purpose equivalent found.

quanttoolbox.sustainable_finance.entropy

Shannon-entropy diversity/dependence measures, and the Israel-Rosenthal-Wei (2001) technique for repairing an estimated Markov generator matrix that fails to satisfy the generator constraints (non-negative off-diagonal entries, zero row sums) because it was estimated from discretely-observed transition data.

Ported from HSF toolbox hsf/{shannon_entropy, shannon_entropy_markov_chain,estimate_markov_generator}.m.

Translation notes:

  • MATLAB's missrv(x, v) (replace non-finite/flagged entries of x with v) is used throughout the originals to implement the 0 * log(0) = 0 convention; translated as scipy.special.xlogy(p, p), which computes p * log(p) and returns exactly 0.0 at p == 0 without ever evaluating log(0) (so, unlike np.where(p > 0, p * np.log(p), 0.0), it raises no divide-by-zero warning).
  • shannon_entropy_markov_chain.m approximates the Markov chain's stationary distribution via expm(Lambda * 1000) (a long-horizon transition-probability matrix, every row of which has converged to the stationary distribution) -- kept as-is via scipy.linalg.expm.

  • estimate_markov_generator.m implements Israel, Rosenthal & Wei (2001)'s two generator-repair methods for a possibly-invalid estimated generator Lambda (e.g. from Lambda = logm(P) / dt on an empirical transition matrix P, which need not itself be a valid generator):

  • Lambda1 ("diagonal adjustment"): zero out negative off-diagonal entries and absorb the removed mass into the diagonal. This exactly preserves each row's original sum (min(x, 0) + max(x, 0) = x), so Lambda1's row sums equal Lambda's row sums (zero, if Lambda already had zero row sums as intended).

  • Lambda2 ("proportional redistribution"): redistribute each row's negative mass proportionally across that row's positive off-diagonal entries, leaving rows with no positive off-diagonal mass (g[i] == 0) unchanged.

Both are standard generator-matrix regularizations in credit-rating transition-intensity estimation (Roncalli's hsf toolbox applies this to rating-migration generators); the row-sum-preservation property of Lambda1 is verified in this module's tests.

estimate_markov_generator(lam)

Repair a possibly-invalid estimated Markov generator lam (e.g. one with negative off-diagonal entries) into two valid generators (non-negative off-diagonal, zero row sums), via Israel, Rosenthal & Wei (2001)'s two methods.

Returns (Lambda1, Lambda2) -- see the module docstring for the method descriptions.

Original: hsf/estimate_markov_generator.m

Source code in src/quanttoolbox/sustainable_finance/entropy.py 
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
def estimate_markov_generator(lam: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Repair a possibly-invalid estimated Markov generator `lam` (e.g. one
    with negative off-diagonal entries) into two valid generators
    (non-negative off-diagonal, zero row sums), via Israel, Rosenthal & Wei
    (2001)'s two methods.

    Returns ``(Lambda1, Lambda2)`` -- see the module docstring for the
    method descriptions.

    Original: hsf/estimate_markov_generator.m
    """
    lam = np.asarray(lam, dtype=float)
    k = lam.shape[0]
    off_diag_mask = ~np.eye(k, dtype=bool)

    neg_off_diag = np.where(off_diag_mask, np.minimum(lam, 0.0), 0.0)
    pos_off_diag = np.where(off_diag_mask, np.maximum(lam, 0.0), 0.0)
    row_neg_sum = neg_off_diag.sum(axis=1)
    row_pos_sum = pos_off_diag.sum(axis=1)

    new_diag = np.diag(lam) + row_neg_sum
    lambda1 = pos_off_diag.copy()
    np.fill_diagonal(lambda1, new_diag)

    g = np.abs(np.diag(lam)) + row_pos_sum
    b = -row_neg_sum
    g_col = g[:, None]
    safe_g = np.where(g_col > 0, g_col, 1.0)
    lambda2 = lam - b[:, None] * np.abs(lam) / safe_g
    lambda2 = np.where(g_col > 0, lambda2, lam)

    neg_off_diag_mask = off_diag_mask & (lam < 0)
    lambda2 = np.where(neg_off_diag_mask, 0.0, lambda2)

    return lambda1, lambda2

shannon_entropy(p_xy)

Shannon entropy/mutual-information decomposition of a discrete distribution: p_xy may be a 1D probability vector (single variable, in which case I_Y = I_XY = 0 and I_X_Y = I_X), or a 2D joint probability matrix (rows = X, columns = Y).

Returns (I_X, I_Y, I_XY, I_X_Y): the marginal entropy of X, the marginal entropy of Y, the mutual information between X and Y, and the joint entropy of (X, Y).

Original: hsf/shannon_entropy.m

Source code in src/quanttoolbox/sustainable_finance/entropy.py 
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
def shannon_entropy(p_xy: np.ndarray) -> tuple[float, float, float, float]:
    """Shannon entropy/mutual-information decomposition of a discrete
    distribution: `p_xy` may be a 1D probability vector (single variable,
    in which case `I_Y = I_XY = 0` and `I_X_Y = I_X`), or a 2D joint
    probability matrix (rows = X, columns = Y).

    Returns ``(I_X, I_Y, I_XY, I_X_Y)``: the marginal entropy of X, the
    marginal entropy of Y, the mutual information between X and Y, and the
    joint entropy of (X, Y).

    Original: hsf/shannon_entropy.m
    """
    p_xy = np.asarray(p_xy, dtype=float)

    if p_xy.ndim == 1:
        p_x = p_xy
        i_x = -float(np.sum(xlogy(p_x, p_x)))
        return i_x, 0.0, 0.0, i_x

    p_x = np.sum(p_xy, axis=1)
    p_y = np.sum(p_xy, axis=0)

    i_x = -float(np.sum(xlogy(p_x, p_x)))
    i_y = -float(np.sum(xlogy(p_y, p_y)))
    i_xy_joint = -float(np.sum(xlogy(p_xy, p_xy)))

    i_xy = i_x + i_y - i_xy_joint
    return i_x, i_y, i_xy, i_xy_joint

shannon_entropy_markov_chain(lambda_, t)

Shannon entropy/mutual-information decomposition of the joint distribution of a continuous-time Markov chain's state at times 0 and t, with generator lambda_ and stationary (long-run, t -> inf) distribution approximated via expm(lambda_ * 1000).

t may be an array of horizons; returns one value per horizon for each of (I_X, I_Y, I_XY, I_X_Y).

Original: hsf/shannon_entropy_markov_chain.m

Source code in src/quanttoolbox/sustainable_finance/entropy.py 
 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
def shannon_entropy_markov_chain(
    lambda_: np.ndarray, t: np.ndarray | float
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Shannon entropy/mutual-information decomposition of the joint
    distribution of a continuous-time Markov chain's state at times 0 and
    `t`, with generator `lambda_` and stationary (long-run, ``t -> inf``)
    distribution approximated via ``expm(lambda_ * 1000)``.

    `t` may be an array of horizons; returns one value per horizon for
    each of ``(I_X, I_Y, I_XY, I_X_Y)``.

    Original: hsf/shannon_entropy_markov_chain.m
    """
    lambda_ = np.asarray(lambda_, dtype=float)
    t_arr = np.atleast_1d(np.asarray(t, dtype=float))
    n = lambda_.shape[0]

    p_inf = expm(lambda_ * 1000.0)
    pi = p_inf[0, :]

    n_t = t_arr.shape[0]
    i_x = np.zeros(n_t)
    i_y = np.zeros(n_t)
    i_xy = np.zeros(n_t)
    i_xy_joint = np.zeros(n_t)

    for it in range(n_t):
        if t_arr[it] == 0.0:
            p_t = np.eye(n)
        else:
            p_t = expm(lambda_ * t_arr[it])
        p_xy = pi[:, None] * p_t
        i_x[it], i_y[it], i_xy[it], i_xy_joint[it] = shannon_entropy(p_xy)

    return i_x, i_y, i_xy, i_xy_joint