Numerical Integration of Time-Varying Kolmogorov Equations

Why Numerical Integration Is Needed

When transition intensities are constant over an interval, matrix methods can give transition probabilities directly. If the intensities change continuously with time, however, the forward equation generally has to be solved as a differential equation. Numerical integration provides an approximation by advancing the state probabilities through a sequence of short time steps.

The calculation is simple in principle, but its accuracy depends on using the intensity matrix at the right time, matching the time unit to the intensity unit, and checking the resulting probabilities. A numerical vector can sum to one and still contain a negative entry; conversely, a vector with all nonnegative entries can still be wrong if the initial distribution or time argument was misread.

Definition: Let $p(t)$ be a row vector of state probabilities at time $t$, and let $Q(t)$ be the time-dependent generator matrix. The forward Kolmogorov equation is $$ \frac{d}{dt}p(t)=p(t)Q(t), \qquad p(0)=p_0. $$ For distinct states $i$ and $j$, $q_{ij}(t)$ is the transition intensity from $i$ to $j$ at time $t$. The diagonal entry is $q_{ii}(t)=-\sum_{j\ne i}q_{ij}(t)$, so each row of $Q(t)$ sums to zero.

This tutorial uses row vectors, so probabilities are multiplied by the generator on the right. The initial vector $p_0$ must have nonnegative entries summing to one. Intensities must use the same time unit as the integration step: an intensity measured per year must be paired with a step measured in years. The method applies when the intensities are specified functions of time and the forward equation is defined over the interval being integrated. If an intensity changes discontinuously, split the calculation at the change point and use the appropriate intensity on each side.

Explicit Euler Integration

The explicit Euler method approximates the derivative over a short interval by its value at the beginning of that interval. For a step from $t_n$ to $t_{n+1}=t_n+h$, it gives

p_{n+1}=p_n+h\,p_nQ(t_n)=p_n\bigl(I+hQ(t_n)\bigr).

Here $h$ is the step length and $I$ is the identity matrix. In a time-varying model, the intensity matrix is evaluated at the current time $t_n$, not held fixed at its value at time zero. Repeating the update approximates the solution over a longer interval.

A useful way to interpret the update is to track probability flow. In state $i$, probability leaves at the total outgoing intensity and enters from other states at their corresponding transition intensities. The vector $p_nQ(t_n)$ gives the instantaneous net change in each state probability. Euler adds $h$ times that change to the current vector.

Step-size check: If $p_n$ is a probability vector and every entry of $I+hQ(t_n)$ is nonnegative, then the Euler update is also a probability vector. Since the rows of $Q(t_n)$ sum to zero, the rows of $I+hQ(t_n)$ sum to one. Nonnegative entries require, in particular, $h\sum_{j\ne i}q_{ij}(t_n)\le 1$ for each state $i$. If this condition fails, Euler can produce negative probabilities even though the underlying model is valid.

Worked Example: Comparing Two Euler Step Sizes

A hypothetical two-state model has a healthy state 0 and a disabled state 1. The intensity from state 0 to state 1 is $q_{01}(t)=0.20+0.10t$ per year, where $t$ is measured in years. The recovery intensity is $q_{10}(t)=0.10$ per year. At time zero the person is certainly healthy, so $p(0)=(1,0)$. Use left-endpoint Euler steps to estimate the distribution at one year, first with $h=0.50$ years and then with $h=0.25$ years. Give final probabilities to four decimal places.

The generator is $$ Q(t)= \begin{pmatrix} -(0.20+0.10t)&0.20+0.10t\\ 0.10&-0.10 \end{pmatrix}. $$ For a probability vector $p=(p_0,p_1)$, the derivative of the disabled-state probability is $p_0q_{01}-p_1q_{10}$.

