Mathematical Demography, Multiple Decrement Models & Population Dynamics

In the single-decrement survival models analyzed in Unit 3, death was the sole terminating event. In real-world actuarial environments—such as defined-benefit corporate pension funds, health insurance schemes, and life insurance portfolios—an individual may exit the active insured population due to one of several mutually exclusive competing causes:

• Cause 1: Normal Age Retirement
• Cause 2: Ill-Health Total Disability
• Cause 3: Death in Active Service
• Cause 4: Voluntary Withdrawal (Resignation / Policy Surrender)

The mathematical theory governing these competing pathways is Multiple Decrement Theory.

The Continuous Multiple Decrement Model

Let $(x)$ be an active life aged $x$. The time until exit is modeled by a continuous random variable $T \ge 0$, and the cause of exit is modeled by a discrete random variable $J \in \{1, 2, \dots, m\}$.

The Cause-Specific Force of Decrement for cause $j$ at attained age $x + t$ is:

The Total Force of Decrement across all operating causes is the direct sum of cause-specific hazards:

The probability that a life aged $x$ survives all $m$ decrements for at least $t$ years is:

Crude Decrement Probabilities: ${}_t q_x^{(j)}$

The marginal probability that $(x)$ leaves the active cohort before duration $t$ specifically due to cause $j$ in the presence of all other competing decrements is the Crude Decrement Probability:

Summing across all $m$ causes reproduces the total decrement probability:

When $t = 1$, the annual crude decrement rate is written simply as $q_x^{(j)}$.

Actuaries frequently confront counterfactual policy questions: What would the mortality rate of factory workers be if industrial toxic exposures were completely eradicated? Or: What would the pension fund retirement rate be if early disability exits were eliminated?

Answering these questions requires isolating cause $j$ from all other competing risks.

The Associated Single Decrement Table

For each cause $j$, actuaries construct a theoretical Associated Single Decrement Table (ASDT), representing an imaginary universe where cause $j$ is the only operating hazard.

The Net Decrement Probability $q_x^{\prime (j)}$ represents the probability that a life aged $x$ exits within one year due to cause $j$ assuming all other decrements are inoperative:

Because an individual cannot be saved from death by withdrawing from a pension scheme, net probabilities satisfy:

The crude probability $q_x^{(j)}$ is strictly smaller than the net probability $q_x^{\prime (j)}$ because competing causes remove lives who would otherwise have exited via cause $j$ later in the year!

Mathematical Conversion Relationships

To translate between crude rates $q_x^{(j)}$ and net rates $q_x^{\prime (j)}$, actuaries apply one of two fractional assumptions:

1. Constant Force of Decrement (CFD)

Assuming the ratio of each cause-specific hazard to total hazard remains constant over the year:

Integrating yields the celebrated exponential conversion:

2. Uniform Distribution of Decrements (UDD in Multiple Decrement Table)

Assuming each crude decrement distributes uniformly over the year (${}_t q_x^{(j)} = t \cdot q_x^{(j)}$), for a two-decrement model ($m = 2$):

The denominator $1 - \frac{1}{2} q_x^{(2)}$ adjusts for the fact that lives exiting via cause 2 on average exit at mid-year ($t = 0.5$), exposing them to cause 1 for only half the year.

Many actuarial contracts involve multiple lives simultaneously (e.g., joint life annuities for married couples, or survivor pensions for widows).

Let $(x)$ and $(y)$ be two lives with remaining future lifetimes $T_x$ and $T_y$.

1. The Joint Life Status: $(xy)$

The Joint Life Status (denoted $(xy)$) is a status that functions as long as both lives survive, failing upon the first death:

Assuming $T_x$ and $T_y$ are statistically independent:

The force of mortality of the joint status is additive:

The joint life annuity-due is:

2. The Last Survivor Status: $(\overline{xy})$

The Last Survivor Status (denoted $(\overline{xy})$) functions as long as at least one life survives, failing only upon the second death:

The probability of surviving past $t$ is:

Duality Identities for Multi-Life Contracts

Taking discounted integrals and sums yields the universal multi-life addition theorems:

These identities permit the immediate pricing of complex joint-and-survivor annuity options from single-life and joint-life tables.

In 1907–1939, Alfred J. Lotka revolutionized mathematical demography by formulating the continuous renewal theory of human population dynamics.

