Numerical Solutions for Time-Inhomogeneous Kolmogorov Equations

Introduction

In our last lesson, we established the fundamental differential equations—Kolmogorov's forward and backward equations—that govern time-inhomogeneous Markov processes. We saw that while these equations perfectly describe the system's evolution, their time-dependent nature means they generally cannot be solved with a simple matrix exponential. This leaves us with a critical question: how do we actually find the transition probabilities needed for actuarial calculations?

This lesson will answer that question directly. We will learn a practical, powerful technique to approximate the solutions to these differential equations: Euler's method. You will learn the intuition behind this numerical approach and, more importantly, how to apply it step-by-step to solve for transition probabilities in the kind of multi-state models you will encounter in the CS2 exam. Your preference for calculation and problem-solving will be well-served today as we move from theoretical equations to concrete numerical results.

The Essence of Euler's Method

At its core, solving a differential equation like $\frac{d y}{d t} = f (t, y)$ means finding a function $y (t)$ whose slope at any point $(t, y)$ is given by $f (t, y)$ . If we can't find a formula for $y (t)$ , we can try to build the function piece by piece.

This is exactly what Euler's method does. Starting from a known point $(t_{0}, y_{0})$ , we calculate the slope $f (t_{0}, y_{0})$ . We then assume the function is a straight line with that slope for a very small step forward in time, $h$ . This allows us to predict the next point, $(t_{1}, y_{1})$ , and repeat the process.

The video below provides an excellent graphical illustration of this idea.

Euler method | Lecture 48 | Numerical Methods for Engineers

This video from Jeffrey Chasnov's "Numerical Methods for Engineers" series introduces the concept of Euler's method for a first-order ordinary differential equation (ODE).

Please watch from the beginning to 05:29. The key aspects to focus on are: The initial setup of the problem: an ODE with an initial condition. The graphical interpretation of "marching forward" along the tangent line. The derivation of the iterative formula: x_{n+1} = x_n + \Delta t \cdot f(t_n, x_n).

As the video explains, this process gives us a simple, iterative formula. If we know the value of our function at time $t_{n}$ , we can approximate its value at time $t_{n + 1} = t_{n} + h$ as:

y (t_{n} + h) \approx y (t_{n}) + h \cdot \frac{d y}{d t}_{t = t_{n}} = y (t_{n}) + h \cdot f (t_{n}, y (t_{n}))

The image below provides a fantastic summary of this process. The blue curve is the true, unknown solution. The black piecewise line is our approximation, created by taking small steps along the tangent at the start of each interval.

A graphical representation of Euler's method. The approximation starts at $(x_0, y_0)$ and takes sequential linear steps, each determined by the slope of the curve at the beginning of the step. The formulas for the ODE and the iterative update are shown on the right.

Applying Euler's Method to Kolmogorov's Equations

Now, let's connect this back to our time-inhomogeneous Markov models. The Kolmogorov forward equation gives us the derivative we need:

\frac{d}{d t} (_{t} P_{x}) = (_{t} P_{x}) Γ_{x + t}

This is a matrix equation, but it represents a system of individual differential equations, one for each transition probability $_{t} p_{x}^{ij}$ . The derivative of each probability is a function of the other probabilities and the transition intensities at that time.

We can apply the Euler approximation directly. The derivative $\frac{d}{d t} p_{ij} (t)$ can be approximated using a finite difference:

\frac{d}{d t}_{t} p_{x}^{ij} \approx \frac{_{t + h} p _{x}^{ij} - _{t} p _{x}^{ij}}{h}

By substituting this into the Kolmogorov forward equation and rearranging, we get the Euler update rule for our transition probabilities:

_{t + h} p_{x}^{ij} \approx_{t} p_{x}^{ij} + h \cdot (\frac{d}{d t}_{t} p_{x}^{ij})

Let's see this in action with a concrete example.

[PDF] Multiple State Models

This document from Michigan State University provides lecture notes on multiple state models and explicitly shows how Euler's method is used to find numerical solutions.

Please read the section titled "Numerical evaluation of transition probabilities". Focus on the text from this paragraph. It shows exactly how the approximation of the derivative leads to the recursive, forward-marching solution.

The process is as follows:

Start at time $t = 0$ with the known initial conditions: $_{0} p_{x}^{ii} = 1$ for all states $i$ , and $_{0} p_{x}^{ij} = 0$ for all $i \neq = j$ .
Choose a small step size $h$ (e.g., $h = 1/12$ for monthly steps, or $h = 1$ for yearly steps).
Calculate the probabilities at time $t = h$ by applying the Euler update rule to the probabilities at $t = 0$ . You will need the values of the transition intensities $μ_{x + t}^{ij}$ at $t = 0$ .
Calculate the probabilities at time $t = 2 h$ by applying the Euler update rule to the probabilities you just found at $t = h$ . You will now need the transition intensities at $t = h$ .
Repeat this process, "marching forward" in time until you reach the desired duration.

A Worked Example: Health-Sickness Model

Let's solidify this with a full example. We'll use the three-state healthy(0)-sick(1)-dead(2) model, which is a classic actuarial application.

[PDF] Multiple State Models

The same lecture notes provide a detailed walk-through of a numerical solution for a health-sickness model.

Please study the material from the section titled "Illustrative example from book" through "Detailed results with step size h = 1/12". Specifically: Examine the numerical process. See how the general Kolmogorov forward equations are written for this specific 3-state model, and then how they are immediately converted into Euler's approximation formulas. Carefully study the table of detailed results. This table is the output of the iterative process.

Let's manually compute the first step of that table to ensure we understand the mechanics. The goal is to find the probabilities at time $t = 1/12$ for a person aged 60, given they start in the Healthy state (state 0).