Method and $h=0.50$. At $t=0$, $q_{01}=0.20$. Starting from $(1,0)$, the derivative is $(-0.20,0.20)$. The first update is $$ (1,0)+0.50(-0.20,0.20)=(0.90,0.10). $$ At $t=0.50$, $q_{01}=0.20+0.10(0.50)=0.25$. The derivative at $(0.90,0.10)$ is $$ (-0.90(0.25)+0.10(0.10),\;0.90(0.25)-0.10(0.10))=(-0.215,0.215). $$ Thus the second update gives $$ (0.90,0.10)+0.50(-0.215,0.215)=(0.7925,0.2075). $$

Method and $h=0.25$. At $t=0$, the derivative is $(-0.20,0.20)$, giving $(0.95,0.05)$. At $t=0.25$, $q_{01}=0.225$, so the derivative is $(-0.95(0.225)+0.05(0.10),\;0.95(0.225)-0.05(0.10))=(-0.20875,0.20875)$. The next vector is $(0.8978125,0.1021875)$.

At $t=0.50$, $q_{01}=0.25$. The derivative is $(-0.8978125(0.25)+0.1021875(0.10),\;0.8978125(0.25)-0.1021875(0.10))=(-0.214234375,0.214234375)$. The updated vector is $(0.84425390625,0.15574609375)$.

At $t=0.75$, $q_{01}=0.275$. The derivative is $(-0.84425390625(0.275)+0.15574609375(0.10),\;0.84425390625(0.275)-0.15574609375(0.10))=(-0.21659521484375,0.21659521484375)$. The one-year estimate is therefore $(0.7901051025390625,0.2098948974609375)$, or $(0.7901,0.2099)$ to four decimals.

Check and interpretation. Both final vectors sum to 1 and have nonnegative entries. The disabled-state estimates are 0.2075 with two steps and 0.2099 with four steps. The smaller-step result is a refinement, not a proof of exactness. The main trap avoided is using the initial intensity $0.20$ for every step instead of evaluating $q_{01}$ at each step's starting time.

Why the Basic Checks Work

The Euler formula follows by a first-order Taylor approximation. Over a small interval of length $h$, the solution satisfies $p(t+h)=p(t)+h p'(t)+O(h^2)$, where $O(h^2)$ denotes terms whose magnitude is of second order in $h$. Substituting the forward equation $p'(t)=p(t)Q(t)$ gives the Euler update, with a local error of order $h^2$. Over a fixed interval, the accumulated error is generally of first order in the step length, provided the equation is sufficiently regular.

The probability-sum check follows directly from the zero row sums of the generator. Let $\mathbf{1}$ be a column vector of ones. Then $Q(t)\mathbf{1}=0$, and an Euler step gives

p_{n+1}\mathbf{1} =p_n\mathbf{1}+h\,p_nQ(t_n)\mathbf{1} =p_n\mathbf{1}.

Thus a correctly formed generator preserves the total probability in exact arithmetic at each Euler step. This argument does not establish nonnegativity: the step-size condition must also be checked. For numerical work, small discrepancies in a sum of one can arise from rounding, but a substantial discrepancy usually signals an arithmetic, matrix-orientation, or generator-entry error.

Specify the model and units.
Write the states, the time-dependent off-diagonal intensities, and the time unit. Form each diagonal entry as the negative sum of the other entries in its row.

Check the initial condition.
Confirm that the initial probabilities are nonnegative and sum to one. Do not replace a stated mixed initial distribution by certainty in one state.

Advance one step at a time.
Evaluate $Q$ at the current time, calculate $p_nQ(t_n)$, and then multiply the derivative by the step length.

Validate and refine.
Check the row sums of $Q$, the sum and signs of the probability vector after each update, and compare with a calculation using smaller steps.

Worked Example: A Three-State Model with a Mixed Initial Distribution

Consider a hypothetical sequential model with states 0, 1, and 2. Transitions are possible only from 0 to 1 and from 1 to 2; state 2 is absorbing. The intensity from 0 to 1 is $0.20+0.10t$ per year, and the intensity from 1 to 2 is $0.40$ per year. At time zero, the distribution is $p(0)=(0.80,0.20,0)$. Use two left-endpoint Euler steps of $h=0.50$ years to estimate $p(1)$.

