kolmogorov_bwd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cc5c1148",
   "metadata": {},
   "source": [
    "# The Kolmogorov Backward Equation\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aa7a437",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "!pip install quantecon\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e72e9303",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "As models become more complex, deriving analytical representations of the\n",
    "Markov semigroup $ (P_t) $ becomes harder.\n",
    "\n",
    "This is analogous to the idea that solutions to continuous time models often\n",
    "lack analytical solutions.\n",
    "\n",
    "For example, when studying deterministic paths in continuous time,\n",
    "infinitesimal descriptions ([ODEs](https://en.wikipedia.org/wiki/Ordinary_differential_equation) and [PDEs](https://en.wikipedia.org/wiki/Partial_differential_equation)) are often more intuitive and easier\n",
    "to write down than the associated solutions.\n",
    "\n",
    "(This is one of the shining insights of mathematics, beginning with the work\n",
    "of great scientists such as Isaac Newton.)\n",
    "\n",
    "We will see in this lecture that the same is true for continuous time Markov\n",
    "chains.\n",
    "\n",
    "To help us focus on intuition in this lecture, rather than technicalities, the state space is assumed to be finite, with $ |S|=n $.\n",
    "\n",
    "Later we will investigate the case where $ |S| = \\infty $.\n",
    "\n",
    "We will use the following imports"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "760be8e7",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import quantecon as qe\n",
    "from numba import njit\n",
    "\n",
    "from scipy.linalg import expm\n",
    "from scipy.stats import binom\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e5d214b",
   "metadata": {},
   "source": [
    "\n",
    "<a id='sdji'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb7176d",
   "metadata": {},
   "source": [
    "## State Dependent Jump Intensities\n",
    "\n",
    "As we have seen, continuous time Markov chains jump between states, and hence can\n",
    "have the form\n",
    "\n",
    "$$\n",
    "X_t = \\sum_{k \\geq 0} Y_k \\mathbb 1\\{J_k \\leq t < J_{k+1}\\}\n",
    "    \\qquad (t \\geq 0)\n",
    "$$\n",
    "\n",
    "where $ (J_k) $ are jump times and $ (Y_k) $ are the states at each jump.\n",
    "\n",
    "(We are assuming that $ J_k \\to \\infty $ with probability one, so that $ X_t $ is well\n",
    "defined for all $ t \\geq 0 $, but this is always true for when holding times are exponential and the state space is finite.)\n",
    "\n",
    "In the [previous lecture](markov_prop.ipynb),\n",
    "\n",
    "- the sequence $ (Y_k) $ was drawn from a Markov matrix $ K $ and called the embedded jump chain, while  \n",
    "- the holding times $ W_k := J_k - J_{k-1} $ were IID and Exp$ (\\lambda) $ for some\n",
    "  constant jump intensity $ \\lambda $.  \n",
    "\n",
    "\n",
    "In this lecture, we will generalize by allowing the jump intensity to vary\n",
    "with the state.\n",
    "\n",
    "This difference sounds minor but in fact it will allow us to reach full generality\n",
    "in our description of continuous time Markov chains, as\n",
    "clarified below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ea1094a",
   "metadata": {},
   "source": [
    "### Motivation\n",
    "\n",
    "As a motivating example, recall [the inventory model](markov_prop.ipynb#inventory-dynam),\n",
    "where we assumed that the wait time for the next customer was equal\n",
    "to the wait time for new inventory.\n",
    "\n",
    "This assumption was made purely for convenience and seems unlikely to hold true.\n",
    "\n",
    "When we relax it, the jump intensities depend on the state.\n",
    "\n",
    "\n",
    "<a id='jumpchainalgo'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5ee12a0",
   "metadata": {},
   "source": [
    "### Jump Chain Algorithm\n",
    "\n",
    "We start with three primitives\n",
    "\n",
    "1. An initial condition $ \\psi $,  \n",
    "1. a Markov matrix $ K $ on $ S $ satisfying $ K(x, x) = 0 $ for all $ x \\in S $\n",
    "  and  \n",
    "1. a function $ \\lambda $ mapping $ S $ to $ (0, \\infty) $.  \n",
    "\n",
    "\n",
    "The process $ (X_t) $\n",
    "\n",
    "- starts at state $ x $, draw from $ \\psi $,  \n",
    "- waits there for an exponential time $ W $ with rate $ \\lambda(x) $ and then  \n",
    "- updates to a new state $ y $ drawn from $ K(x, \\cdot) $.  \n",
    "\n",
    "\n",
    "Now we take $ y $ as the new state for the process and repeat.\n",
    "\n",
    "Here is the same algorithm written more explicitly:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf15f21a",
   "metadata": {},
   "source": [
    "###  (Jump Chain Algorithm)\n",
    "\n",
    "**Inputs** $ \\psi \\in \\dD $, rate function $ \\lambda $, Markov matrix $ K $\n",
    "\n",
    "**Outputs** Markov chain $ (X_t) $\n",
    "\n",
    "1. Draw $ Y_0 $ from $ \\psi $, set $ J_0 = 0 $ and $ k=1 $.  \n",
    "1. Draw $ W_k $ independently from Exp$ (\\lambda(Y_{k-1})) $.  \n",
    "1. Set $ J_k = J_{k-1} + W_k $.  \n",
    "1. Set $ X_t = Y_{k-1} $ for $ t $ in $ [J_{k-1}, J_k) $.  \n",
    "1. Draw $ Y_k $ from $ K(Y_{k-1}, \\cdot) $.  \n",
    "1. Set $ k = k+1 $ and go to step 2.  \n",
    "\n",
    "\n",
    "The sequence $ (W_k) $ is drawn as an IID sequence and $ (W_k) $ and $ (Y_k) $ are\n",
    "drawn independently.\n",
    "\n",
    "The restriction $ K(x,x) = 0 $ for all $ x $ implies that $ (X_t) $ actually jumps at each jump time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a32faf6b",
   "metadata": {},
   "source": [
    "## Computing the Semigroup\n",
    "\n",
    "For the jump process $ (X_t) $ with time varying intensities described in the\n",
    "jump chain algorithm, calculating the Markov semigroup is not a trivial exercise.\n",
    "\n",
    "The approach we adopt is\n",
    "\n",
    "1. Use probabilistic reasoning to obtain an integral equation that the\n",
    "  semigroup must satisfy.  \n",
    "1. Convert the integral equation into a differential equation that is easier\n",
    "  to work with.  \n",
    "1. Solve this differential equation to obtain the Markov semigroup $ (P_t) $.  \n",
    "\n",
    "\n",
    "The differential equation in question has a special name:  the Kolmogorov backward equation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca7b4bfd",
   "metadata": {},
   "source": [
    "### An Integral Equation\n",
    "\n",
    "Here is the first step in the sequence listed above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8749265",
   "metadata": {},
   "source": [
    "###  (An Integral Equation)\n",
    "\n",
    "The semigroup $ (P_t) $ of the jump chain with rate function $ \\lambda $ and Markov matrix $ K $ obeys the integral equation\n",
    "\n",
    "\n",
    "<a id='equation-kbinteg'></a>\n",
    "$$\n",
    "P_t(x, y) = e^{-t \\lambda(x)} I(x, y)\n",
    "    + \\lambda(x) \n",
    "      \\int_0^t (K P_{t-\\tau})(x, y) e^{- \\tau \\lambda(x)} d \\tau \\tag{4.1}\n",
    "$$\n",
    "\n",
    "for all $ t \\geq 0 $ and $ x, y $ in $ S $.\n",
    "\n",
    "Here $ (P_t) $ is the Markov semigroup of $ (X_t) $, the process constructed via\n",
    "[Algorithm 4.1](#ejc_algo), while $ K P_{t-\\tau} $ is the matrix product of $ K $ and\n",
    "$ P_{t-\\tau} $.\n",
    "\n",
    "Proof. Conditioning implicitly on $ X_0 = x $, the semigroup $ (P_t) $ must satisfy\n",
    "\n",
    "\n",
    "<a id='equation-pt-split'></a>\n",
    "$$\n",
    "P_t(x, y) \n",
    "    = \\PP\\{X_t = y\\}\n",
    "    = \\PP\\{X_t = y, \\; J_1 > t \\}\n",
    "        + \\PP\\{X_t = y, \\; J_1 \\leq t \\} \\tag{4.2}\n",
    "$$\n",
    "\n",
    "Regarding the first term on the right hand side of [(4.2)](#equation-pt-split), we have\n",
    "\n",
    "\n",
    "<a id='equation-pt-first'></a>\n",
    "$$\n",
    "\\PP\\{X_t = y, \\; J_1 > t \\}\n",
    "        = I(x, y) P\\{J_1 > t \\}\n",
    "        = I(x, y) e^{- t \\lambda(x)} \\tag{4.3}\n",
    "$$\n",
    "\n",
    "where $ I(x, y) = \\mathbb 1\\{x = y\\} $.\n",
    "\n",
    "For the second term on the right hand side of [(4.2)](#equation-pt-split), we have\n",
    "\n",
    "$$\n",
    "\\PP\\{X_t = y, \\; J_1 \\leq t \\}\n",
    "    = \\EE \n",
    "        \\left[\n",
    "            \\mathbb 1\\{J_1 \\leq t\\} \\PP\\{X_t = y \\,|\\, W_1, Y_1\\}\n",
    "        \\right]\n",
    "    = \\EE \n",
    "        \\left[\n",
    "            \\mathbb 1\\{J_1 \\leq t\\} P_{t - J_1} (Y_1, y) \n",
    "        \\right]\n",
    "$$\n",
    "\n",
    "Evaluating the expectation and using the independence of $ J_1 $ and $ Y_1 $, this becomes\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\PP\\{X_t = y, \\; J_1 \\leq t \\}\n",
    "    & = \\int_0^\\infty\n",
    "            \\mathbb 1\\{\\tau \\leq t\\}\n",
    "            \\sum_z K(x, z) P_{t - \\tau} (z, y)  \\lambda(x) e^{-\\tau \\lambda(x)} \n",
    "            d \\tau\n",
    "        \\\\\n",
    "    & = \\lambda(x)\n",
    "            \\int_0^t\n",
    "            \\sum_z K(x, z) P_{t - \\tau} (z, y)  e^{-\\tau \\lambda(x)} \n",
    "            d \\tau\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Combining this result with [(4.2)](#equation-pt-split) and [(4.3)](#equation-pt-first) gives\n",
    "[(4.1)](#equation-kbinteg)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "081f5532",
   "metadata": {},
   "source": [
    "### Kolmogorov’s Differential Equation\n",
    "\n",
    "We have now confirmed that the semigroup $ (P_t) $ associated with the jump\n",
    "chain process $ (X_t) $ satisfies [(4.1)](#equation-kbinteg).\n",
    "\n",
    "Equation [(4.1)](#equation-kbinteg) is important but we can simplify it further without\n",
    "losing information by taking the time derivative.\n",
    "\n",
    "This leads to our main result for the lecture"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca7b406d",
   "metadata": {},
   "source": [
    "###  (Kolmogorov Backward Equation)\n",
    "\n",
    "The semigroup $ (P_t) $ of the jump chain with rate function $ \\lambda $ and Markov matrix $ K $ satisfies the **Kolmogorov backward equation**\n",
    "\n",
    "\n",
    "<a id='equation-kolbackeq'></a>\n",
    "$$\n",
    "P'_t = Q P_t \n",
    "    \\quad \\text{where } \\;\n",
    "    Q(x, y) := \\lambda(x) (K(x, y) - I(x, y)) \\tag{4.4}\n",
    "$$\n",
    "\n",
    "The derivative on the left hand side of [(4.4)](#equation-kolbackeq) is taken element by\n",
    "element, with respect to $ t $, so that\n",
    "\n",
    "$$\n",
    "P'_t(x, y) = \\left( \\frac{d}{dt} P_t(x, y) \\right)\n",
    "    \\qquad ((x, y) \\in S \\times S)\n",
    "$$\n",
    "\n",
    "The proof that differentiating [(4.1)](#equation-kbinteg) yields [(4.4)](#equation-kolbackeq) is an\n",
    "important exercise (see below)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b172c7f",
   "metadata": {},
   "source": [
    "### Exponential Solution\n",
    "\n",
    "The Kolmogorov backward equation is a matrix-valued differential equation.\n",
    "\n",
    "Recall that, for a scalar differential equation $ y'_t = a y_t $ with constant\n",
    "$ a $ and initial condition $ y_0 $, the solution is $ y_t = e^{ta} y_0 $.\n",
    "\n",
    "This, along with $ P_0 = I $, encourages us to guess that the solution to\n",
    "Kolmogorov’s backward equation [(4.4)](#equation-kolbackeq) is\n",
    "\n",
    "\n",
    "<a id='equation-expsol'></a>\n",
    "$$\n",
    "P_t = e^{t Q} \\tag{4.5}\n",
    "$$\n",
    "\n",
    "where the right hand side is the [matrix exponential](https://en.wikipedia.org/wiki/Matrix_exponential), with definition\n",
    "\n",
    "\n",
    "<a id='equation-expofun'></a>\n",
    "$$\n",
    "e^{tQ} \n",
    "    = \\sum_{k \\geq 0} \\frac{1}{k!} (tQ)^k\n",
    "    = I + tQ + \\frac{t^2}{2!} Q^2 + \\cdots \\tag{4.6}\n",
    "$$\n",
    "\n",
    "Working element by element, it is not difficult to confirm that\n",
    "the derivative of the exponential function $ t \\mapsto e^{tQ} $ is\n",
    "\n",
    "\n",
    "<a id='equation-expoderiv'></a>\n",
    "$$\n",
    "\\frac{d}{dt} e^{t Q} = Q e^{t Q} = e^{t Q} Q \\tag{4.7}\n",
    "$$\n",
    "\n",
    "Hence, differentiating [(4.5)](#equation-expsol) gives $ P'_t = Q e^{t Q} = Q P_t $, which\n",
    "convinces us that the exponential solution satisfies [(4.4)](#equation-kolbackeq).\n",
    "\n",
    "Notice that our solution\n",
    "\n",
    "\n",
    "<a id='equation-psolq'></a>\n",
    "$$\n",
    "P_t = e^{t Q} \n",
    "    \\quad \\text{where } \\;\n",
    "    Q(x, y) := \\lambda(x) (K(x, y) - I(x, y)) \\tag{4.8}\n",
    "$$\n",
    "\n",
    "for the semigroup of the jump process $ (X_t) $ associated with the jump matrix\n",
    "$ K $ and the jump intensity function $ \\lambda \\colon S \\to (0, \\infty) $ is\n",
    "consistent with our earlier result.\n",
    "\n",
    "In particular, we [showed](markov_prop.ipynb#consjumptransemi) that, for the model with\n",
    "constant jump intensity $ \\lambda $, we have $ P_t = e^{t \\lambda (K - I)} $.\n",
    "\n",
    "This is obviously a special case of [(4.8)](#equation-psolq)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddad6783",
   "metadata": {},
   "source": [
    "## Properties of the Solution\n",
    "\n",
    "Let’s investigate further the properties of the exponential solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a03662f5",
   "metadata": {},
   "source": [
    "### Checking the Transition Semigroup Properties\n",
    "\n",
    "While we have confirmed that $ P_t = e^{t Q} $ solves the Kolmogorov backward\n",
    "equation, we still need to check that this solution is a Markov semigroup."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbf5868b",
   "metadata": {},
   "source": [
    "###  (From Jump Chain to Semigroup)\n",
    "\n",
    "Let $ \\lambda $ map $ S $ to $ \\RR_+ $ and let $ K $ be a Markov matrix on $ S $.\n",
    "If $ P_t = e^{t Q} $ for all $ t \\geq 0 $, where\n",
    "$ Q(x, y) = \\lambda(x) (K(x, y) - I(x, y)) $, then $ (P_t) $ is a Markov\n",
    "semigroup on $ S $.\n",
    "\n",
    "Proof. Observe first that $ Q $ has zero row sums, since\n",
    "\n",
    "$$\n",
    "\\sum_y Q(x, y) \n",
    "    = \\lambda(x) \\sum_y (K(x, y) - I(x, y))\n",
    "    = 0\n",
    "$$\n",
    "\n",
    "As a small exercise, you can check that, with $ 1 $\n",
    "representing a column vector of ones, the following is true\n",
    "\n",
    "\n",
    "<a id='equation-zrsnec'></a>\n",
    "$$\n",
    "Q \\text{ has zero row sums }\n",
    "    \\iff\n",
    "    Q^k 1 = 0 \\text{ for all } k \\geq 1 \\tag{4.9}\n",
    "$$\n",
    "\n",
    "This implies that $ Q^k 1 = 0 $ for all $ k $ and, as a result, for any $ t \\geq\n",
    "0 $,\n",
    "\n",
    "$$\n",
    "P_t 1\n",
    "    = e^{tQ} 1\n",
    "    = I1 + tQ1 + \\frac{t^2}{2!} Q^2 1 + \\cdots\n",
    "    = I1 = 1\n",
    "$$\n",
    "\n",
    "In other words, each $ P_t $ has unit row sums.\n",
    "\n",
    "Next we check nonnegativity of all elements of $ P_t $ (which can easily fail for\n",
    "matrix exponentials).\n",
    "\n",
    "To this end, adopting an argument from [[Stroock, 2013](zreferences.ipynb#id14)], we\n",
    "set $ m := \\max_x \\lambda(x) $ and $ \\hat P := I + Q / m $.\n",
    "\n",
    "It is not difficult to check that $ \\hat P $ is a Markov matrix and $ Q = m( \\hat\n",
    "P - I) $.\n",
    "\n",
    "Recalling that, for matrix exponentials, $ e^{A+B} = e^A e^B $ whenever $ AB =\n",
    "BA $, we have\n",
    "\n",
    "$$\n",
    "e^{tQ} \n",
    "    = e^{tm (\\hat P - I)} \n",
    "    = e^{-tm I} e^{tm \\hat P}\n",
    "    = e^{-tm} \n",
    "        \\left( \n",
    "            I + tm \\hat P + \\frac{(tm)^2}{2!} \\hat P^2 + \\cdots\n",
    "        \\right)\n",
    "$$\n",
    "\n",
    "It is clear from this representation that all entries of $ e^{tQ} $ are\n",
    "nonnegative.\n",
    "\n",
    "Finally, we need to check the continuity condition $ P_t(x, y) \\to I(x,y) $ as\n",
    "$ t \\to 0 $, which is also part of the definition of a Markov semigroup.\n",
    "This is immediate, in the present case, because the exponential function is\n",
    "continuous, and hence $ P_t = e^{tQ} \\to e^0 = I $.\n",
    "\n",
    "We can now be reassured that our solution to the Kolmogorov backward equation\n",
    "is indeed a Markov semigroup."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adbdc58f",
   "metadata": {},
   "source": [
    "### Uniqueness\n",
    "\n",
    "Might there be another, entirely different Markov semigroup that also\n",
    "satisfies the Kolmogorov backward equation?\n",
    "\n",
    "The answer is no:  linear ODEs in finite dimensional vector space with\n",
    "constant coefficients and fixed initial conditions (in this case $ P_0 = I $)\n",
    "have unique solutions.\n",
    "\n",
    "In fact it’s not hard to supply a proof — see the exercises."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7a8ea91",
   "metadata": {},
   "source": [
    "## Application: The Inventory Model\n",
    "\n",
    "Let us look at a modified version of the inventory model where jump\n",
    "intensities depend on the state.\n",
    "\n",
    "In particular, the wait time for new inventory will now be exponential at rate\n",
    "$ \\gamma $.\n",
    "\n",
    "The arrival rate for customers will still be denoted by $ \\lambda $ and allowed\n",
    "to differ from $ \\gamma $.\n",
    "\n",
    "For parameters we take"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82fb13c4",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "α = 0.6\n",
    "λ = 0.5\n",
    "γ = 0.1\n",
    "b = 10\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7c1e806",
   "metadata": {},
   "source": [
    "Our plan is to investigate the distribution $ \\psi_T $ of $ X_T $ at $ T=30 $.\n",
    "\n",
    "We will do this by simulating many independent draws of $ X_T $ and\n",
    "histogramming them.\n",
    "\n",
    "(In the exercises you are asked to calculate $ \\psi_T $ a different way, via\n",
    "[(4.8)](#equation-psolq).)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff5a0ff",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "@njit\n",
    "def draw_X(T, X_0, max_iter=5000):\n",
    "    \"\"\"\n",
    "    Generate one draw of X_T given X_0.\n",
    "    \"\"\"\n",
    "\n",
    "    J, Y = 0, X_0\n",
    "    m = 0\n",
    "\n",
    "    while m < max_iter:\n",
    "        s = 1/γ if Y == 0 else 1/λ\n",
    "        W = np.random.exponential(scale=s)  # W ~ E(λ)\n",
    "        J += W\n",
    "        if J >= T:\n",
    "            return Y\n",
    "        # Otherwise update Y\n",
    "        if Y == 0:\n",
    "            Y = b\n",
    "        else:\n",
    "            U = np.random.geometric(α)\n",
    "            Y = Y - min(Y, U)\n",
    "        m += 1\n",
    "\n",
    "\n",
    "@njit\n",
    "def independent_draws(T=10, num_draws=100):\n",
    "    \"Generate a vector of independent draws of X_T.\"\n",
    "\n",
    "    draws = np.empty(num_draws, dtype=np.int64)\n",
    "\n",
    "    for i in range(num_draws):\n",
    "        X_0 = np.random.binomial(b+1, 0.25)\n",
    "        draws[i] = draw_X(T, X_0)\n",
    "\n",
    "    return draws\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa768099",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "T = 30\n",
    "n = b + 1\n",
    "draws = independent_draws(T, num_draws=100_000)\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.bar(range(n), [np.mean(draws == i) for i in range(n)], width=0.8, alpha=0.6)\n",
    "ax.set_xlabel(\"inventory\", fontsize=14)\n",
    "\n",
    "plt.show()\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2858e03f",
   "metadata": {},
   "source": [
    "If you experiment with the code above, you will see that the large amount of\n",
    "mass on zero is due to the low arrival rate $ \\gamma $ for inventory."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11c25de8",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a08f2d45",
   "metadata": {},
   "source": [
    "## Exercise 4.1\n",
    "\n",
    "In the discussion above, we generated an approximation of $ \\psi_T $ when\n",
    "$ T=30 $, the initial condition is Binomial$ (n, 0.25) $ and parameters\n",
    "are set to"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d258f582",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "α = 0.6\n",
    "λ = 0.5\n",
    "γ = 0.1\n",
    "b = 10\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4343b166",
   "metadata": {},
   "source": [
    "The calculation was done by simulating independent draws and histogramming.\n",
    "\n",
    "Try to generate the same figure using [(4.8)](#equation-psolq) instead, modifying code from\n",
    "[our lecture](markov_prop.ipynb) on the Markov property."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bd274b9",
   "metadata": {},
   "source": [
    "## Exercise 4.2\n",
    "\n",
    "Prove that differentiating [(4.1)](#equation-kbinteg) at each $ (x, y) $ yields [(4.4)](#equation-kolbackeq)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2699d11f",
   "metadata": {},
   "source": [
    "## Exercise 4.3\n",
    "\n",
    "We claimed above that the solution $ P_t = e^{t Q} $ is the unique\n",
    "Markov semigroup satisfying the backward equation $ P'_t = Q P_t $.\n",
    "\n",
    "Try to supply a proof.\n",
    "\n",
    "(This is not an easy exercise but worth thinking about in any case.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d49981e",
   "metadata": {},
   "source": [
    "## Solutions\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">code is currently not supported in `sphinx-exercise`\n",
    "so code-cell solutions are immediately after this\n",
    "solution block."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ac7869d",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 4.1](#kolmogorov-bwd-1)\n",
    "\n",
    "Here is one solution:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04cd6292",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```ipython3\n",
    "α = 0.6\n",
    "λ = 0.5\n",
    "γ = 0.1\n",
    "b = 10\n",
    "\n",
    "states = np.arange(n)\n",
    "I = np.identity(n)\n",
    "\n",
    "# Embedded jump chain matrix\n",
    "K = np.zeros((n, n))\n",
    "K[0, -1] = 1\n",
    "for i in range(1, n):\n",
    "    for j in range(0, i):\n",
    "        if j == 0:\n",
    "            K[i, j] = (1 - α)**(i-1)\n",
    "        else:\n",
    "            K[i, j] = α * (1 - α)**(i-j-1)\n",
    "\n",
    "# Jump intensities as a function of the state\n",
    "r = np.ones(n) * λ\n",
    "r[0] = γ\n",
    "\n",
    "# Q matrix\n",
    "Q = np.empty_like(K)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        Q[i, j] = r[i] * (K[i, j] - I[i, j])\n",
    "\n",
    "def P_t(ψ, t):\n",
    "    return ψ @ expm(t * Q)\n",
    "\n",
    "ψ_0 = binom.pmf(states, n, 0.25)\n",
    "ψ_T = P_t(ψ_0, T)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.bar(range(n), ψ_T, width=0.8, alpha=0.6)\n",
    "ax.set_xlabel(\"inventory\", fontsize=14)\n",
    "\n",
    "plt.show()\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cefe760a",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 4.2](#kolmogorov-bwd-2)\n",
    "\n",
    "One can easily verify that, when $ f $ is a differentiable function and $ \\alpha >\n",
    "0 $, we have\n",
    "\n",
    "\n",
    "<a id='equation-gdiff'></a>\n",
    "$$\n",
    "g(t) = e^{- t \\alpha} f(t)\n",
    "    \\quad \\implies \\quad\n",
    "    g'(t) = e^{- t \\alpha} f'(t) - \\alpha g(t) \\tag{4.10}\n",
    "$$\n",
    "\n",
    "Note also that, with the change of variable $ s = t - \\tau $, we can rewrite\n",
    "[(4.1)](#equation-kbinteg) as\n",
    "\n",
    "\n",
    "<a id='equation-kbinteg2'></a>\n",
    "$$\n",
    "P_t(x, y) =\n",
    "    e^{-t \\lambda(x)}\n",
    "    \\left\\{\n",
    "        I(x, y)\n",
    "        + \\lambda(x)\n",
    "        \\int_0^t (K P_s)(x, y) e^{s \\lambda(x)} d s\n",
    "    \\right\\} \\tag{4.11}\n",
    "$$\n",
    "\n",
    "Applying [(4.10)](#equation-gdiff) yields\n",
    "\n",
    "$$\n",
    "P'_t(x, y)\n",
    "    = e^{-t \\lambda(x)}\n",
    "        \\left\\{ \n",
    "             \\lambda(x)\n",
    "             (K P_t)(x, y) e^{t \\lambda(x)}\n",
    "        \\right\\}\n",
    "        - \\lambda(x) P_t(x, y)\n",
    "$$\n",
    "\n",
    "After minor rearrangements this becomes\n",
    "\n",
    "$$\n",
    "P'_t(x, y)\n",
    "    = \\lambda(x) [ (K - I)  P_t](x, y)\n",
    "$$\n",
    "\n",
    "which is identical to [(4.4)](#equation-kolbackeq)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73c6e51a",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 4.3](#kolmogorov-bwd-3)\n",
    "\n",
    "Here is one proof of uniqueness.\n",
    "\n",
    "Suppose that $ (\\hat P_t) $ is another Markov semigroup satisfying\n",
    "$ P'_t = Q P_t $.\n",
    "\n",
    "Fix $ t > 0 $ and let $ V_s $ be defined by $ V_s = P_s \\hat P_{t-s} $ for all $ s\n",
    "\\geq 0 $.\n",
    "\n",
    "Note that $ V_0 = \\hat P_t $ and $ V_t = P_t $.\n",
    "\n",
    "Note also that $ s \\mapsto V_s $ is differentiable, with derivative\n",
    "\n",
    "$$\n",
    "V'_s \n",
    "    = P'_s \\hat P_{t-s} - P_s \\hat P'_{t-s}\n",
    "    = P_s Q \\hat P_{t-s} - P_s Q \\hat P_{t-s}\n",
    "    = 0\n",
    "$$\n",
    "\n",
    "where, in the second last equality, we used [(4.7)](#equation-expoderiv).\n",
    "\n",
    "Hence $ V_s $ is constant, so our previous observations $ V_0 = \\hat P_t $ and $ V_t = P_t $\n",
    "now yield $ \\hat P_t = P_t $.\n",
    "\n",
    "Since $ t $ was arbitrary, the proof is now done."
   ]
  }
 ],
 "metadata": {
  "date": 1732753873.910782,
  "filename": "kolmogorov_bwd.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "The Kolmogorov Backward Equation"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}