Actuarial Validation & Filing Automation

Vectorized Sensitivity Analysis for Reserve Assumptions

Vectorized sensitivity analysis applies a whole matrix of assumption bumps to one base cash-flow projection and returns the reserve delta for every shock in a single NumPy pass, replacing the Python loop that revalues the block once per shock. This focused walkthrough builds that shock cube one array axis at a time. It is the compute kernel behind Sensitivity Analysis & Stress Testing and reuses the broadcasting habits taught in Optimizing Pandas DataFrames for Actuarial Cash Flow Projections.

The Problem

A sensitivity exhibit that shocks mortality, lapse, and interest each at several magnitudes needs the reserve recomputed dozens of times, and the obvious implementation — a for shock in shocks loop that calls the valuation engine and collects each result — pays the full projection cost once per iteration and spends most of its wall-clock time in Python-level dispatch rather than arithmetic. Because the base cash flows are identical across every shock and only a scalar multiplier or an additive rate shift changes, the entire grid is really a single tensor operation waiting to be written as one: stack the shocks on a new leading array axis, broadcast the base projection against them, and let NumPy compute every reserve at C speed in one pass. The task is to construct that broadcast correctly so the vectorized answer is bit-for-bit equal to the loop it replaces.

A Minimal Working Example

The example below revalues a small in-force block against a matrix of six shocks. The base projection is a decrement-and-discount present value; the shock matrix supplies, per shock, a mortality multiplier, a lapse multiplier, and an additive rate shift. Everything happens in one call, and the result is a reserve delta per shock. It runs as-is.

from __future__ import annotations

import numpy as np

rng = np.random.default_rng(20260717)

# --- base in-force block: 5 policies over 30 projection years ---------------
n_pol, n_yr = 5, 30
face_amount = np.array([250_000, 500_000, 100_000, 750_000, 300_000], dtype=np.float64)
base_qx = np.full((n_pol, n_yr), 0.008)          # base mortality_table rates
base_wx = np.full((n_pol, n_yr), 0.06)           # base lapse_rate vector
base_spot = np.linspace(0.045, 0.055, n_yr)      # base discount_curve spot rates
benefit = face_amount[:, None] * np.ones((n_pol, n_yr))

# --- shock matrix: (n_shocks, 3) = [mort_mult, lapse_mult, rate_shift] -------
shock_names = ["base", "mort+15%", "mort-15%", "lapse+25%", "rates-100bp", "rates+100bp"]
shocks = np.array(
    [
        [1.00, 1.00,  0.000],   # base
        [1.15, 1.00,  0.000],   # mortality up 15%
        [0.85, 1.00,  0.000],   # mortality down 15%
        [1.00, 1.25,  0.000],   # lapse up 25%
        [1.00, 1.00, -0.010],   # rates down 100bp
        [1.00, 1.00,  0.010],   # rates up 100bp
    ],
    dtype=np.float64,
)


def revalue_all_shocks(base_qx, base_wx, base_spot, benefit, shocks):
    """Return the reserve for every shock in one vectorized pass. Shape (n_shocks,)."""
    m = shocks[:, 0][:, None, None]                       # (S,1,1) mortality mult
    l = shocks[:, 1][:, None, None]                       # (S,1,1) lapse mult
    d = shocks[:, 2][:, None]                             # (S,1)   additive shift

    qx = np.clip(base_qx[None] * m, 0.0, 1.0)             # (S,P,Y)
    wx = np.clip(base_wx[None] * l, 0.0, 1.0)             # (S,P,Y)

    surv = np.cumprod(1.0 - qx - wx, axis=2)              # survival through year t
    bop = np.concatenate([np.ones_like(surv[:, :, :1]), surv[:, :, :-1]], axis=2)

    yr = np.arange(n_yr)
    disc = (1.0 + (base_spot[None] + d)) ** -(yr + 1)     # (S,Y) discount factors
    death_pv = bop * qx * benefit[None] * disc[:, None, :]
    return death_pv.sum(axis=(1, 2))                      # (S,) reserve per shock


reserves = revalue_all_shocks(base_qx, base_wx, base_spot, benefit, shocks)
base = reserves[shock_names.index("base")]
deltas = reserves - base

for name, r, delta in zip(shock_names, reserves, deltas):
    print(f"{name:>12}  reserve={r:14,.2f}  delta={delta:+13,.2f}")

The present value each shocked reserve computes is the discounted expected death benefit

Rs  =  ∑p∑t=1T  t−1px(s)⋅qx+t−1(s)⋅Bp,t⋅(1+it+δs)−t,R_s \;=\; \sum_{p}\sum_{t=1}^{T}\; {}_{t-1}p^{(s)}_{x} \cdot q^{(s)}_{x+t-1}\cdot B_{p,t}\cdot \bigl(1 + i_t + \delta_s\bigr)^{-t},

where the shock index ss scales the decrement rates and shifts the spot curve by δs\delta_s, and the reserve delta reported for the exhibit is simply Rs−RbaseR_s - R_{\text{base}}.

Vectorized reserve sensitivity: base projection broadcast against a shock matrix into a reserve delta vector A base cash-flow projection box and a shock matrix box of shape shocks by assumptions both feed a highlighted vectorized revaluation box that runs in one pass with no loop, which in turn produces a reserve delta vector of one delta per shock. Base cash-flow projection (P × Y) Shock matrix shape (S × A) Vectorized revalue one pass, no loop Reserve Δ vector one delta per shock (S,)
The base projection and the shock matrix broadcast into one vectorized revaluation, yielding a reserve delta for every shock.

