"""GMM fast-path kernel codegen -- per-topology numba kernel generation, compile, cache.
For the headline-only fast path (``full=False``) the multi-state kernel is
GENERATED per state-machine topology (state count, edge topology, lump-sum
flags, which states pay premium / benefit are baked into the emitted Python
source) so numba's compiled inner loop has only scalar arithmetic, no array
indirection. Generated source is cached on disk by a hash of its text; numba
caches the native code beside it -- compile-once per topology across processes.
``_scalar_kernel`` is the directly-defined scalar fallback. This is
metaprogramming (code that writes kernels), kept out of the measurement logic in
``_engine``; it imports only stdlib + numba + numpy (no domain objects).
"""
from __future__ import annotations
import hashlib
import importlib.util
import os
import sys
from pathlib import Path
import numpy as np
from numba import njit, prange
# ---------------------------------------------------------------------------
# Codegen specialisation -- the multi-state fast kernel
# ---------------------------------------------------------------------------
#
# The multi-state CPU fast-path (full=False) kernel is generated per Model: the state
# count, the full edge topology, the lump-sum flags, and which states pay
# premium or a benefit all become part of the generated Python source so
# numba's compiled inner loop has no array indirections left, only scalar
# arithmetic on register-resident occupancy. One source file lives on disk
# per unique topology (see _cache_dir), and numba caches its native
# code next to it -- compile-once per topology across Python processes.
#
# Earlier iterations of this engine kept a Markov-only closure factory plus
# hand-unrolled n_states=2 and n_states=3 kernels alongside the codegen
# path. They became dead code once the codegen path was extended to all
# n_states>=2 (commit e5e8e83) and have been removed; the git history
# carries them for anyone who wants to read the earlier shape.
def _markov_kernel_source(n_states, edge_from, edge_to, edge_lump_sum,
state_pays_premium, state_pays_benefit,
premium_term_to=None,
use_morbidity=True, use_annuity=True,
use_disability=True, use_lae=True,
use_surrender=True,
surrender_is_amount=False) -> str:
"""Generate the Python source of a fully-specialised fast kernel.
All structural parameters (n_states, edge topology, lump-sum flags,
premium- and benefit-paying states) are baked into the source as
literals. The returned text is intended for ``exec`` in a namespace
that exposes ``np``, ``njit`` and ``prange``.
The ``use_*`` flags specialise the inner loop at source-generation time:
a cash-flow stream the portfolio does not use (morbidity, annuity,
disability income, LAE, surrender) has its per-cell lines omitted from
the generated source entirely -- no runtime branch, the absent term is
simply not there. Each flag widens the codegen cache key, so a book that
does use the stream gets its own kernel. Skipped terms are zero, so the
result is identical to the all-streams kernel.
"""
n_edges = len(edge_from)
edge_from = [int(x) for x in edge_from]
edge_to = [int(x) for x in edge_to]
edge_lump_sum = [bool(x) for x in edge_lump_sum]
state_pays_premium = [bool(x) for x in state_pays_premium]
state_pays_benefit = [bool(x) for x in state_pays_benefit]
# Calendar-keyed (at_premium_term) deterministic transitions: per-state
# destination (-1 none, -2 exit, >=0 dest state). The movers relabel
# occupancy at the model point's premium_term -- emitted inside every
# occupancy-advance pass so each pass's in-force trajectory switches at
# premium_term consistently. A self destination (dto == s) is a no-op.
if premium_term_to is None:
premium_term_to = [-1] * n_states
else:
premium_term_to = [int(x) for x in premium_term_to]
_movers = [(s, premium_term_to[s]) for s in range(n_states)
if premium_term_to[s] != -1 and premium_term_to[s] != s]
sum_all = " + ".join(f"occ_{i}" for i in range(n_states))
sum_prem = " + ".join(f"occ_{i}" for i in range(n_states)
if state_pays_premium[i]) or "0.0"
sum_ben = " + ".join(f"occ_{i}" for i in range(n_states)
if state_pays_benefit[i]) or "0.0"
# State-machine lapse exits: occupancy on each state times that state's own
# lapse rate, so surrender follows the actual lapse (a non-lapsing state
# contributes nothing; a paid-up state lapses at its own rate). For a single
# active state this equals ``inforce_t * lapse``, the historical formula.
sum_lapse = " + ".join(
f"occ_{i} * state_lapse[{i}, sx, age_idx, year]" for i in range(n_states))
L: list[str] = []
def line(indent: int, text: str) -> None:
L.append(" " * indent + text)
# Module-level prologue. The generated text is written to a real .py
# file under the on-disk cache so numba's @njit(cache=True) can anchor
# its own compile cache to that file -- without a source file numba's
# cache silently no-ops and every Python process pays the JIT cost.
line(0, '"""Auto-generated by fastcashflow._measurement.gmm.codegen.'
'_markov_kernel_source -- do not edit."""')
line(0, "import numpy as np")
line(0, "from numba import njit, prange")
line(0, "")
line(0, "")
def emit_init(indent: int) -> None:
for i in range(n_states):
line(indent, f"occ_{i} = 0.0")
line(indent, "if ss == 0:")
line(indent + 4, "occ_0 = cnt")
for i in range(1, n_states):
line(indent, f"elif ss == {i}:")
line(indent + 4, f"occ_{i} = cnt")
def emit_edge_step(indent: int, scale: str = "",
include_lump: bool = True) -> None:
for i in range(n_states):
line(indent, f"occ_next_{i} = 0.0")
for e in range(n_edges):
ef, et, ls = edge_from[e], edge_to[e], edge_lump_sum[e]
line(indent,
f"flow_{e} = occ_{ef}{scale} * edge_prob[sx, age_idx, year, {e}]")
line(indent, f"occ_next_{et} += flow_{e}")
if include_lump and ls:
line(indent,
f"pv_disability += flow_{e} * disability_benefit[mp] * discount_factor_mid_t")
# Calendar-keyed move (active -> paid-up at premium_term): snapshot each
# mover's pre-move occupancy, then relabel -- order-independent, safe for
# chained / self destinations (mirrors the detailed kernel).
if _movers:
line(indent, "if premium_term > 0 and t + 1 == premium_term:")
for s, _dto in _movers:
line(indent + 4, f"_pt_src_{s} = occ_next_{s}")
for s, dto in _movers:
line(indent + 4, f"occ_next_{s} -= _pt_src_{s}")
if dto >= 0:
line(indent + 4, f"occ_next_{dto} += _pt_src_{s}")
# dto == -2: exit the in-force set (not added anywhere)
for i in range(n_states):
line(indent, f"occ_{i} = occ_next_{i}")
line(0, "@njit(parallel=True, cache=True)")
line(0, "def kernel(edge_from, edge_to, edge_prob, edge_lump_sum,")
line(0, " state_pays_premium, state_pays_benefit, start_state, "
"issue_index, sex,")
line(0, " term_months, contract_boundary_months, count, premium,")
line(0, " premium_term_months, premium_frequency_months, "
"annuity_frequency_months,")
line(0, " coverage_index, coverage_amount, coverage_offset, coverage_rates, "
"premium_factor, annuity_factor, coverage_risk,")
line(0, " coverage_is_diagnosis, maturity_benefit, "
"annuity_payment,")
line(0, " disability_income, disability_benefit,")
line(0, " acquisition_premium, acquisition_per_policy, maintenance_premium,")
line(0, " maintenance_per_policy, lae,")
line(0, " discount_factor_bom, discount_factor_mid, mortality_factor,")
line(0, " morbidity_factor, longevity_factor, "
"disability_factor,")
line(0, " coverage_waiting, coverage_reduction_end, "
"coverage_reduction_factor,")
line(0, " lapse_monthly, state_lapse, surrender_curve, surrender_base):")
line(4, "n_mp = issue_index.shape[0]")
line(4, "bel = np.empty(n_mp)")
line(4, "ra = np.empty(n_mp)")
line(4, "csm = np.empty(n_mp)")
line(4, "loss_component = np.empty(n_mp)")
line(4, "for mp in prange(n_mp):")
line(8, "term = term_months[mp]")
line(8, "boundary = contract_boundary_months[mp]")
line(8, "premium_term = premium_term_months[mp]")
line(8, "prem_freq = premium_frequency_months[mp]")
line(8, "ann_freq = annuity_frequency_months[mp]")
line(8, "age_idx = issue_index[mp]")
line(8, "sx = sex[mp]")
line(8, "cnt = count[mp]")
line(8, "prem = premium[mp]")
line(8, "annuity = annuity_payment[mp]")
line(8, "c_start = coverage_offset[mp]")
line(8, "c_end = coverage_offset[mp + 1]")
line(8, "ss = start_state[mp]")
line(8, "pf = 1.0")
line(8, "ann_prem = prem * 12.0 / prem_freq")
emit_init(8)
line(8, "pv_mortality = 0.0")
line(8, "pv_morbidity = 0.0")
line(8, "pv_disability = 0.0")
line(8, "pv_premium = 0.0")
line(8, "pv_expense = 0.0")
line(8, "pv_annuity = 0.0")
line(8, "pv_surrender = 0.0")
line(8, "cum_premium = 0.0")
line(8, "last_year = -1")
line(8, "claim_rate = 0.0")
line(8, "morb_rate = 0.0")
line(8, "prem_due = 0")
line(8, "ann_due = 0")
line(8, "prem_left = premium_term")
# Main t loop
line(8, "for t in range(boundary):")
line(12, "year = t // 12")
line(12, "if year != last_year:")
line(16, "claim_rate = 0.0")
line(16, "morb_rate = 0.0")
line(16, "for k in range(c_start, c_end):")
line(20, "cov_idx = coverage_index[k]")
line(20, "if coverage_is_diagnosis[cov_idx]:")
line(24, "continue")
line(20, "if coverage_waiting[k] != 0 or coverage_reduction_end[k] != 0:")
line(24, "continue")
line(20, "rate = coverage_rates[cov_idx, sx, age_idx, year] * coverage_amount[k]")
line(20, "if coverage_risk[cov_idx] == 0:")
line(24, "claim_rate += rate")
line(20, "else:")
line(24, "morb_rate += rate")
line(16, "last_year = year")
line(16, "pf = premium_factor[sx, age_idx, year]")
line(16, "ann_prem = prem * pf * 12.0 / prem_freq")
line(12, f"inforce_t = {sum_all}")
line(12, f"prem_occ = {sum_prem}")
if use_disability:
line(12, f"benefit_occ = {sum_ben}")
line(12, "discount_factor_bom_t = discount_factor_bom[t]")
line(12, "discount_factor_mid_t = discount_factor_mid[t]")
line(12, "if prem_due == 0 and prem_left > 0:")
line(16, "level = prem_occ * prem * pf")
line(16, "prem_due = prem_freq - 1")
line(12, "else:")
line(16, "level = 0.0")
line(16, "prem_due -= 1")
line(12, "prem_left -= 1")
line(12, "pv_premium += level * discount_factor_bom_t")
if use_surrender and not surrender_is_amount:
line(12, "cum_premium += level")
line(12, "pv_mortality += inforce_t * claim_rate * discount_factor_mid_t")
if use_morbidity:
line(12, "pv_morbidity += inforce_t * morb_rate * discount_factor_mid_t")
if use_annuity:
line(12, "if ann_due == 0:")
line(16, "pv_annuity += inforce_t * annuity * annuity_factor[sx, age_idx, year] * discount_factor_bom_t")
line(16, "ann_due = ann_freq - 1")
line(12, "else:")
line(16, "ann_due -= 1")
if use_disability:
line(12, "pv_disability += benefit_occ * disability_income[mp] * discount_factor_mid_t")
if use_surrender:
line(12, f"lapse_flow = {sum_lapse}")
if surrender_is_amount:
# amount_per_policy / amount_per_unit: surrender_curve[t] is the
# surrender amount at duration t (per policy, or per unit of
# surrender_base[mp]); lapse_flow is the number lapsing.
# surrender_base is 1.0 for amount_per_policy.
line(12, "pv_surrender += (lapse_flow * surrender_curve[t]")
line(12, " * surrender_base[mp] * discount_factor_mid_t)")
else:
# cum_premium_factor: cum_premium aggregates inforce * premium; the
# effective lapse fraction is lapse_flow / inforce_t (the raw rate for a
# single state). inforce_t carries the count, so dividing it out avoids a
# cnt^2 scaling.
line(12, "eff_lapse = lapse_flow / inforce_t if inforce_t > 0.0 else 0.0")
line(12, "pv_surrender += (eff_lapse")
line(12, " * cum_premium * surrender_curve[t] * discount_factor_mid_t)")
line(12, "acquisition_expense = cnt * (acquisition_premium * ann_prem + acquisition_per_policy) if t == 0 else 0.0")
line(12, "maintenance_premium_expense = inforce_t * maintenance_premium[t] * ann_prem / 12.0 if t < premium_term else 0.0")
line(12, "maintenance_per_policy_expense = inforce_t * maintenance_per_policy[t]")
if use_lae:
line(12, "lae_expense = lae[t] * "
"inforce_t * (claim_rate + morb_rate)")
line(12, "pv_expense += (acquisition_expense + maintenance_premium_expense + maintenance_per_policy_expense + lae_expense) * discount_factor_mid_t")
else:
line(12, "pv_expense += (acquisition_expense + maintenance_premium_expense + maintenance_per_policy_expense) * discount_factor_mid_t")
emit_edge_step(12, scale="", include_lump=True)
line(8, f"total = {sum_all}")
line(8, "pm = (total * maturity_benefit[mp] * discount_factor_bom[boundary]) "
"if boundary == term else 0.0")
# Coverage-rule pass
line(8, "for k in range(c_start, c_end):")
line(12, "cov_idx = coverage_index[k]")
line(12, "if coverage_is_diagnosis[cov_idx]:")
line(16, "continue")
line(12, "wait = coverage_waiting[k]")
line(12, "red_end = coverage_reduction_end[k]")
line(12, "if wait == 0 and red_end == 0:")
line(16, "continue")
line(12, "benefit = coverage_amount[k]")
line(12, "red_factor = coverage_reduction_factor[k]")
line(12, "mortality_risk = coverage_risk[cov_idx] == 0")
emit_init(12)
line(12, "for t in range(boundary):")
line(16, "year = t // 12")
line(16, "if t >= wait:")
line(20, "mult = red_factor if t < red_end else 1.0")
line(20, f"inforce_t = {sum_all}")
line(20, "contrib = (inforce_t * coverage_rates[cov_idx, sx, age_idx, year]")
line(20, " * benefit * mult * discount_factor_mid[t])")
line(20, "if mortality_risk:")
line(24, "pv_mortality += contrib")
line(20, "else:")
line(24, "pv_morbidity += contrib")
emit_edge_step(16, scale="", include_lump=False)
# Diagnosis pass
line(8, "for k in range(c_start, c_end):")
line(12, "cov_idx = coverage_index[k]")
line(12, "if not coverage_is_diagnosis[cov_idx]:")
line(16, "continue")
line(12, "benefit = coverage_amount[k]")
line(12, "wait = coverage_waiting[k]")
line(12, "red_end = coverage_reduction_end[k]")
line(12, "red_factor = coverage_reduction_factor[k]")
emit_init(12)
line(12, "d_year = -1")
line(12, "d_rate = 0.0")
line(12, "for t in range(boundary):")
line(16, "year = t // 12")
line(16, "if year != d_year:")
line(20, "d_rate = coverage_rates[cov_idx, sx, age_idx, year]")
line(20, "d_year = year")
line(16, "if t >= wait:")
line(20, "mult = red_factor if t < red_end else 1.0")
line(20, f"healthy = {sum_all}")
line(20, "pv_morbidity += healthy * d_rate * benefit * mult * discount_factor_mid[t]")
line(16, "undiagnosed = 1.0 - d_rate")
emit_edge_step(16, scale=" * undiagnosed", include_lump=False)
line(8, "bel_mp = (pv_mortality + pv_morbidity + pv_disability + pm")
line(8, " + pv_annuity + pv_expense + pv_surrender - pv_premium)")
line(8, "ra_mp = (mortality_factor * pv_mortality + morbidity_factor * pv_morbidity")
line(8, " + disability_factor * pv_disability "
"+ longevity_factor * (pm + pv_annuity))")
line(8, "fcf = bel_mp + ra_mp")
line(8, "bel[mp] = bel_mp")
line(8, "ra[mp] = ra_mp")
line(8, "csm[mp] = max(0.0, -fcf)")
line(8, "loss_component[mp] = max(0.0, fcf)")
line(4, "return bel, ra, csm, loss_component")
return "\n".join(L)
_MARKOV_KERNEL_CACHE: dict = {}
def _cache_dir() -> Path:
"""Return the on-disk directory holding generated kernel source files.
Honours ``XDG_CACHE_HOME`` and falls back to ``~/.cache``; the
fastcashflow-private subdirectory is created on demand.
"""
base = os.environ.get("XDG_CACHE_HOME")
root = Path(base) if base else Path.home() / ".cache"
cache = root / "fastcashflow" / "codegen"
cache.mkdir(parents=True, exist_ok=True)
return cache
def _atomic_write_text(path: Path, content: str) -> None:
"""Write ``content`` to ``path`` atomically.
Two processes racing on the same topology hash produce byte-identical
sources (codegen is deterministic), so the final file content is the
same regardless of who wins -- the atomic replace just avoids a
partially-written file being observed.
"""
tmp = path.with_name(f"{path.name}.{os.getpid()}.tmp")
tmp.write_text(content)
os.replace(tmp, path)
def _get_markov_kernel(n_states, edge_from, edge_to, edge_lump_sum,
state_pays_premium, state_pays_benefit,
premium_term_to=None,
use_morbidity=True, use_annuity=True,
use_disability=True, use_lae=True,
use_surrender=True, surrender_is_amount=False):
"""Return a codegen-specialised kernel for the given state machine.
Two-level cache:
* Process-local dict keyed on the topology -- a repeat lookup in the
same Python process returns the already-imported function with no
filesystem touch.
* On-disk ``.py`` file under the codegen cache directory, named by a
hash of the generated source. The first call for a given topology
generates the source and writes the file atomically; ``@njit(cache
=True)`` inside the generated module then persists the compiled
bytecode in its ``__pycache__`` so subsequent Python processes pay
no JIT cost. Hashing the *source* automatically invalidates the
cache when the codegen logic itself changes.
"""
if premium_term_to is None:
premium_term_to = [-1] * int(n_states)
key = (
int(n_states),
tuple(int(x) for x in edge_from),
tuple(int(x) for x in edge_to),
tuple(bool(x) for x in edge_lump_sum),
tuple(bool(x) for x in state_pays_premium),
tuple(bool(x) for x in state_pays_benefit),
tuple(int(x) for x in premium_term_to),
bool(use_morbidity), bool(use_annuity), bool(use_disability),
bool(use_lae), bool(use_surrender), bool(surrender_is_amount),
)
cached = _MARKOV_KERNEL_CACHE.get(key)
if cached is not None:
return cached
src = _markov_kernel_source(
n_states, edge_from, edge_to, edge_lump_sum,
state_pays_premium, state_pays_benefit, premium_term_to=premium_term_to,
use_morbidity=use_morbidity, use_annuity=use_annuity,
use_disability=use_disability, use_lae=use_lae,
use_surrender=use_surrender, surrender_is_amount=surrender_is_amount,
)
digest = hashlib.sha256(src.encode("utf-8")).hexdigest()[:16]
cache_path = _cache_dir() / f"fast_kernel_{digest}.py"
if not cache_path.exists():
_atomic_write_text(cache_path, src)
module_name = f"_fastcashflow_codegen_{digest}"
module = sys.modules.get(module_name)
if module is None:
spec = importlib.util.spec_from_file_location(module_name, cache_path)
module = importlib.util.module_from_spec(spec)
sys.modules[module_name] = module
spec.loader.exec_module(module)
kernel = module.kernel
_MARKOV_KERNEL_CACHE[key] = kernel
return kernel
[문서]
def clear_codegen_cache(*, prune_older_than_days: float | None = None) -> int:
"""Remove generated kernel sources from the on-disk codegen cache.
Each unique state-machine topology persists a ``fast_kernel_*.py``
(or ``fast_kernel_sm_*.py``) source file plus its ``__pycache__``
sidecars. Over a session that explores many topologies these
accumulate. This helper sweeps them.
Parameters
----------
prune_older_than_days
Only drop files whose mtime is older than this many days. ``None``
(the default) drops everything.
Returns
-------
int
Number of files removed (sources and ``__pycache__`` artefacts).
Notes
-----
Disk-only: modules already imported in the current process continue to
work via :data:`sys.modules`. Codegen for a still-needed topology will
transparently regenerate on the next call.
"""
import time
cache_dir = _cache_dir()
cutoff = (
time.time() - 86400.0 * prune_older_than_days
if prune_older_than_days is not None else None
)
removed = 0
for entry in cache_dir.iterdir():
if entry.is_file() and entry.name.startswith("fast_kernel_"):
if cutoff is None or entry.stat().st_mtime < cutoff:
entry.unlink()
removed += 1
elif entry.is_dir() and entry.name == "__pycache__":
for sub in entry.iterdir():
if cutoff is None or sub.stat().st_mtime < cutoff:
sub.unlink()
removed += 1
return removed
# ---------------------------------------------------------------------------
# Semi-Markov codegen (Phase (c) prototype -- cancer reincidence)
# ---------------------------------------------------------------------------
#
# When the Model declares any state with ``sojourn_tracking_months > 0`` the engine
# tracks per-cohort occupancy in that state -- ``occ[state, tau]`` where
# ``tau`` is the sojourn time (months since entering the state). Transitions
# marked ``sojourn_dependent=True`` then carry per-cohort rates, the natural
# way to express recovery, reincidence and exclusion-period effects.
#
# The semi-Markov codegen extends the existing codegen pattern: every cohort
# of every state becomes a scalar local (``occ_<s>_<tau>``), every edge
# unrolls per-cohort (a residual stay advances the cohort with absorbing
# semantics at the last cohort, a transient transition enters the destination
# state's cohort 0). The generated source goes through the same disk-cached
# numba compile path as the Markov codegen -- compile-once per topology.
#
# Coverage-rule and diagnosis-coverage passes are emitted as no-op stubs when
# the Model is semi-Markov: this prototype targets the cancer-reincidence
# product, which has no rule-bearing or diagnosis coverages on the model
# points themselves -- the reincidence benefit rides the transition lump-sum.
# Full coverage-rule / diagnosis support for semi-Markov is a follow-up.
def _semi_markov_kernel_source(
n_states, state_duration_max, periodic_benefit_term_months,
edge_from, edge_to, edge_lump_sum,
state_pays_premium, state_pays_benefit,
use_annuity=True, use_lae=True, use_surrender=True,
surrender_is_amount=False,
) -> str:
"""Generate the Python source of a semi-Markov-aware fast kernel.
Same structure as :func:`_markov_kernel_source` but with per-state
cohort scalars and per-cohort edge processing. State ``s`` with
``state_duration_max[s] = D`` produces ``occ_s_0, occ_s_1, ..., occ_s_{D-1}``
local scalars; the last cohort absorbs everyone with sojourn time at
least ``D - 1`` months.
"""
n_edges = len(edge_from)
edge_from = [int(x) for x in edge_from]
edge_to = [int(x) for x in edge_to]
edge_lump_sum = [bool(x) for x in edge_lump_sum]
state_pays_premium = [bool(x) for x in state_pays_premium]
state_pays_benefit = [bool(x) for x in state_pays_benefit]
D = [int(x) for x in state_duration_max]
cap = [int(x) for x in periodic_benefit_term_months]
def occ(s, tau):
return f"occ_{s}_{tau}"
def occ_next(s, tau):
return f"occ_next_{s}_{tau}"
# All cohort scalars across all states -- for sums and resets.
all_occ = [occ(s, tau) for s in range(n_states) for tau in range(D[s])]
all_next = [occ_next(s, tau) for s in range(n_states) for tau in range(D[s])]
sum_all = " + ".join(all_occ)
sum_prem = " + ".join(occ(s, tau) for s in range(n_states)
if state_pays_premium[s]
for tau in range(D[s])) or "0.0"
# ``periodic_benefit_term_months[s] > 0`` caps the paid cohorts: only sojourn
# ``tau < cap`` is paid (lives past the cap stay in force, see
# ``periodic_benefit_term_months`` on State). ``cap`` is a compile-time loop bound,
# so the hot kernel carries no extra branch. Validation guarantees
# ``cap < D``, so ``min`` is belt-and-suspenders.
sum_ben = " + ".join(
occ(s, tau) for s in range(n_states)
if state_pays_benefit[s]
for tau in range(min(D[s], cap[s]) if cap[s] else D[s])) or "0.0"
# State-machine lapse exits: each state's cohort-summed occupancy times its
# own lapse rate (cohort-independent), so surrender follows the actual
# lapse. Equals ``inforce_t * lapse`` for a single lapsing state.
sum_lapse = " + ".join(
f"({' + '.join(occ(s, tau) for tau in range(D[s]))})"
f" * state_lapse[{s}, sx, age_idx, year]"
for s in range(n_states) if D[s] > 0) or "0.0"
L: list[str] = []
def line(indent: int, text: str) -> None:
L.append(" " * indent + text)
def emit_init(indent: int) -> None:
"""Initialize all cohort scalars to 0, then seat ``cnt`` based on ss."""
for s in range(n_states):
for tau in range(D[s]):
line(indent, f"{occ(s, tau)} = 0.0")
# Seating: when entering a state via seating, always cohort 0.
line(indent, "if ss == 0:")
line(indent + 4, f"{occ(0, 0)} = cnt")
for s in range(1, n_states):
line(indent, f"elif ss == {s}:")
line(indent + 4, f"{occ(s, 0)} = cnt")
def emit_edge_step(indent: int, include_lump: bool = True,
scale: str = "") -> None:
"""Emit a per-timestep advance: reset next-buffer, process every edge
across every source cohort, then copy back.
Edge semantics:
- residual (edge_from[e] == edge_to[e]): cohort tau advances to tau+1,
with cohort D-1 absorbing (long-tail).
- transient (edge_from[e] != edge_to[e]): enters destination state's
cohort 0.
- exits are not represented as edges; occupancy that doesn't go to
any next-buffer slot simply leaves the in-force set.
``scale`` (e.g. ``" * undiagnosed"``) multiplies every flow -- used by
the diagnosis-coverage pass to deplete the not-yet-diagnosed
in-force fraction along the state machine.
"""
for s in range(n_states):
for tau in range(D[s]):
line(indent, f"{occ_next(s, tau)} = 0.0")
for e in range(n_edges):
s_from = edge_from[e]
s_to = edge_to[e]
is_residual = (s_from == s_to)
ls = edge_lump_sum[e]
for tau in range(D[s_from]):
line(indent,
f"flow_{e}_{tau} = {occ(s_from, tau)}{scale} "
f"* edge_prob[sx, age_idx, year, {e}, {tau}]")
if is_residual:
# Cohort advance; cap at D-1 (absorbing).
next_tau = tau + 1 if tau + 1 < D[s_from] else D[s_from] - 1
line(indent,
f"{occ_next(s_from, next_tau)} += flow_{e}_{tau}")
else:
# Transient: enter destination's cohort 0.
line(indent,
f"{occ_next(s_to, 0)} += flow_{e}_{tau}")
if include_lump and ls:
line(indent,
f"pv_disability += flow_{e}_{tau} "
f"* disability_benefit[mp] * discount_factor_mid_t")
# Copy back.
for s in range(n_states):
for tau in range(D[s]):
line(indent, f"{occ(s, tau)} = {occ_next(s, tau)}")
# --- Prologue: imports + decorator ----------------------------------
line(0, '"""Auto-generated by fastcashflow._measurement.gmm.codegen.'
'_semi_markov_kernel_source -- do not edit."""')
line(0, "import numpy as np")
line(0, "from numba import njit, prange")
line(0, "")
line(0, "")
line(0, "@njit(parallel=True, cache=True)")
line(0, "def kernel(edge_prob, start_state, issue_index, sex,")
line(0, " term_months, contract_boundary_months, count, premium,")
line(0, " premium_term_months, premium_frequency_months, "
"annuity_frequency_months,")
line(0, " coverage_index, coverage_amount, coverage_offset, coverage_rates, "
"premium_factor, annuity_factor, coverage_risk,")
line(0, " coverage_is_diagnosis, maturity_benefit, "
"annuity_payment,")
line(0, " disability_income, disability_benefit,")
line(0, " acquisition_premium, acquisition_per_policy, maintenance_premium,")
line(0, " maintenance_per_policy, lae,")
line(0, " discount_factor_bom, discount_factor_mid, mortality_factor,")
line(0, " morbidity_factor, longevity_factor, "
"disability_factor,")
line(0, " coverage_waiting, coverage_reduction_end, "
"coverage_reduction_factor,")
line(0, " lapse_monthly, state_lapse, surrender_curve, surrender_base):")
line(4, "n_mp = issue_index.shape[0]")
line(4, "bel = np.empty(n_mp)")
line(4, "ra = np.empty(n_mp)")
line(4, "csm = np.empty(n_mp)")
line(4, "loss_component = np.empty(n_mp)")
line(4, "for mp in prange(n_mp):")
line(8, "term = term_months[mp]")
line(8, "boundary = contract_boundary_months[mp]")
line(8, "premium_term = premium_term_months[mp]")
line(8, "prem_freq = premium_frequency_months[mp]")
line(8, "ann_freq = annuity_frequency_months[mp]")
line(8, "age_idx = issue_index[mp]")
line(8, "sx = sex[mp]")
line(8, "cnt = count[mp]")
line(8, "prem = premium[mp]")
line(8, "annuity = annuity_payment[mp]")
line(8, "c_start = coverage_offset[mp]")
line(8, "c_end = coverage_offset[mp + 1]")
line(8, "ss = start_state[mp]")
line(8, "pf = 1.0")
line(8, "ann_prem = prem * 12.0 / prem_freq")
emit_init(8)
line(8, "pv_mortality = 0.0")
line(8, "pv_morbidity = 0.0")
line(8, "pv_disability = 0.0")
line(8, "pv_premium = 0.0")
line(8, "pv_expense = 0.0")
line(8, "pv_annuity = 0.0")
line(8, "pv_surrender = 0.0")
line(8, "cum_premium = 0.0")
line(8, "last_year = -1")
line(8, "claim_rate = 0.0")
line(8, "morb_rate = 0.0")
line(8, "prem_due = 0")
line(8, "ann_due = 0")
line(8, "prem_left = premium_term")
# Saved per-month total in-force (sum across all cohorts). The rule
# and diagnosis passes below reuse this trajectory instead of
# re-running the state machine -- a 5x-10x win for semi-Markov
# contracts that combine a deep cohort grid with rule-bearing or
# diagnosis coverages. Same trick the detailed _project_kernel_semi_markov
# has been using since the cohort-aware projection landed.
line(8, "inforce_traj = np.empty(term)")
# --- Main t loop ---------------------------------------------------
line(8, "for t in range(boundary):")
line(12, "year = t // 12")
line(12, "if year != last_year:")
line(16, "claim_rate = 0.0")
line(16, "morb_rate = 0.0")
line(16, "for k in range(c_start, c_end):")
line(20, "cov_idx = coverage_index[k]")
line(20, "if coverage_is_diagnosis[cov_idx]:")
line(24, "continue # diagnosis coverages run separately")
line(20, "if coverage_waiting[k] != 0 or coverage_reduction_end[k] != 0:")
line(24, "continue # rule-bearing coverages run separately")
line(20, "rate = coverage_rates[cov_idx, sx, age_idx, year] * coverage_amount[k]")
line(20, "if coverage_risk[cov_idx] == 0:")
line(24, "claim_rate += rate")
line(20, "else:")
line(24, "morb_rate += rate")
line(16, "last_year = year")
line(16, "pf = premium_factor[sx, age_idx, year]")
line(16, "ann_prem = prem * pf * 12.0 / prem_freq")
line(12, f"inforce_t = {sum_all}")
line(12, "inforce_traj[t] = inforce_t")
line(12, f"prem_occ = {sum_prem}")
line(12, f"benefit_occ = {sum_ben}")
line(12, "discount_factor_bom_t = discount_factor_bom[t]")
line(12, "discount_factor_mid_t = discount_factor_mid[t]")
line(12, "if prem_due == 0 and prem_left > 0:")
line(16, "level = prem_occ * prem * pf")
line(16, "prem_due = prem_freq - 1")
line(12, "else:")
line(16, "level = 0.0")
line(16, "prem_due -= 1")
line(12, "prem_left -= 1")
line(12, "pv_premium += level * discount_factor_bom_t")
if use_surrender and not surrender_is_amount:
line(12, "cum_premium += level")
line(12, "pv_mortality += inforce_t * claim_rate * discount_factor_mid_t")
line(12, "pv_morbidity += inforce_t * morb_rate * discount_factor_mid_t")
if use_annuity:
line(12, "if ann_due == 0:")
line(16, "pv_annuity += inforce_t * annuity * annuity_factor[sx, age_idx, year] * discount_factor_bom_t")
line(16, "ann_due = ann_freq - 1")
line(12, "else:")
line(16, "ann_due -= 1")
line(12, "pv_disability += benefit_occ * disability_income[mp] * discount_factor_mid_t")
if use_surrender:
line(12, f"lapse_flow = {sum_lapse}")
if surrender_is_amount:
# amount_per_policy / amount_per_unit: surrender_curve[t] is the
# surrender amount at duration t (per policy, or per unit of
# surrender_base[mp]); lapse_flow is the number lapsing.
# surrender_base is 1.0 for amount_per_policy.
line(12, "pv_surrender += (lapse_flow * surrender_curve[t]")
line(12, " * surrender_base[mp] * discount_factor_mid_t)")
else:
# cum_premium_factor: cum_premium aggregates inforce * premium; the
# effective lapse fraction is lapse_flow / inforce_t (the raw rate for a
# single state). inforce_t carries the count, so dividing it out avoids a
# cnt^2 scaling.
line(12, "eff_lapse = lapse_flow / inforce_t if inforce_t > 0.0 else 0.0")
line(12, "pv_surrender += (eff_lapse")
line(12, " * cum_premium * surrender_curve[t] * discount_factor_mid_t)")
line(12, "acquisition_expense = cnt * (acquisition_premium * ann_prem + acquisition_per_policy) if t == 0 else 0.0")
line(12, "maintenance_premium_expense = inforce_t * maintenance_premium[t] * ann_prem / 12.0 if t < premium_term else 0.0")
line(12, "maintenance_per_policy_expense = inforce_t * maintenance_per_policy[t]")
if use_lae:
line(12, "lae_expense = lae[t] * "
"inforce_t * (claim_rate + morb_rate)")
line(12, "pv_expense += (acquisition_expense + maintenance_premium_expense + maintenance_per_policy_expense + lae_expense) * discount_factor_mid_t")
else:
line(12, "pv_expense += (acquisition_expense + maintenance_premium_expense + maintenance_per_policy_expense) * discount_factor_mid_t")
emit_edge_step(12, include_lump=True)
line(8, f"total = {sum_all}")
line(8, "pm = (total * maturity_benefit[mp] * discount_factor_bom[boundary]) "
"if boundary == term else 0.0")
# --- Coverage-rule pass --------------------------------------------
# Rule-bearing non-diagnosis coverages: reuse the per-month total
# in-force trajectory the main pass saved instead of re-running the
# cohort-aware state machine. Each coverage becomes an O(term) scalar
# loop -- the same trick the detailed _project_kernel_semi_markov
# has been using all along.
line(8, "for k in range(c_start, c_end):")
line(12, "cov_idx = coverage_index[k]")
line(12, "if coverage_is_diagnosis[cov_idx]:")
line(16, "continue")
line(12, "wait = coverage_waiting[k]")
line(12, "red_end = coverage_reduction_end[k]")
line(12, "if wait == 0 and red_end == 0:")
line(16, "continue # rule-free -- already in the main pass")
line(12, "benefit = coverage_amount[k]")
line(12, "red_factor = coverage_reduction_factor[k]")
line(12, "mortality_risk = coverage_risk[cov_idx] == 0")
line(12, "for t in range(wait, boundary):")
line(16, "year = t // 12")
line(16, "mult = red_factor if t < red_end else 1.0")
line(16, "contrib = (inforce_traj[t] * coverage_rates[cov_idx, sx, age_idx, year]")
line(16, " * benefit * mult * discount_factor_mid[t])")
line(16, "if mortality_risk:")
line(20, "pv_mortality += contrib")
line(16, "else:")
line(20, "pv_morbidity += contrib")
# --- Diagnosis-coverage pass ---------------------------------------
# Diagnosis coverages run off a depleting not-yet-diagnosed pool.
# The pool's depletion (`undiagnosed`) is a scalar that multiplies the
# saved cohort-aware in-force trajectory -- mathematically equivalent
# to running the state machine with each flow scaled by (1 - d_rate),
# but a single scalar loop per coverage rather than a full cohort walk.
line(8, "for k in range(c_start, c_end):")
line(12, "cov_idx = coverage_index[k]")
line(12, "if not coverage_is_diagnosis[cov_idx]:")
line(16, "continue")
line(12, "benefit = coverage_amount[k]")
line(12, "wait = coverage_waiting[k]")
line(12, "red_end = coverage_reduction_end[k]")
line(12, "red_factor = coverage_reduction_factor[k]")
line(12, "undiagnosed = 1.0")
line(12, "d_year = -1")
line(12, "d_rate = 0.0")
line(12, "for t in range(boundary):")
line(16, "year = t // 12")
line(16, "if year != d_year:")
line(20, "d_rate = coverage_rates[cov_idx, sx, age_idx, year]")
line(20, "d_year = year")
line(16, "if t >= wait:")
line(20, "mult = red_factor if t < red_end else 1.0")
line(20, "pv_morbidity += (inforce_traj[t] * undiagnosed * d_rate * benefit")
line(20, " * mult * discount_factor_mid[t])")
line(16, "undiagnosed *= (1.0 - d_rate)")
# --- Final output --------------------------------------------------
line(8, "bel_mp = (pv_mortality + pv_morbidity + pv_disability + pm")
line(8, " + pv_annuity + pv_expense + pv_surrender - pv_premium)")
line(8, "ra_mp = (mortality_factor * pv_mortality + morbidity_factor * pv_morbidity")
line(8, " + disability_factor * pv_disability "
"+ longevity_factor * (pm + pv_annuity))")
line(8, "fcf = bel_mp + ra_mp")
line(8, "bel[mp] = bel_mp")
line(8, "ra[mp] = ra_mp")
line(8, "csm[mp] = max(0.0, -fcf)")
line(8, "loss_component[mp] = max(0.0, fcf)")
line(4, "return bel, ra, csm, loss_component")
return "\n".join(L)
_SEMI_MARKOV_KERNEL_CACHE: dict = {}
def _get_semi_markov_kernel(
n_states, state_duration_max, periodic_benefit_term_months,
edge_from, edge_to, edge_lump_sum,
state_pays_premium, state_pays_benefit,
use_annuity=True, use_lae=True, use_surrender=True,
surrender_is_amount=False,
):
"""Return a semi-Markov codegen-specialised kernel for the given
topology + per-state cohort counts.
Same two-level cache pattern as :func:`_get_markov_kernel` --
process-local dict for in-memory hits, content-addressed ``.py`` file
on disk so ``@njit(cache=True)`` persists the compiled native code
across processes.
"""
key = (
int(n_states),
tuple(int(x) for x in state_duration_max),
# ``periodic_benefit_term_months`` changes which cohorts ``sum_ben`` pays, so
# two models that differ only by the cap MUST NOT share a kernel --
# omitting it here would serve a stale kernel and silently mis-state
# the disability cash flow / BEL.
tuple(int(x) for x in periodic_benefit_term_months),
tuple(int(x) for x in edge_from),
tuple(int(x) for x in edge_to),
tuple(bool(x) for x in edge_lump_sum),
tuple(bool(x) for x in state_pays_premium),
tuple(bool(x) for x in state_pays_benefit),
bool(use_annuity), bool(use_lae), bool(use_surrender),
bool(surrender_is_amount),
)
cached = _SEMI_MARKOV_KERNEL_CACHE.get(key)
if cached is not None:
return cached
src = _semi_markov_kernel_source(
n_states, state_duration_max, periodic_benefit_term_months,
edge_from, edge_to, edge_lump_sum,
state_pays_premium, state_pays_benefit,
use_annuity=use_annuity, use_lae=use_lae, use_surrender=use_surrender,
surrender_is_amount=surrender_is_amount,
)
digest = hashlib.sha256(src.encode("utf-8")).hexdigest()[:16]
cache_path = _cache_dir() / f"fast_kernel_sm_{digest}.py"
if not cache_path.exists():
_atomic_write_text(cache_path, src)
module_name = f"_fastcashflow_codegen_sm_{digest}"
module = sys.modules.get(module_name)
if module is None:
spec = importlib.util.spec_from_file_location(module_name, cache_path)
module = importlib.util.module_from_spec(spec)
sys.modules[module_name] = module
spec.loader.exec_module(module)
kernel = module.kernel
_SEMI_MARKOV_KERNEL_CACHE[key] = kernel
return kernel
@njit(parallel=True, cache=True)
def _scalar_kernel(issue_index, sex, term_months, contract_boundary_months,
count, premium,
premium_term_months, premium_frequency_months,
annuity_frequency_months, coverage_index, coverage_amount, coverage_offset,
coverage_rates, premium_factor, annuity_factor, coverage_risk, coverage_is_diagnosis,
maturity_benefit, annuity_payment,
acquisition_premium, acquisition_per_policy, maintenance_premium,
maintenance_per_policy, lae,
discount_factor_bom, discount_factor_mid,
mortality_factor, morbidity_factor, longevity_factor,
coverage_waiting, coverage_reduction_end, coverage_reduction_factor,
coverage_pays_account_balance,
survival_monthly, lapse_monthly, surrender_curve,
use_morbidity, use_annuity, use_lae, use_surrender,
surrender_is_amount, surrender_base,
mortality_monthly,
has_account, mp_account, account_value0, account_face,
account_prem_to_av, account_coi_rate, account_admin_fee,
account_credit, account_charge, account_surr_charge,
account_gmab,
account_mortality_cv, account_expense_cv, account_z):
"""Scalar-inforce fast path of the general codegen fast kernel
(:func:`_markov_kernel_source`).
Used when the in-force projection collapses to a single survival track --
no user-supplied Model, no waiver inception, every model point
seated in the active state. The in-force is carried as a scalar; the
monthly decay is one multiply against the precomputed
``survival_monthly[sex, age, year] = (1 - q_monthly) * (1 - l_monthly)``
table. Numerically identical to the general codegen kernel for this configuration
-- the disability income, disability lump-sum and benefit-state pieces
of the general kernel evaluate to zero here -- and recovers the
pre-Phase(b) speed (see ``docs/tutorial/13-engine-design.md``).
"""
n_mp = issue_index.shape[0]
bel = np.empty(n_mp)
ra = np.empty(n_mp)
csm = np.empty(n_mp)
loss_component = np.empty(n_mp)
for mp in prange(n_mp):
term = term_months[mp]
boundary = contract_boundary_months[mp]
premium_term = premium_term_months[mp]
prem_freq = premium_frequency_months[mp]
ann_freq = annuity_frequency_months[mp]
age_idx = issue_index[mp]
sx = sex[mp]
cnt = count[mp]
prem = premium[mp]
annuity = annuity_payment[mp]
c_start = coverage_offset[mp]
c_end = coverage_offset[mp + 1]
inforce_t = cnt
pv_mortality = 0.0 # mortality-risk claim PV (death claims)
pv_morbidity = 0.0 # morbidity-risk claim PV (health claims)
pv_premium = 0.0
pv_expense = 0.0
pv_annuity = 0.0
pv_surrender = 0.0 # PV of surrender value
cum_premium = 0.0 # cumulative premium paid -- surrender basis
last_year = -1
claim_rate = 0.0
morb_rate = 0.0
# Universal-life account-value carrier (only rolled for account rows).
# The account is a single policy's balance, rolled with the within-month
# order verbatim from the standalone UL kernel; in-force weighting enters
# at the death / surrender / maturity / fund aggregation, not in the roll.
roll_av = has_account and mp_account[mp]
a_av = account_value0[mp] if roll_av else 0.0
face_av = account_face[mp] if roll_av else 0.0
cr_av = account_credit[mp] if roll_av else 0.0
half_credit_av = (1.0 + cr_av) ** 0.5
full_credit_av = 1.0 + cr_av
pv_account_death = 0.0 # PV of account death claim, inforce_t*q*max(av_mid,face)
pv_account_surr = 0.0 # PV of account surrender, inforce_t*l*(1-q)*av_mid
pv_account_nar = 0.0 # PV of net-amount-at-risk death (the UL RA base)
av_term = 0.0 # month-end account value at the contract boundary
# Counters replace modulo / less-than checks in the inner loop --
# ``prem_due`` ticks down to the next premium-paying month,
# ``ann_due`` to the next annuity month, and ``prem_left`` to the
# end of the premium-paying term. Profiling shows the modulo /
# comparison form costs ~2/3 of the inner-loop time at large
# portfolios -- the counter form lets the compiler keep the loop
# branch-light and 1M MP runs in ~50 ms again.
prem_due = 0
ann_due = 0
prem_left = premium_term
for t in range(boundary):
year = t // 12
if year != last_year:
claim_rate = 0.0
morb_rate = 0.0
for k in range(c_start, c_end):
cov_idx = coverage_index[k]
if coverage_is_diagnosis[cov_idx]:
continue
if coverage_pays_account_balance[cov_idx]:
continue # account-backed death pays from the AV
if coverage_waiting[k] != 0 or coverage_reduction_end[k] != 0:
continue
rate = coverage_rates[cov_idx, sx, age_idx, year] * coverage_amount[k]
if coverage_risk[cov_idx] == 0:
claim_rate += rate
else:
morb_rate += rate
last_year = year
discount_factor_bom_t = discount_factor_bom[t]
discount_factor_mid_t = discount_factor_mid[t]
if prem_due == 0 and prem_left > 0:
level = inforce_t * prem * premium_factor[sx, age_idx, year]
prem_due = prem_freq - 1
else:
level = 0.0
prem_due -= 1
prem_left -= 1
pv_premium += level * discount_factor_bom_t
pv_mortality += inforce_t * claim_rate * discount_factor_mid_t
if roll_av:
# Universal-life within-month account roll (the account-roll
# within-month order): premium in, COI on the net amount at risk,
# admin fee, every cost-deducting rider's fixed charge, floor at
# zero, half / full crediting. The account
# death / surrender / NAR claims settle on the half-credited
# av_mid, the maturity on the month-end balance.
a_av += account_prem_to_av[mp, t]
coi_nar = face_av - a_av
if coi_nar < 0.0:
coi_nar = 0.0
a_av -= (account_admin_fee[t] + account_coi_rate[mp, t] * coi_nar
+ account_charge[mp, t])
if a_av < 0.0:
a_av = 0.0
av_mid_t = a_av * half_credit_av
a_av = a_av * full_credit_av
av_term = a_av
# Account death pays max(av_mid, face) on the month's deaths
# (inforce_t * q); the pays_account_balance coverage was excluded
# from claim_rate above, so this is written ONCE here.
q_m = mortality_monthly[sx, age_idx, year]
deaths_m = inforce_t * q_m
best = av_mid_t if av_mid_t > face_av else face_av
pv_account_death += deaths_m * best * discount_factor_mid_t
# The net amount at risk (face above the account) is the only
# insurance-risk exposure -- the UL RA prices mortality on it.
nar = face_av - av_mid_t
if nar < 0.0:
nar = 0.0
pv_account_nar += deaths_m * nar * discount_factor_mid_t
# Surrender pays the account value on the mid-month lapse exits.
# The non-maturity, non-death exit count in the single-survival
# track is inforce_t * l * (1 - q) (the lapses net of the competing
# death decrement) -- matching the full path's exits - deaths.
l_m = lapse_monthly[sx, age_idx, year]
# Surrender pays the account net of the surrender charge.
surr_av = av_mid_t * (1.0 - account_surr_charge[mp, t])
if surr_av < 0.0:
surr_av = 0.0
pv_account_surr += inforce_t * l_m * (1.0 - q_m) * surr_av * discount_factor_mid_t
# Streams the portfolio does not use are skipped wholesale: the
# guards are loop-invariant per call, so the branch predicts
# perfectly and the per-cell array gathers / multiplies of an
# absent stream are not paid. A death-only book pays for none of
# surrender, annuity, morbidity or LAE.
if use_morbidity:
pv_morbidity += inforce_t * morb_rate * discount_factor_mid_t
if use_annuity:
if ann_due == 0:
pv_annuity += inforce_t * annuity * annuity_factor[sx, age_idx, year] * discount_factor_bom_t
ann_due = ann_freq - 1
else:
ann_due -= 1
ann_prem = prem * premium_factor[sx, age_idx, year] * 12.0 / prem_freq
acquisition_expense = (cnt * (acquisition_premium * ann_prem + acquisition_per_policy)
if t == 0 else 0.0)
maintenance_premium_expense = (inforce_t * maintenance_premium[t] * ann_prem / 12.0
if t < premium_term else 0.0)
maintenance_per_policy_expense = inforce_t * maintenance_per_policy[t]
if use_lae:
lae_expense = lae[t] * inforce_t * (claim_rate + morb_rate)
pv_expense += (acquisition_expense + maintenance_premium_expense + maintenance_per_policy_expense + lae_expense) * discount_factor_mid_t
else:
pv_expense += (acquisition_expense + maintenance_premium_expense + maintenance_per_policy_expense) * discount_factor_mid_t
if use_surrender:
if surrender_is_amount:
# amount_per_policy / amount_per_unit: surrender_curve[t] is
# the surrender amount at duration t (per policy, or per
# unit of surrender_base[mp]); inforce_t * lapse_rate is the
# number lapsing. surrender_base is 1.0 for amount_per_policy.
pv_surrender += (lapse_monthly[sx, age_idx, year]
* inforce_t * surrender_curve[t]
* surrender_base[mp] * discount_factor_mid_t)
else:
# cum_premium_factor: cum_premium aggregates inforce_t *
# premium and is the surrender basis; multiplying by
# lapse_rate alone gives the per-month surrender outflow
# (the count is already in cum_premium).
cum_premium += level
pv_surrender += (lapse_monthly[sx, age_idx, year]
* cum_premium * surrender_curve[t] * discount_factor_mid_t)
inforce_t *= survival_monthly[sx, age_idx, year]
if roll_av:
# Account maturity: survivors reaching the boundary take
# max(matured account value, GMAB) (the maturity benefit doubles as
# the guaranteed accumulation-benefit floor). Seeded at the boundary
# and discounted with the boundary start-of-month factor, mirroring
# the full path's _roll_forward_kernel maturity seed.
gmab = account_gmab[mp]
mat_av = av_term if av_term > gmab else gmab
pm = (inforce_t * mat_av * discount_factor_bom[boundary]
if boundary == term else 0.0)
else:
pm = (inforce_t * maturity_benefit[mp] * discount_factor_bom[boundary]
if boundary == term else 0.0)
# Non-diagnosis coverages with a waiting or reduced-benefit rule:
# rerun the survival on the same scalar track so the benefit
# multiplier (which can change mid-year) applies cleanly.
for k in range(c_start, c_end):
cov_idx = coverage_index[k]
if coverage_is_diagnosis[cov_idx]:
continue
if coverage_pays_account_balance[cov_idx]:
continue # account-backed death pays from the AV (above)
wait = coverage_waiting[k]
red_end = coverage_reduction_end[k]
if wait == 0 and red_end == 0:
continue
benefit = coverage_amount[k]
red_factor = coverage_reduction_factor[k]
mortality_risk = coverage_risk[cov_idx] == 0
inforce_t = cnt
for t in range(boundary):
year = t // 12
if t >= wait:
mult = red_factor if t < red_end else 1.0
contrib = (inforce_t * coverage_rates[cov_idx, sx, age_idx, year]
* benefit * mult * discount_factor_mid[t])
if mortality_risk:
pv_mortality += contrib
else:
pv_morbidity += contrib
inforce_t *= survival_monthly[sx, age_idx, year]
# Diagnosis coverages: claims run off a depleting "not yet diagnosed"
# pool, which depletes both by survival and by the diagnosis rate.
for k in range(c_start, c_end):
cov_idx = coverage_index[k]
if not coverage_is_diagnosis[cov_idx]:
continue
benefit = coverage_amount[k]
wait = coverage_waiting[k]
red_end = coverage_reduction_end[k]
red_factor = coverage_reduction_factor[k]
healthy = cnt
d_year = -1
d_rate = 0.0
for t in range(boundary):
year = t // 12
if year != d_year:
d_rate = coverage_rates[cov_idx, sx, age_idx, year]
d_year = year
if t >= wait:
mult = red_factor if t < red_end else 1.0
pv_morbidity += healthy * d_rate * benefit * mult * discount_factor_mid[t]
healthy *= survival_monthly[sx, age_idx, year] * (1.0 - d_rate)
if roll_av:
# Universal-life BEL: the slot mortality / surrender streams are
# empty for the account leg (the pays_account_balance coverage was
# excluded from claim_rate, and account surrender bypasses the slot
# surrender path), so the account death / surrender / maturity PVs
# enter here instead. The account value the entity holds (fund) is
# netted ONCE -- only the inception fund (inforce[0] * av0 = cnt*av0)
# survives into the as-of BEL, the premium that builds it counted
# once via pv_premium. The RA prices mortality on the net amount at
# risk plus expense risk, bypassing the slot RA factors.
bel_mp = (pv_account_death + pv_account_surr + pm
+ pv_morbidity + pv_annuity + pv_expense
- pv_premium - cnt * account_value0[mp])
# A cost-deducting rider's health benefit (pv_morbidity) bears
# morbidity risk -- priced alongside the at-risk mortality and
# expense (morbidity_factor == account_z * morbidity_cv). Zero for an
# account book with no such rider, so byte-identical.
ra_mp = (account_z * (account_mortality_cv * pv_account_nar
+ account_expense_cv * pv_expense)
+ morbidity_factor * pv_morbidity)
else:
bel_mp = pv_mortality + pv_morbidity + pm + pv_annuity + pv_expense + pv_surrender - pv_premium
ra_mp = (mortality_factor * pv_mortality + morbidity_factor * pv_morbidity
+ longevity_factor * (pm + pv_annuity))
fcf = bel_mp + ra_mp
bel[mp] = bel_mp
ra[mp] = ra_mp
csm[mp] = max(0.0, -fcf)
loss_component[mp] = max(0.0, fcf)
return bel, ra, csm, loss_component