The generator and flow equations are $$ Q(t)= \begin{pmatrix} -(0.20+0.10t)&0.20+0.10t&0\\ 0&-0.40&0.40\\ 0&0&0 \end{pmatrix}, \qquad p'(t)=\bigl(-p_0q_{01},\;p_0q_{01}-0.40p_1,\;0.40p_1\bigr). $$

At $t=0$, $q_{01}=0.20$. From $(0.80,0.20,0)$, the derivative is $(-0.80(0.20),\;0.80(0.20)-0.20(0.40),\;0.20(0.40))=(-0.16,0.08,0.08)$. Multiplying by $0.50$ and adding gives $$ (0.80,0.20,0)+0.50(-0.16,0.08,0.08)=(0.72,0.24,0.04). $$

At $t=0.50$, $q_{01}=0.20+0.10(0.50)=0.25$. The derivative at $(0.72,0.24,0.04)$ is $(-0.72(0.25),\;0.72(0.25)-0.24(0.40),\;0.24(0.40))=(-0.18,0.084,0.096)$. The second update is $$ (0.72,0.24,0.04)+0.50(-0.18,0.084,0.096)=(0.63,0.282,0.088). $$

Check and interpretation. The entries sum to $0.63+0.282+0.088=1$, and all are nonnegative. State 2 has probability $0.088$ under this approximation. The mixed initial vector matters: the initial probability in state 1 contributes to transitions into state 2 during the first step. Treating the starting population as entirely in state 0 would change every subsequent flow.

Improving an Approximation and Recognizing Failure

A smaller Euler step usually reduces discretization error, but it requires more calculations. A predictor-corrector method, also called Heun's method, uses the derivative at both ends of a step. It first predicts the endpoint with Euler, then averages the derivative at the start and the predicted endpoint:

\begin{aligned} f(t,p)&=pQ(t),\\ p^*&=p_n+h f(t_n,p_n),\\ p_{n+1}&=p_n+\frac{h}{2}\left[f(t_n,p_n)+f(t_n+h,p^*)\right]. \end{aligned}

For a sufficiently smooth differential equation, this method has second-order global accuracy, compared with first-order global accuracy for Euler. The predictor $p^*$ is an intermediate estimate used to evaluate the endpoint derivative; it is not automatically a valid probability vector. The corrected result must still be checked for nonnegative entries and a sum of one. Higher accuracy does not remove the need for probability checks.

Worked Example: Detecting an Invalid Euler Step

A hypothetical two-state model has intensity $q_{01}(t)=1.50+0.50t$ per year and recovery intensity $q_{10}(t)=0.50$ per year. The individual starts in state 0 with probability one. First attempt one Euler step of $h=0.75$ years from time zero; then use three steps of $h=0.25$ years to estimate the distribution at time $0.75$. Give final values to four decimals.

At time zero, the outgoing intensity from state 0 is $1.50$. Its derivative vector from $(1,0)$ is $(-1.50,1.50)$. The single coarse step gives $$ (1,0)+0.75(-1.50,1.50)=(-0.125,1.125). $$ The entries sum to one, but the negative state-0 probability makes this an invalid probability approximation.

Now use $h=0.25$. At $t=0$, the update is $(1,0)+0.25(-1.50,1.50)=(0.625,0.375)$.

At $t=0.25$, the forward intensity is $1.50+0.50(0.25)=1.625$. The derivative at $(0.625,0.375)$ is $$ (-0.625(1.625)+0.375(0.50),\;0.625(1.625)-0.375(0.50))=(-0.828125,0.828125). $$ The next vector is $(0.625,0.375)+0.25(-0.828125,0.828125)=(0.41796875,0.58203125)$.