Let $B(t)$ denote the continuous birth rate (births per unit time) in a closed female population at calendar time $t$. Let $p(a) = {}_a p_0$ denote the probability that a female newborn survives to exact age $a$. Let $m(a)$ denote the age-specific fertility rate (the annual rate of female childbearing for women aged $a$, restricted to reproductive ages $[\alpha, \beta] \approx [15, 49]$).

The Lotka Renewal Integral Equation

Females born at time $t - a$ who survive to time $t$ are aged $a$. Their number is $B(t - a) p(a)$. Each woman of age $a$ gives birth at rate $m(a)$.

Summing births across all maternal ages:

where $G(t)$ represents births contributed by females who were already alive at initial time $t = 0$. For $t > \beta \approx 50\text{ years}$, the initial population has aged past reproduction ($G(t) = 0$), and the renewal equation becomes purely endogenous:

The Characteristic Equation for the Intrinsic Growth Rate ($r$)

Lotka sought trial exponential solutions of the form $B(t) = C e^{r t}$. Substituting into the renewal equation:

Dividing both sides by $C e^{r t}$ yields Lotka's Characteristic Equation:

The unique real root $r$ of this equation is the Intrinsic Rate of Natural Increase (Lotka's $r$).

• If $\int_\alpha^\beta p(a) m(a)\,da > 1.0$: $r > 0$ (the population expands exponentially).
• If $\int_\alpha^\beta p(a) m(a)\,da = 1.0$: $r = 0$ (the population reaches stationary equilibrium).
• If $\int_\alpha^\beta p(a) m(a)\,da < 1.0$: $r < 0$ (the population declines exponentially).

Net and Gross Reproduction Rates

1. Gross Reproduction Rate ($GRR$): The average number of female daughters a woman would bear if she survived to the end of her reproductive years ($\beta$):

2. Net Reproduction Rate ($NRR$, $R_0$): The average number of female daughters a newborn girl will bear, accounting for cohort mortality:

Notice that $R_0$ is the characteristic equation evaluated at $r = 0$.

The Asymptotic Stable Age Distribution

Lotka proved that regardless of the initial age distribution of a population, if age-specific fertility $m(a)$ and survival $p(a)$ remain constant over time, the population will inexorably converge to a fixed, time-invariant Stable Age Distribution $c(a)$:

where $b = \left[ \int_0^\omega e^{-r a} p(a)\,da \right]^{-1}$ is the stable crude birth rate. This result—the Strong Ergodic Theorem of Demography—establishes that a population completely forgets its historical demographic shocks!

In 1945, P. H. Leslie discretized Lotka’s continuous renewal theory into an elegant matrix algebra formulation widely applied in official national demographic forecasting.

Divide the female population into $k$ equal age groups of width $h$ years (typically $h = 5\text{ years}$, giving age brackets $0-4, 5-9, \dots, 80-84, 85+$). Let the state vector at time step $t$ be:

where $n_{i, t}$ is the female population count in age group $i$ at time $t$.

The Leslie Projection Matrix ($\mathbf{L}$)

The population evolves from time step $t$ to $t + 1$ (after $h$ years) via the linear matrix system:

The Leslie Matrix $\mathbf{L}$ is a $k \times k$ square matrix with a canonical sparse structure:

where:

• Top Row ($F_i$): The net effective reproduction coefficient for age group $i$, measuring daughters born during the 5-year interval who survive to enter age group 1 ($0-4$):

• Sub-Diagonal ($P_i$): The cohort survival probability from age group $i$ to age group $i + 1$:

• All other matrix entries are zero.

Spectral Analysis and Dominant Eigenvalue

Because $\mathbf{L}$ is a non-negative, irreducible matrix, the Perron-Frobenius Theorem guarantees:

1. $\mathbf{L}$ has a unique, positive real eigenvalue $\lambda_1 > 0$ that strictly exceeds the modulus of all other complex eigenvalues ($|\lambda_j| < \lambda_1$ for $j \ge 2$).
2. The long-term asymptotic growth rate is governed entirely by $\lambda_1$:

where $\lambda_1 = e^{r h}$ is directly related to Lotka's intrinsic rate of natural increase $r$.

1. The right eigenvector $\mathbf{v}_1$ corresponding to $\lambda_1$ defines the asymptotic stable age structure.
2. The left eigenvector $\mathbf{u}_1^T$ defines Fisher's Reproductive Value across age brackets.

In Section 8.1, the aggregate claims random variable was defined as the random sum:

where $N$ is the annual claim frequency random variable and $\{X_i\}_{i=1}^\infty$ are independent, identically distributed positive claim severity random variables with probability distribution $f_X(x) = P(X = x)$ for discrete claim sizes $x \in \{1, 2, \dots, m\}$.

Calculating the probability mass function of aggregate claims $g_k = P(S = k)$ via direct multi-fold convolutions requires $\mathcal{O}(k^2)$ operations per fold, becoming computationally intractable when $N$ can assume hundreds or thousands of claims. In 1981, Harry Panjer introduced a revolutionary recursive formula that computes the exact distribution of $S$ in linear time $\mathcal{O}(k)$.

1. The Panjer $(a, b, 0)$ Class of Frequency Distributions

Panjer identified that the primary counting distributions in actuarial risk theory satisfy a common linear recursive relation:

The distributions satisfying this relation form the Panjer $(a, b, 0)$ Class:

1. Poisson($\lambda$): $a = 0, \; b = \lambda, \; p_0 = e^{-\lambda}$.

2. Binomial($m, p$): $a = -\frac{p}{1-p}, \; b = \frac{(m+1)p}{1-p}, \; p_0 = (1-p)^m$.

3. Negative Binomial($r, p$): $a = 1 - p, \; b = (r-1)(1-p), \; p_0 = p^r$.

4. Geometric($p$): Special case of Negative Binomial with $r = 1$: $a = 1 - p, \; b = 0$.

2. Statement and Proof of the Panjer Recursion

Let $N$ belong to the $(a, b, 0)$ class, and let severity $X$ take values on the positive integers $\{1, 2, \dots\}$ with $f_j = P(X = j)$. The aggregate claims distribution $g_k = P(S = k)$ satisfies:

When severity is strictly positive ($f_0 = P(X = 0) = 0$), the leading factor $1 - a f_0 = 1$.

Mathematical Derivation: The probability generating function of $S$ is $G_S(z) = G_N(G_X(z))$. From the definition of the $(a, b, 0)$ class:

Differentiating $G_S(z) = G_N(G_X(z))$ with respect to $z$:

Substituting $G_N'(G_X(z)) = \frac{a + b}{1 - a G_X(z)} G_N(G_X(z))$ and equating coefficient series powers of $z^{k-1}$ produces the exact Panjer recursive identity.

Because each term $g_k$ depends only on previously computed values $\{g_0, g_1, \dots, g_{k-1}\}$, the entire probability distribution of aggregate claims can be generated in a single sequential sweep, enabling real-time pricing of aggregate stop-loss reinsurance contracts.

In classical continuous-time ruin theory, the insurer's surplus process $\{U(t)\}_{t \ge 0}$ evolves according to:

where:

• $u = U(0) \ge 0$ is the initial surplus capital.
• $c > 0$ is the continuous premium collection rate per unit time.
• $N(t)$ is a homogeneous Poisson process with claim arrival intensity $\lambda > 0$.
• $X_i$ are independent, identically distributed positive claim amounts with cumulative distribution function $F_X(x)$ and mean $\mu_X = \mathbb{E}[X] < \infty$.

The continuous premium collection rate includes a safety loading $\theta > 0$:

If $\theta \le 0$, ruin occurs with probability $1.0$ as $t \to \infty$ by the Strong Law of Large Numbers. Thus, the positive safety loading condition $\theta > 0$ is strictly required for long-term survival.

1. The Probability of Ruin

Let $T = \inf\{t > 0 : U(t) < 0\}$ denote the time of ruin (with $T = \infty$ if $U(t) \ge 0$ for all $t \ge 0$). The ultimate ruin probability regarded as a function of initial capital $u$ is:

The non-ruin (survival) probability is $\phi(u) = 1 - \psi(u)$.

2. Derivation of the Integro-Differential Equation

Consider a minuscule elapsed time interval $h > 0$. During $[0, h]$:

• With probability $1 - \lambda h + o(h)$, no claims occur, and the surplus grows deterministically to $u + c h$.
• With probability $\lambda h + o(h)$, exactly one claim of size $X$ occurs.
• The probability of two or more claims is $o(h)$.

By conditioning on the events occurring during $[0, h]$:

Rearranging terms, dividing by $c h$, and taking the limit as $h \to 0^+$:

This is the celebrated Integro-Differential Equation of Ruin Theory.

3. Exact Solution for Exponential Claims

When individual claim severities are exponentially distributed with parameter $\alpha > 0$ ($f_X(x) = \alpha e^{-\alpha x}$ for $x > 0$), so that $\mu_X = 1/\alpha$, the integro-differential equation can be solved in closed analytical form.

Differentiating the equation and substituting the exponential density reduces the integro-differential equation to a second-order linear ordinary differential equation with constant coefficients:

Using the boundary conditions:

1. $\lim_{u \to \infty} \psi(u) = 0$ (infinite capital eliminates ruin).
2. Initial surplus condition at $u = 0$: $\psi(0) = \frac{\lambda}{c \alpha} = \frac{\lambda \mu_X}{(1 + \theta) \lambda \mu_X} = \frac{1}{1 + \theta}$.

The exact closed-form solution for the ultimate ruin probability is:

This analytical result demonstrates that:

• Even with zero initial capital ($u = 0$), the insurer has a survival probability of $\phi(0) = 1 - \frac{1}{1 + \theta} = \frac{\theta}{1 + \theta} > 0$.
• The probability of ruin decays exponentially fast with respect to initial capital $u$, where the exponential decay parameter $R = \frac{\theta \alpha}{1 + \theta}$ is the exact Lundberg adjustment coefficient.

Defined benefit (DB) pension plans promise retirement annuities defined by a benefit formula: $$\text{Annual Pension} = k \cdot n \cdot S_{\text{final}}$$ where $k$ is the annual accrual rate (e.g., $1.5\% - 2.0\%$), $n$ is years of credited service, and $S_{\text{final}}$ is the final average salary (e.g., average of the final three years). To determine annual employer contributions, consulting actuaries utilize cost allocation methods governed by the International Accounting Standard (IAS 19) and the US Financial Accounting Standards Board (ASC 715). ### 1. The Projected Unit Credit (PUC) Method The **Projected Unit Credit** method allocates total projected retirement benefits across years of service accrued to date: - At attained age $x$ with entry age $e$, elapsed service is $y = x - e$, and total service at normal retirement age $r$ is $n = r - e$. - The **Actuarial Accrued Liability (AAL)** is the present value of the projected retirement benefit earned for service rendered prior to the valuation date: $$\text{AAL}_x^{\text{PUC}} = \left( \frac{x - e}{r - e} \right) \cdot PVFB_x$$ where $PVFB_x$ is the Present Value of Future Benefits at age $x$, incorporating salary growth projections $s_r / s_x$. - The **Normal Cost (NC)** is the present value of the additional benefit unit allocated to the current year: $$\text{NC}_x^{\text{PUC}} = \frac{1}{r - e} \cdot PVFB_x = \frac{\text{AAL}_x^{\text{PUC}}}{x - e}$$ *Characteristic Behavior*: Because $PVFB_x$ increases exponentially with age as the employee nears retirement (due to fewer years of interest discounting and higher survival probabilities), the PUC normal cost per employee rises rapidly with age. ### 2. The Entry Age Normal (EAN) Method The **Entry Age Normal** method treats the pension funding cost as a level percentage of employee salary $C$ payable from initial employment age $e$ to retirement age $r$: - The contribution rate $C$ is determined at entry age $e$ to equate the present value of future normal costs to the present value of future benefits: $$C \cdot \ddot{a}_{e:\overline{r-e}|}^{s} = PVFB_e \implies C = \frac{PVFB_e}{\ddot{a}_{e:\overline{r-e}|}^{s}}$$ where $\ddot{a}_{e:\overline{r-e}|}^{s}$ is a salary-weighted annuity factor. - At attained age $x$, the **Normal Cost** is simply: $$\text{NC}_x^{\text{EAN}} = C \cdot s_x$$ - The **Actuarial Accrued Liability** under EAN is computed prospectively: $$\text{AAL}_x^{\text{EAN}} = PVFB_x - C \cdot s_x \cdot \ddot{a}_{x:\overline{r-x}|}^{s}$$ or retrospectively: $$\text{AAL}_x^{\text{EAN}} = C \sum_{j=e}^{x-1} s_j \frac{{}_j p_e (1 + i)^{x-j}}{{}_{x-e} p_e}$$ The difference between the Actuarial Accrued Liability and the market value of plan assets ($AAL - \text{Assets}$) represents the **Unfunded Actuarial Liability (UAL)**, which employers amortize over statutory periods (typically 15 to 30 years).

While Panjer's recursion provides an exact method for computing aggregate claims distributions when claim frequency belongs to the $(a, b, 0)$ or $(a, b, 1)$ classes, many empirical insurance portfolios feature claim frequency distributions outside these families (e.g., Delaporte distributions, zero-modified Poisson-inverse Gaussian mixtures).

To compute the aggregate claims distribution $g_k = P(S = k)$ for arbitrary frequency and continuous severity distributions, modern actuarial science utilizes Characteristic Function Inversion via the Fast Fourier Transform (FFT).

1. The Characteristic Function of Aggregate Claims

Recall that the aggregate claims random variable is $S = \sum_{j=1}^N X_j$. Let $\phi_X(t) = \mathbb{E}[e^{itX}]$ denote the characteristic function of individual claim severity $X$, and let $P_N(z) = \mathbb{E}[z^N]$ denote the probability generating function of claim count $N$.

By the tower property of conditional expectation:

For a Compound Poisson process with frequency parameter $\lambda$:

For a Compound Negative Binomial process with parameters $(r, p)$:

2. Discretization and Discrete Fourier Transform (DFT)

Let claim severity be discretized on a regular grid of size $h > 0$ with $M = 2^m$ points ($M \in \{1024, 2048, 4096, 8192\}$):

The Discrete Fourier Transform of the severity vector $\mathbf{f}$ is:

This is computed in $\mathcal{O}(M \log M)$ operations using the Cooley-Tukey FFT algorithm.

Next, the aggregate characteristic values $\hat{g}_k$ are evaluated directly through the theoretical generating relation:

Finally, the probability mass function of aggregate claims $\mathbf{g} = [g_0, g_1, \dots, g_{M-1}]^T$ is recovered via the Inverse Fast Fourier Transform (IFFT):

3. Mitigating Aliasing and Wrapping Errors: Exponential Tilting

Because the discrete Fourier transform implicitly assumes periodic boundary conditions with period $M$, aggregate claim probabilities beyond $M$ wrap around and contaminate early probabilities $g_0, g_1$, creating Aliasing / Wrapping Error.

To eliminate aliasing, actuaries deploy Exponential Tilting (Esscher Transform):

1. Define the tilted severity distribution:

where $\theta > 0$ is a damping parameter.

1. Run the FFT and IFFT on the tilted distribution to obtain the tilted aggregate probabilities $\tilde{g}_j$.
2. Untilt the resulting aggregate distribution:

Exponential damping suppresses tail values prior to periodicity wrapping, reducing aliasing error below machine precision ($10^{-14}$) and ensuring pristine aggregate stop-loss reinsurance pricing.

In defined benefit pension management, traditional asset allocation models that focus purely on maximizing asset returns (e.g., $60\%$ equities / $40\%$ bonds) expose pension plans to catastrophic funding ratio collapses when interest rates drop.

To guarantee solvency, pension funds deploy Liability-Driven Investment (LDI).

1. The Pension Funding Ratio and Surplus Volatility

Let $A$ be the market value of pension plan assets, and let $L$ be the Present Value of Defined Benefit Liabilities ($PVFB$ or $AAL$). The Funding Ratio is:

The pension Surplus is:

The variance of changes in pension surplus $\Delta S = \Delta A - \Delta L$ is:

where $\sigma_A$ is asset volatility, $\sigma_L$ is liability volatility, and $\rho_{AL}$ is the asset-liability correlation.

If assets are invested in equities while liabilities consist of long-duration guaranteed bonds, the asset-liability correlation $\rho_{AL}$ is close to zero or negative. Consequently, surplus volatility is massive ($\sigma_S^2 \approx \sigma_A^2 + \sigma_L^2$). If interest rates drop by 100 basis points, liabilities surge by $20\%$, causing the funding ratio to collapse into deep deficit.

2. The LDI Portfolio Bifurcation

LDI decomposes the pension portfolio into two dedicated structural allocations:

1. The Liability Hedging Portfolio (LHP): Designed strictly to achieve $\rho_{AL} \approx 1.0$ and match the key rate duration and inflation sensitivity of liabilities:

• Physical long-dated sovereign bonds and index-linked gilts.
• Derivative overlays: Interest Rate Swaps, Inflation Swaps, and Swaptions. In a receiver swap, the pension fund receives fixed interest and pays floating, synthesizing long bond duration with minimal cash capital outlay.

2. The Return-Seeking Portfolio (RSP): Invested in high-expected-return growth assets (global equities, private infrastructure, private credit, and real estate) to systematically eliminate any existing unfunded actuarial liability ($UAL$).

By matching liability duration via synthetic swaps in the LHP, $\rho_{AL} \to 1.0$ and $\sigma_L \approx \sigma_A^{\text{hedging}}$, surplus volatility $\sigma_S^2 \to 0$. This stabilizes employer pension contribution rates, insulates corporate sponsors from balance sheet volatility, and secures the retirement pensions of millions of workers.