How It Works, Block by Block

The shock matrix is the whole interface. Every shock is one row of three numbers — a mortality multiplier, a lapse multiplier, and an additive rate shift — so the grid of scenarios is a plain (n_shocks, 3) array. Keeping the shock specification as data means the base row lives alongside the perturbations and the delta is always computed against a reserve produced by the identical arithmetic, never a separately cached figure.

Reshaping is what creates the broadcast. The three columns are lifted onto new axes — the decrement multipliers become (S, 1, 1) so they broadcast across the policy and year axes of the base rates, and the rate shift becomes (S, 1) so it broadcasts across the year axis of the spot curve. Indexing base_qx[None] promotes the (P, Y) base to (1, P, Y), and multiplying by an (S, 1, 1) multiplier yields the full (S, P, Y) shocked-rate cube without a single Python-level iteration over shocks.

The survival and discount factors follow the same shape logic. np.cumprod along the year axis turns the per-year survival probabilities into cumulative in-force factors for every shock at once, and the concatenate prepends a column of ones so year one starts fully in force. The discount factors are computed per shock because the additive shift d moves the spot curve, giving a (S, Y) discount grid that broadcasts against the (S, P, Y) expected benefits.

The reduction collapses the cube to a vector. Summing the discounted expected death benefits over the policy and year axes leaves one reserve per shock, an (S,) vector. Subtracting the base entry produces the delta vector that becomes the sensitivity exhibit — the entire grid delivered by one revalue_all_shocks call rather than a loop.

Edge Cases and Production Hardening

Non-linear interactions are invisible to single-axis shocks. Shocking mortality and lapse in separate rows and adding their deltas is not the same as shocking both together when a product has a lapse-supported guarantee — the second-order interaction term is lost. When that interaction is material, add an explicit joint row to the matrix rather than inferring the combined effect by superposition:

# Joint stress captures the interaction a sum of single-axis deltas misses.
joint = np.array([[1.15, 1.25, -0.010]])          # mort+15%, lapse+25%, rates-100bp
shocks_ext = np.vstack([shocks, joint])
reserves_ext = revalue_all_shocks(base_qx, base_wx, base_spot, benefit, shocks_ext)
interaction = reserves_ext[-1] - reserves[1] - reserves[3] - reserves[4] + 2 * base

Large face amounts can overflow a narrow dtype. If the base rates or benefits are stored in float32 to save memory, a multi-million-dollar face amount accumulated across thirty discounted years can lose precision or, in an integer intermediate, overflow outright. Keep the projection in float64 and validate the input dtypes before the broadcast; the failure mode and its detection are covered in Preventing NumPy dtype Overflow in Reserve Projections.

The full shock cube can exhaust memory. The intermediate (S, P, Y) arrays scale with the product of shocks, policies, and years; a hundred shocks against a million policies over forty years is billions of float64 cells. Because the reserve sums over policies, stream the block in row chunks and accumulate:

def revalue_streaming(qx, wx, spot, benefit, shocks, chunk=100_000):
    total = np.zeros(shocks.shape[0])
    for s in range(0, qx.shape[0], chunk):
        sl = slice(s, s + chunk)
        total += revalue_all_shocks(qx[sl], wx[sl], spot, benefit[sl], shocks)
    return total

This caps peak memory at n_shocks × chunk × n_years × 8 bytes while producing a bit-identical total, because summation over policies is associative.

Compliance Note

This kernel is the computational evidence behind the sensitivity testing that ASOP No. 56 (Modeling) expects an actuary to perform and document — it calls for understanding how a model responds to changes in its significant assumptions. Under VM-20 the same deltas support the margin the appointed actuary holds for each material assumption and feed the sensitivity disclosures of the VM-31 PBR Actuarial Report. Because the shock matrix is data, the exact grid that produced a filed exhibit can be serialized and replayed, which is precisely the reproducibility a reviewer needs to reconstruct the numbers. For the surrounding process — cataloging shocks, streaming the block, and recording the run — see the parent guide on Sensitivity Analysis & Stress Testing.

Frequently Asked Questions

Why is the vectorized version faster than a loop over shocks?

The loop pays Python-level dispatch and re-traverses the base projection once per shock, spending most of its time outside the arithmetic. The vectorized version stacks the shocks on a leading array axis and lets NumPy compute every reserve in one C-level pass, so the cost of a dozen shocks approaches the cost of one full projection rather than a dozen.

Does broadcasting change the reserve compared with revaluing one shock at a time?

No. The broadcast performs the identical arithmetic the loop would, just organized as one tensor operation, so the vectorized reserves are bit-for-bit equal to a per-shock loop up to floating-point associativity. Asserting that equality in a test is worthwhile because any divergence points to a reshape bug rather than a modeling difference.

How do I capture the interaction between two assumptions?

Add an explicit joint row to the shock matrix that bumps both axes at once, then compare its delta to the sum of the two single-axis deltas. The difference is the interaction term, which matters for products like lapse-supported guarantees where mortality and lapse effects are not additive.

What keeps the shock cube from exhausting memory on a large block?

Stream the in-force block in row chunks and accumulate the per-shock reserve across chunks. Because the reserve sums over policies, chunked accumulation is exact, and peak memory stays proportional to the chunk size times the number of shocks rather than the whole portfolio.

Up a level: Sensitivity Analysis & Stress Testing · Actuarial Model Ingestion & Testing Workflows

Regulatory references: Actuarial Standards Board ASOP No. 56 (Modeling); NAIC Valuation Manual VM-20 and VM-31.