4. The Kolmogorov Backward Equation#
In addition to what’s in Anaconda, this lecture will need the following libraries:
!pip install quantecon
Show code cell output
Requirement already satisfied: quantecon in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (0.7.2)
Requirement already satisfied: numba>=0.49.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (0.60.0)
Requirement already satisfied: numpy>=1.17.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (1.26.4)
Requirement already satisfied: requests in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (2.32.3)
Requirement already satisfied: scipy>=1.5.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (1.13.1)
Requirement already satisfied: sympy in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (1.13.2)
Requirement already satisfied: llvmlite<0.44,>=0.43.0dev0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from numba>=0.49.0->quantecon) (0.43.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (3.7)
Requirement already satisfied: urllib3<3,>=1.21.1 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (2.2.3)
Requirement already satisfied: certifi>=2017.4.17 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (2024.8.30)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from sympy->quantecon) (1.3.0)
4.1. Overview#
As models become more complex, deriving analytical representations of the
Markov semigroup
This is analogous to the idea that solutions to continuous time models often lack analytical solutions.
For example, when studying deterministic paths in continuous time, infinitesimal descriptions (ODEs and PDEs) are often more intuitive and easier to write down than the associated solutions.
(This is one of the shining insights of mathematics, beginning with the work of great scientists such as Isaac Newton.)
We will see in this lecture that the same is true for continuous time Markov chains.
To help us focus on intuition in this lecture, rather than technicalities, the state space is assumed to be finite, with
Later we will investigate the case where
We will use the following imports
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import quantecon as qe
from numba import njit
from scipy.linalg import expm
from scipy.stats import binom
4.2. State Dependent Jump Intensities#
As we have seen, continuous time Markov chains jump between states, and hence can have the form
where
(We are assuming that
In the previous lecture,
the sequence
was drawn from a Markov matrix and called the embedded jump chain, whilethe holding times
were IID and Exp for some constant jump intensity .
In this lecture, we will generalize by allowing the jump intensity to vary with the state.
This difference sounds minor but in fact it will allow us to reach full generality in our description of continuous time Markov chains, as clarified below.
4.2.1. Motivation#
As a motivating example, recall the inventory model, where we assumed that the wait time for the next customer was equal to the wait time for new inventory.
This assumption was made purely for convenience and seems unlikely to hold true.
When we relax it, the jump intensities depend on the state.
4.2.2. Jump Chain Algorithm#
We start with three primitives
An initial condition
,a Markov matrix
on satisfying for all anda function
mapping to .
The process
starts at state
, draw from ,waits there for an exponential time
with rate and thenupdates to a new state
drawn from .
Now we take
Here is the same algorithm written more explicitly:
(Jump Chain Algorithm)
Inputs
Outputs Markov chain
Draw
from , set and .Draw
independently from Exp .Set
.Set
for in .Draw
from .Set
and go to step 2.
The sequence
The restriction
4.3. Computing the Semigroup#
For the jump process
The approach we adopt is
Use probabilistic reasoning to obtain an integral equation that the semigroup must satisfy.
Convert the integral equation into a differential equation that is easier to work with.
Solve this differential equation to obtain the Markov semigroup
.
The differential equation in question has a special name: the Kolmogorov backward equation.
4.3.1. An Integral Equation#
Here is the first step in the sequence listed above.
(An Integral Equation)
The semigroup
for all
Here
4.3.2. Kolmogorov’s Differential Equation#
We have now confirmed that the semigroup
Equation (4.1) is important but we can simplify it further without losing information by taking the time derivative.
This leads to our main result for the lecture
(Kolmogorov Backward Equation)
The semigroup
The derivative on the left hand side of (4.4) is taken element by
element, with respect to
The proof that differentiating (4.1) yields (4.4) is an important exercise (see below).
4.3.3. Exponential Solution#
The Kolmogorov backward equation is a matrix-valued differential equation.
Recall that, for a scalar differential equation
This, along with
where the right hand side is the matrix exponential, with definition
Working element by element, it is not difficult to confirm that
the derivative of the exponential function
Hence, differentiating (4.5) gives
Notice that our solution
for the semigroup of the jump process
In particular, we showed that, for the model with
constant jump intensity
This is obviously a special case of (4.8).
4.4. Properties of the Solution#
Let’s investigate further the properties of the exponential solution.
4.4.1. Checking the Transition Semigroup Properties#
While we have confirmed that
(From Jump Chain to Semigroup)
Let
Proof. Observe first that
As a small exercise, you can check that, with
This implies that
In other words, each
Next we check nonnegativity of all elements of
To this end, adopting an argument from [Stroock, 2013], we
set
It is not difficult to check that
Recalling that, for matrix exponentials,
It is clear from this representation that all entries of
Finally, we need to check the continuity condition
We can now be reassured that our solution to the Kolmogorov backward equation is indeed a Markov semigroup.
4.4.2. Uniqueness#
Might there be another, entirely different Markov semigroup that also satisfies the Kolmogorov backward equation?
The answer is no: linear ODEs in finite dimensional vector space with
constant coefficients and fixed initial conditions (in this case
In fact it’s not hard to supply a proof — see the exercises.
4.5. Application: The Inventory Model#
Let us look at a modified version of the inventory model where jump intensities depend on the state.
In particular, the wait time for new inventory will now be exponential at rate
The arrival rate for customers will still be denoted by
For parameters we take
α = 0.6
λ = 0.5
γ = 0.1
b = 10
Our plan is to investigate the distribution
We will do this by simulating many independent draws of
(In the exercises you are asked to calculate
@njit
def draw_X(T, X_0, max_iter=5000):
"""
Generate one draw of X_T given X_0.
"""
J, Y = 0, X_0
m = 0
while m < max_iter:
s = 1/γ if Y == 0 else 1/λ
W = np.random.exponential(scale=s) # W ~ E(λ)
J += W
if J >= T:
return Y
# Otherwise update Y
if Y == 0:
Y = b
else:
U = np.random.geometric(α)
Y = Y - min(Y, U)
m += 1
@njit
def independent_draws(T=10, num_draws=100):
"Generate a vector of independent draws of X_T."
draws = np.empty(num_draws, dtype=np.int64)
for i in range(num_draws):
X_0 = np.random.binomial(b+1, 0.25)
draws[i] = draw_X(T, X_0)
return draws
T = 30
n = b + 1
draws = independent_draws(T, num_draws=100_000)
fig, ax = plt.subplots()
ax.bar(range(n), [np.mean(draws == i) for i in range(n)], width=0.8, alpha=0.6)
ax.set_xlabel("inventory", fontsize=14)
plt.show()
If you experiment with the code above, you will see that the large amount of
mass on zero is due to the low arrival rate
4.6. Exercises#
In the discussion above, we generated an approximation of
α = 0.6
λ = 0.5
γ = 0.1
b = 10
The calculation was done by simulating independent draws and histogramming.
Try to generate the same figure using (4.8) instead, modifying code from our lecture on the Markov property.
We claimed above that the solution
Try to supply a proof.
(This is not an easy exercise but worth thinking about in any case.)
4.7. Solutions#
Note
code is currently not supported in sphinx-exercise
so code-cell solutions are immediately after this
solution block.
Solution to Exercise 4.1
Here is one solution:
α = 0.6
λ = 0.5
γ = 0.1
b = 10
states = np.arange(n)
I = np.identity(n)
# Embedded jump chain matrix
K = np.zeros((n, n))
K[0, -1] = 1
for i in range(1, n):
for j in range(0, i):
if j == 0:
K[i, j] = (1 - α)**(i-1)
else:
K[i, j] = α * (1 - α)**(i-j-1)
# Jump intensities as a function of the state
r = np.ones(n) * λ
r[0] = γ
# Q matrix
Q = np.empty_like(K)
for i in range(n):
for j in range(n):
Q[i, j] = r[i] * (K[i, j] - I[i, j])
def P_t(ψ, t):
return ψ @ expm(t * Q)
ψ_0 = binom.pmf(states, n, 0.25)
ψ_T = P_t(ψ_0, T)
fig, ax = plt.subplots()
ax.bar(range(n), ψ_T, width=0.8, alpha=0.6)
ax.set_xlabel("inventory", fontsize=14)
plt.show()
Solution to Exercise 4.2
One can easily verify that, when
Note also that, with the change of variable
Applying (4.10) yields
After minor rearrangements this becomes
which is identical to (4.4).
Solution to Exercise 4.3
Here is one proof of uniqueness.
Suppose that
Fix
Note that
Note also that
where, in the second last equality, we used (4.7).
Hence
Since