Initial Conditions (t=0):

$_{0} p_{60}^{00} = 1$
$_{0} p_{60}^{01} = 0$
$_{0} p_{60}^{02} = 0$

Intensities at age 60+0 = 60 (from the table):

$μ_{60}^{01} = 0.01420$
$μ_{60}^{02} = 0.01495$
$μ_{60}^{10} = 0.00142$
$μ_{60}^{12} = 0.01495$

Step size: $h = 1/12$

Now, we apply the Euler formulas from the PDF:

1. Calculate $_{1/12} p_{60}^{00}$ :

_{1/12} p_{60}^{00} \approx_{0} p_{60}^{00} + h [_{0} p_{60}^{01} μ_{60}^{10} -_{0} p_{60}^{00} (μ_{60}^{01} + μ_{60}^{02})]

\approx 1 + \frac{1}{12} [0 \cdot (0.00142) - 1 \cdot (0.01420 + 0.01495)]

\approx 1 + \frac{1}{12} [- 0.02915] \approx 1 - 0.002429 = 0.997571

This matches the value 0.99757 in the table for $t = 1/12$ .

2. Calculate $_{1/12} p_{60}^{01}$ :

_{1/12} p_{60}^{01} \approx_{0} p_{60}^{01} + h [_{0} p_{60}^{00} μ_{60}^{01} -_{0} p_{60}^{01} (μ_{60}^{10} + μ_{60}^{12})]

\approx 0 + \frac{1}{12} [1 \cdot (0.01420) - 0 \cdot (0.00142 + 0.01495)]

\approx \frac{1}{12} [0.01420] \approx 0.001183

This matches the value 0.00118 in the table.

3. Calculate $_{1/12} p_{60}^{02}$ :
We could use the Euler formula, but it's simpler to use the fact that probabilities must sum to 1:

_{1/12} p_{60}^{02} \approx 1 -_{1/12} p_{60}^{00} -_{1/12} p_{60}^{01} \approx 1 - 0.997571 - 0.001183 = 0.001246

This matches the value 0.00125 in the table.

To get the probabilities for $t = 2/12$ , you would repeat this exact process, but using the probabilities at $t = 1/12$ as your starting point and the intensities for age $60 + 1/12$ .

Practice on an Exam-Style Question

Euler's method is not just a theoretical tool; it is directly tested on the CS2 exam. The following problem demonstrates a typical format.

[PDF] Advanced Long-Term Actuarial Mathematics Sample Questions

This document contains sample questions from the Society of Actuaries. We will focus on Question 7, which models a television warranty.

Please read the setup for Question 7. Note the 3-state model, the time-dependent transition forces, and the instruction that the probabilities are calculated using the Forward Euler approximation with a step size of h=0.5. Then, review part (b)(i), which asks you to calculate a specific probability.

Problem: Using the student's approach (Euler's method with $h = 0.5$ ), calculate $_{0.5} p_{0}^{01}$ .

Solution:
The Kolmogorov forward equation for the probability of being in state 1 (Under Repair) at time $t$ , given starting in state 0 (Working), is:

\frac{d}{d t}_{t} p_{0}^{01} =_{t} p_{0}^{00} \cdot μ_{t}^{01} -_{t} p_{0}^{01} \cdot (μ_{t}^{10} + μ_{t}^{12})

We apply the Euler approximation for one step from $t = 0$ to $t = 0.5$ :

_{0.5} p_{0}^{01} \approx_{0} p_{0}^{01} + h \cdot [_{0} p_{0}^{00} \cdot μ_{0}^{01} -_{0} p_{0}^{01} \cdot (μ_{0}^{10} + μ_{0}^{12})]

From the problem, we have:

Initial conditions: $_{0} p_{0}^{01} = 0$ , $_{0} p_{0}^{00} = 1$
Step size: $h = 0.5$
Intensities at $t = 0$ : $μ_{0}^{01} = 0.5$ , $μ_{0}^{10} = 0.6$ , $μ_{0}^{12} = 0.2$

Substituting these values:

_{0.5} p_{0}^{01} \approx 0 + 0.5 \cdot [1 \cdot (0.5) - 0 \cdot (0.6 + 0.2)]

_{0.5} p_{0}^{01} \approx 0.5 \cdot 0.5 = 0.25

This calculation demonstrates how to apply the method to a specific question. Part (d) of the question asks how to improve accuracy. The simplest answer is to decrease the step size $h$ . A smaller step size means our linear approximations track the true curve more closely, reducing the error at each step.

Conclusion

Today we have bridged the gap between the theory of time-inhomogeneous processes and their practical application. While the governing Kolmogorov equations can be complex, we have seen that a straightforward numerical technique, Euler's method, allows us to approximate their solutions effectively.

Key Takeaways:

Euler's method provides a numerical solution to differential equations by taking small, linear steps forward in time.
The update rule is $p (t + h) \approx p (t) + h \cdot (rate of change at t)$ .
To apply this to Kolmogorov's equations, we need the initial conditions ( $p (0)$ ), the step size ( $h$ ), and the formulas for the time-dependent transition intensities ( $μ_{x + t}^{ij}$ ).
The process is iterative: we calculate probabilities for each state at $t + h$ using the known values at $t$ , and repeat this until we reach our target duration.
The accuracy of the approximation is improved by using a smaller step size $h$ .

In our next lesson, we will briefly step back from the full multi-state model to explore the core concepts of survival analysis in greater detail. We will solidify the fundamental relationships between the survival function, the hazard function, and the cumulative hazard, which are the building blocks for all the models we have been discussing.

Recommended resources