At $t=0.50$, the forward intensity is $1.75$. The derivative is $$ (-0.41796875(1.75)+0.58203125(0.50),\;0.41796875(1.75)-0.58203125(0.50))=(-0.4404296875,0.4404296875). $$ The last update gives $(0.41796875,0.58203125)+0.25(-0.4404296875,0.4404296875)=(0.307861328125,0.692138671875)$, or $(0.3079,0.6921)$ to four decimals.

Check and interpretation. The refined result is nonnegative and sums to one. At each left endpoint in the refined calculation, $h$ times the largest outgoing intensity is at most $0.25(1.75)=0.4375<1$, consistent with the Euler nonnegativity check. The coarse-step error is not fixed by renormalizing $(-0.125,1.125)$; the negative value indicates that the step should be reduced.

Common Mistakes and Practical Checks

The most important checks are complementary: none should be used alone as proof that the calculation is correct.

Check	What to verify	What a failure may indicate
Initial condition	Entries are nonnegative and sum to one	A missing state, transcription error, or misunderstood starting distribution
Generator	Off-diagonal entries are nonnegative intensities; each row sums to zero	A wrong diagonal or misplaced transition intensity
Time and rate units	Step length and intensity use matching units	For example, using a monthly step with annual rates without converting units
Time argument	Use $Q(t_n)$ for a left-endpoint Euler step	Holding rates fixed or evaluating them at the wrong time
Updated probabilities	Entries are nonnegative and sum to one, up to rounding	A step that is too large, a sign error, or an arithmetic error

A common exam trap is to see a vector whose entries sum to one and conclude that the numerical step is valid. In the coarse step above, $-0.125+1.125=1$, but the negative entry is impossible as a probability. Another trap is to confuse a mixed initial distribution with a conditional starting state. If a problem gives $p(0)=(0.80,0.20,0)$, the forward equation must be initialized with that entire vector.

Key takeaway: A numerical solution is credible only when its setup, time units, generator, initial condition, and resulting probabilities are all checked. Refine the step size when results change materially or Euler produces a negative entry; do not repair an invalid vector by forcing it to sum to one.

Check Your Understanding

For each problem, use the stated row-vector convention. The final answers are provided for checking.

A two-state model has intensities $q_{01}(t)=0.20+0.20t$ per year and $q_{10}(t)=0.10$ per year. Starting from $p(0)=(0.60,0.40)$, use one left-endpoint Euler step of $0.50$ years. State the resulting distribution. (Answer: $(0.5600,0.4400)$)
A sequential three-state model has $q_{01}(t)=0.30+0.20t$ per year and $q_{12}=0.50$ per year; state 2 is absorbing. Starting from $(1,0,0)$, use two Euler steps of $0.50$ years. What is the approximate distribution at year 1? (Answer: $(0.6800,0.2825,0.0375)$)
A constant two-state generator has $q_{01}=1.20$ and $q_{10}=0.40$ per year. Starting in state 0, calculate the Euler result for $h=1$ year and decide whether it is a valid probability vector. Then state the largest $h$ satisfying the one-step nonnegativity condition for both states. (Answer: the $h=1$ result is $(-0.2000,1.2000)$, invalid; the largest permitted step is $1/1.20=0.8333$ years.)
Use one Heun step with $h=1$ year for a two-state model with $q_{01}(t)=0.20+0.10t$, $q_{10}=0.10$, and $p(0)=(1,0)$. What is the corrected estimate at year 1? (Answer: $(0.7900,0.2100)$)
A generator is $Q=\begin{pmatrix}-0.30&0.30\\0.20&-0.20\end{pmatrix}$, with rates per year. Starting from $(0.25,0.75)$, apply one Euler step of $0.50$ years. Check both the generator row sums and the resulting probability vector. (Answer: both row sums are zero; the updated vector is $(0.2875,0.7125)$, which is nonnegative and sums to one.)

What You'll Learn