markov_chains_I.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e0ea056b",
   "metadata": {},
   "source": [
    "# Markov Chains: Basic Concepts\n",
    "\n",
    "\n",
    "<a id='index-0'></a>\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3a0cedd",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54d5b59c",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Markov chains provide  a way to model situations in which  the past casts shadows on the future.\n",
    "\n",
    "By this we mean that observing measurements about a present situation can help us forecast future situations.\n",
    "\n",
    "This can be possible when there are statistical dependencies among measurements of something taken at different points of time.\n",
    "\n",
    "For example,\n",
    "\n",
    "- inflation next year might co-vary  with inflation this year  \n",
    "- unemployment next month might co-vary with unemployment this month  \n",
    "\n",
    "\n",
    "Markov chains are a workhorse for economics and finance.\n",
    "\n",
    "The theory of Markov chains is beautiful and provides many insights into\n",
    "probability and dynamics.\n",
    "\n",
    "In this  lecture, we will\n",
    "\n",
    "- review some of the key ideas from the theory of Markov chains and  \n",
    "- show how Markov chains appear in some economic applications.  \n",
    "\n",
    "\n",
    "Let’s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a8b87e4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import quantecon as qe\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "from matplotlib import cm\n",
    "import matplotlib as mpl\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "from matplotlib.patches import Polygon\n",
    "from mpl_toolkits.mplot3d.art3d import Poly3DCollection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f1b0d89",
   "metadata": {},
   "source": [
    "## Definitions and examples\n",
    "\n",
    "In this section we provide some definitions and  elementary examples.\n",
    "\n",
    "\n",
    "<a id='finite-dp-stoch-mat'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c1610c",
   "metadata": {},
   "source": [
    "### Stochastic matrices\n",
    "\n",
    "Recall that a **probability mass function** over $ n $ possible outcomes is a\n",
    "nonnegative $ n $-vector $ p $ that sums to one.\n",
    "\n",
    "For example, $ p = (0.2, 0.2, 0.6) $ is a probability mass function over $ 3 $ outcomes.\n",
    "\n",
    "A **stochastic matrix** (or **Markov matrix**)  is an $ n \\times n $ square matrix $ P $\n",
    "such that each row of $ P $ is a probability mass function over $ n $ outcomes.\n",
    "\n",
    "In other words,\n",
    "\n",
    "1. each element of $ P $ is nonnegative, and  \n",
    "1. each row of $ P $ sums to one  \n",
    "\n",
    "\n",
    "If $ P $ is a stochastic matrix, then so is the $ k $-th power $ P^k $ for all $ k \\in \\mathbb N $.\n",
    "\n",
    "You are asked to check this in [an exercise](#mc1_ex_3) below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9e9affe",
   "metadata": {},
   "source": [
    "### Markov chains\n",
    "\n",
    "Now we can introduce Markov chains.\n",
    "\n",
    "Before defining a Markov chain rigorously, we’ll  give some examples.\n",
    "\n",
    "\n",
    "<a id='mc-eg2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf4f1532",
   "metadata": {},
   "source": [
    "#### Example 1\n",
    "\n",
    "From  US unemployment data, Hamilton [[Hamilton, 2005](https://intro.quantecon.org/zreferences.html#id187)] estimated the following dynamics.\n",
    "\n",
    "![https://intro.quantecon.org/_static/lecture_specific/markov_chains_I/Hamilton.png](https://intro.quantecon.org/_static/lecture_specific/markov_chains_I/Hamilton.png)\n",
    "\n",
    "Here there are three **states**\n",
    "\n",
    "- “ng” represents normal growth  \n",
    "- “mr” represents mild recession  \n",
    "- “sr” represents severe recession  \n",
    "\n",
    "\n",
    "The arrows represent transition probabilities over one month.\n",
    "\n",
    "For example, the arrow from mild recession to normal growth has 0.145 next to it.\n",
    "\n",
    "This tells us that, according to past data, there is a 14.5% probability of transitioning from mild recession to normal growth in one month.\n",
    "\n",
    "The arrow from normal growth back to normal growth tells us that there is a\n",
    "97% probability of transitioning from normal growth to normal growth (staying\n",
    "in the same state).\n",
    "\n",
    "Note that these are conditional probabilities — the probability of\n",
    "transitioning from one state to another (or staying at the same one) conditional on the\n",
    "current state.\n",
    "\n",
    "To make the problem easier to work with numerically, let’s convert states to\n",
    "numbers.\n",
    "\n",
    "In particular, we agree that\n",
    "\n",
    "- state 0 represents normal growth  \n",
    "- state 1 represents mild recession  \n",
    "- state 2 represents severe recession  \n",
    "\n",
    "\n",
    "Let $ X_t $ record the value of the state at time $ t $.\n",
    "\n",
    "Now we can write the statement “there is a 14.5% probability of transitioning from mild recession to normal growth in one month” as\n",
    "\n",
    "$$\n",
    "\\mathbb P\\{X_{t+1} = 0 \\,|\\, X_t = 1\\} = 0.145\n",
    "$$\n",
    "\n",
    "We can collect all of these conditional probabilities into a matrix, as follows\n",
    "\n",
    "$$\n",
    "P =\n",
    "\\begin{bmatrix}\n",
    "0.971 & 0.029 & 0 \\\\\n",
    "0.145 & 0.778 & 0.077 \\\\\n",
    "0 & 0.508 & 0.492\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notice that $ P $ is a stochastic matrix.\n",
    "\n",
    "Now we have the following relationship\n",
    "\n",
    "$$\n",
    "P(i,j)\n",
    "    = \\mathbb P\\{X_{t+1} = j \\,|\\, X_t = i\\}\n",
    "$$\n",
    "\n",
    "This holds for any $ i,j $ between 0 and 2.\n",
    "\n",
    "In particular, $ P(i,j) $ is the\n",
    "probability of transitioning from state $ i $ to state $ j $ in one month.\n",
    "\n",
    "\n",
    "<a id='mc-eg1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "758fddfd",
   "metadata": {},
   "source": [
    "#### Example 2\n",
    "\n",
    "Consider a worker who, at any given time $ t $, is either unemployed (state 0)\n",
    "or employed (state 1).\n",
    "\n",
    "Suppose that, over a one-month period,\n",
    "\n",
    "1. the unemployed worker finds a job with probability $ \\alpha \\in (0, 1) $.  \n",
    "1. the employed worker loses her job and becomes unemployed with probability $ \\beta \\in (0, 1) $.  \n",
    "\n",
    "\n",
    "Given the above information, we can write out the transition probabilities in matrix form as\n",
    "\n",
    "\n",
    "<a id='equation-p-unempemp'></a>\n",
    "$$\n",
    "P =\n",
    "\\begin{bmatrix}\n",
    "    1 - \\alpha & \\alpha \\\\\n",
    "    \\beta & 1 - \\beta\n",
    "\\end{bmatrix} \\tag{33.1}\n",
    "$$\n",
    "\n",
    "For example,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    P(0,1)\n",
    "        & =\n",
    "        \\text{ probability of transitioning from state \\$0\\$ to state \\$1\\$ in one month}\n",
    "        \\\\\n",
    "        & =\n",
    "        \\text{ probability finding a job next month}\n",
    "        \\\\\n",
    "        & = \\alpha\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Suppose we can estimate the values $ \\alpha $ and $ \\beta $.\n",
    "\n",
    "Then we can address a range of questions, such as\n",
    "\n",
    "- What is the average duration of unemployment?  \n",
    "- Over the long-run, what fraction of the time does a worker find herself unemployed?  \n",
    "- Conditional on employment, what is the probability of becoming unemployed at least once over the next 12 months?  \n",
    "\n",
    "\n",
    "We’ll cover some of these applications below.\n",
    "\n",
    "\n",
    "<a id='mc-eg3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05c23b7c",
   "metadata": {},
   "source": [
    "#### Example 3\n",
    "\n",
    "Imam and Temple [[Imam and Temple, 2023](https://intro.quantecon.org/zreferences.html#id289)] categorize political institutions into\n",
    "three types: democracy $ \\text{(D)} $, autocracy $ \\text{(A)} $, and an intermediate\n",
    "state called anocracy $ \\text{(N)} $.\n",
    "\n",
    "Each institution can have two potential development regimes: collapse $ \\text{(C)} $ and growth $ \\text{(G)} $. This results in six possible states: $ \\text{DG, DC, NG, NC, AG} $ and $ \\text{AC} $.\n",
    "\n",
    "Imam and Temple [[Imam and Temple, 2023](https://intro.quantecon.org/zreferences.html#id289)] estimate the following transition\n",
    "probabilities:\n",
    "\n",
    "$$\n",
    "P :=\n",
    "\\begin{bmatrix}\n",
    "0.86 & 0.11 & 0.03 & 0.00 & 0.00 & 0.00 \\\\\n",
    "0.52 & 0.33 & 0.13 & 0.02 & 0.00 & 0.00 \\\\\n",
    "0.12 & 0.03 & 0.70 & 0.11 & 0.03 & 0.01 \\\\\n",
    "0.13 & 0.02 & 0.35 & 0.36 & 0.10 & 0.04 \\\\\n",
    "0.00 & 0.00 & 0.09 & 0.11 & 0.55 & 0.25 \\\\\n",
    "0.00 & 0.00 & 0.09 & 0.15 & 0.26 & 0.50\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60668461",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "nodes = ['DG', 'DC', 'NG', 'NC', 'AG', 'AC']\n",
    "P = [[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],\n",
    "     [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],\n",
    "     [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],\n",
    "     [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],\n",
    "     [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],\n",
    "     [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f046c51",
   "metadata": {},
   "source": [
    "Here is a visualization, with darker colors indicating higher probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1f8e9a5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "G = nx.MultiDiGraph()\n",
    "\n",
    "for start_idx, node_start in enumerate(nodes):\n",
    "    for end_idx, node_end in enumerate(nodes):\n",
    "        value = P[start_idx][end_idx]\n",
    "        if value != 0:\n",
    "            G.add_edge(node_start,node_end, weight=value)\n",
    "\n",
    "pos = nx.spring_layout(G, seed=10)\n",
    "fig, ax = plt.subplots()\n",
    "nx.draw_networkx_nodes(G, pos, node_size=600, edgecolors='black', node_color='white')\n",
    "nx.draw_networkx_labels(G, pos)\n",
    "\n",
    "arc_rad = 0.2\n",
    "\n",
    "edges = nx.draw_networkx_edges(G, pos, ax=ax, connectionstyle=f'arc3, rad = {arc_rad}', edge_cmap=cm.Blues, width=2,\n",
    "    edge_color=[G[nodes[0]][nodes[1]][0]['weight'] for nodes in G.edges])\n",
    "\n",
    "pc = mpl.collections.PatchCollection(edges, cmap=cm.Blues)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_axis_off()\n",
    "plt.colorbar(pc, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7a86d81",
   "metadata": {},
   "source": [
    "Looking at the data, we see that democracies tend to have longer-lasting growth\n",
    "regimes compared to autocracies (as indicated by the lower probability of\n",
    "transitioning from growth to growth in autocracies).\n",
    "\n",
    "We can also find a higher probability from collapse to growth in democratic regimes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b640ee5",
   "metadata": {},
   "source": [
    "### Defining Markov chains\n",
    "\n",
    "So far we’ve given examples of Markov chains but we haven’t defined them.\n",
    "\n",
    "Let’s do that now.\n",
    "\n",
    "To begin, let $ S $ be a finite set $ \\{x_1, \\ldots, x_n\\} $ with $ n $ elements.\n",
    "\n",
    "The set $ S $ is called the **state space** and $ x_1, \\ldots, x_n $ are the **state values**.\n",
    "\n",
    "A **distribution** $ \\psi $ on $ S $ is a probability mass function of length $ n $, where $ \\psi(i) $ is the amount of probability allocated to state $ x_i $.\n",
    "\n",
    "A **Markov chain** $ \\{X_t\\} $ on $ S $ is a sequence of random variables taking values in $ S $\n",
    "that have the **Markov property**.\n",
    "\n",
    "This means that, for any date $ t $ and any state $ y \\in S $,\n",
    "\n",
    "\n",
    "<a id='equation-fin-markov-mp'></a>\n",
    "$$\n",
    "\\mathbb P \\{ X_{t+1} = y  \\,|\\, X_t \\}\n",
    "= \\mathbb P \\{ X_{t+1}  = y \\,|\\, X_t, X_{t-1}, \\ldots \\} \\tag{33.2}\n",
    "$$\n",
    "\n",
    "This means that once we know the current state $ X_t $,  adding knowledge of earlier states $ X_{t-1}, X_{t-2} $ provides no additional information about probabilities of *future* states.\n",
    "\n",
    "Thus, the dynamics of a Markov chain are fully determined by the set of **conditional probabilities**\n",
    "\n",
    "\n",
    "<a id='equation-mpp'></a>\n",
    "$$\n",
    "P(x, y) := \\mathbb P \\{ X_{t+1} = y \\,|\\, X_t = x \\}\n",
    "\\qquad (x, y \\in S) \\tag{33.3}\n",
    "$$\n",
    "\n",
    "By construction,\n",
    "\n",
    "- $ P(x, y) $ is the probability of going from $ x $ to $ y $ in one unit of time (one step)  \n",
    "- $ P(x, \\cdot) $ is the conditional distribution of $ X_{t+1} $ given $ X_t = x $  \n",
    "\n",
    "\n",
    "We can view $ P $ as a stochastic matrix where\n",
    "\n",
    "$$\n",
    "P_{ij} = P(x_i, x_j)\n",
    "    \\qquad 1 \\leq i, j \\leq n\n",
    "$$\n",
    "\n",
    "Going the other way, if we take a stochastic matrix $ P $, we can generate a Markov\n",
    "chain $ \\{X_t\\} $ as follows:\n",
    "\n",
    "- draw $ X_0 $ from a distribution $ \\psi_0 $ on $ S $  \n",
    "- for each $ t = 0, 1, \\ldots $, draw $ X_{t+1} $ from $ P(X_t,\\cdot) $  \n",
    "\n",
    "\n",
    "By construction, the resulting process satisfies [(33.3)](#equation-mpp)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7aad296f",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "A good way to study Markov chains is to simulate them.\n",
    "\n",
    "Let’s start by doing this ourselves and then look at libraries that can help\n",
    "us.\n",
    "\n",
    "In these exercises, we’ll take the state space to be $ S = 0,\\ldots, n-1 $.\n",
    "\n",
    "(We start at $ 0 $ because Python arrays are indexed from $ 0 $.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fafd940a",
   "metadata": {},
   "source": [
    "### Writing our own simulation code\n",
    "\n",
    "To simulate a Markov chain, we need\n",
    "\n",
    "1. a stochastic matrix $ P $ and  \n",
    "1. a probability mass function $ \\psi_0 $ of length $ n $ from which to draw an initial realization of $ X_0 $.  \n",
    "\n",
    "\n",
    "The Markov chain is then constructed as follows:\n",
    "\n",
    "1. At time $ t=0 $, draw a realization of $ X_0 $ from the distribution $ \\psi_0 $.  \n",
    "1. At each subsequent time $ t $, draw a realization of the new state $ X_{t+1} $ from $ P(X_t, \\cdot) $.  \n",
    "\n",
    "\n",
    "(That is, draw from row $ X_t $ of $ P $.)\n",
    "\n",
    "To implement this simulation procedure, we need a method for generating draws\n",
    "from a discrete distribution.\n",
    "\n",
    "For this task, we’ll use `random.draw` from [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
    "\n",
    "To use `random.draw`, we first need to convert the probability mass function\n",
    "to a cumulative distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2f7d74c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ_0 = (0.3, 0.7)           # probabilities over {0, 1}\n",
    "cdf = np.cumsum(ψ_0)       # convert into cumulative distribution\n",
    "qe.random.draw(cdf, 5)   # generate 5 independent draws from ψ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51801458",
   "metadata": {},
   "source": [
    "We’ll write our code as a function that accepts the following three arguments\n",
    "\n",
    "- A stochastic matrix `P`.  \n",
    "- An initial distribution `ψ_0`.  \n",
    "- A positive integer `ts_length` representing the length of the time series the function should return.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d826eaa",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def mc_sample_path(P, ψ_0=None, ts_length=1_000):\n",
    "\n",
    "    # set up\n",
    "    P = np.asarray(P)\n",
    "    X = np.empty(ts_length, dtype=int)\n",
    "\n",
    "    # Convert each row of P into a cdf\n",
    "    P_dist = np.cumsum(P, axis=1)  # Convert rows into cdfs\n",
    "\n",
    "    # draw initial state, defaulting to 0\n",
    "    if ψ_0 is not None:\n",
    "        X_0 = qe.random.draw(np.cumsum(ψ_0))\n",
    "    else:\n",
    "        X_0 = 0\n",
    "\n",
    "    # simulate\n",
    "    X[0] = X_0\n",
    "    for t in range(ts_length - 1):\n",
    "        X[t+1] = qe.random.draw(P_dist[X[t], :])\n",
    "\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42dc79ac",
   "metadata": {},
   "source": [
    "Let’s see how it works using the small matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4164687",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[0.4, 0.6],\n",
    "     [0.2, 0.8]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "090d5d2a",
   "metadata": {},
   "source": [
    "Here’s a short time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d90f6fb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc_sample_path(P, ψ_0=(1.0, 0.0), ts_length=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4532e9c3",
   "metadata": {},
   "source": [
    "It can be shown that for a long series drawn from `P`, the fraction of the\n",
    "sample that takes value 0 will be about 0.25.\n",
    "\n",
    "(We will explain why [later](https://intro.quantecon.org/markov_chains_II.html#ergodicity).)\n",
    "\n",
    "Moreover, this is true regardless of the initial distribution from which\n",
    "$ X_0 $ is drawn.\n",
    "\n",
    "The following code illustrates this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35fc0aca",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = mc_sample_path(P, ψ_0=(0.1, 0.9), ts_length=1_000_000)\n",
    "np.mean(X == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8f6dc2a",
   "metadata": {},
   "source": [
    "You can try changing the initial distribution to confirm that the output is\n",
    "always close to 0.25 (for the `P` matrix above)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4683e533",
   "metadata": {},
   "source": [
    "### Using QuantEcon’s routines\n",
    "\n",
    "[QuantEcon.py](http://quantecon.org/quantecon-py) has routines for handling Markov chains, including simulation.\n",
    "\n",
    "Here’s an illustration using the same $ P $ as the preceding example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb6a13fb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc = qe.MarkovChain(P)\n",
    "X = mc.simulate(ts_length=1_000_000)\n",
    "np.mean(X == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d62caf3b",
   "metadata": {},
   "source": [
    "The `simulate` routine is faster (because it is [JIT compiled](https://python-programming.quantecon.org/numba.html#numba-link))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed1fe454",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%time mc_sample_path(P, ts_length=1_000_000) # Our homemade code version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "658c538a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%time mc.simulate(ts_length=1_000_000) # qe code version"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb100248",
   "metadata": {},
   "source": [
    "#### Adding state values and initial conditions\n",
    "\n",
    "If we wish to, we can provide a specification of state values to `MarkovChain`.\n",
    "\n",
    "These state values can be integers, floats, or even strings.\n",
    "\n",
    "The following code illustrates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53df5b51",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc = qe.MarkovChain(P, state_values=('unemployed', 'employed'))\n",
    "mc.simulate(ts_length=4, init='employed')  # Start at employed initial state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e084e3b0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.simulate(ts_length=4, init='unemployed')  # Start at unemployed initial state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5067e7b5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.simulate(ts_length=4)  # Start at randomly chosen initial state"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ee0e209",
   "metadata": {},
   "source": [
    "If we want to see indices rather than state values as outputs as  we can use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "467c4520",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.simulate_indices(ts_length=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df8edc9",
   "metadata": {},
   "source": [
    "\n",
    "<a id='mc-md'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22afd8a8",
   "metadata": {},
   "source": [
    "## Distributions over time\n",
    "\n",
    "We learned that\n",
    "\n",
    "1. $ \\{X_t\\} $ is a Markov chain with stochastic matrix $ P $  \n",
    "1. the distribution of $ X_t $ is known to be $ \\psi_t $  \n",
    "\n",
    "\n",
    "What then is the distribution of $ X_{t+1} $, or, more generally, of $ X_{t+m} $?\n",
    "\n",
    "To answer this, we let $ \\psi_t $ be the distribution of $ X_t $ for $ t = 0, 1, 2, \\ldots $.\n",
    "\n",
    "Our first aim is to find $ \\psi_{t + 1} $ given $ \\psi_t $ and $ P $.\n",
    "\n",
    "To begin, pick any $ y \\in S $.\n",
    "\n",
    "To get the probability of being at $ y $ tomorrow (at $ t+1 $), we account for\n",
    "all ways this can happen and sum their probabilities.\n",
    "\n",
    "This leads to\n",
    "\n",
    "$$\n",
    "\\mathbb P \\{X_{t+1} = y \\}\n",
    "   = \\sum_{x \\in S} \\mathbb P \\{ X_{t+1} = y \\, | \\, X_t = x \\}\n",
    "               \\cdot \\mathbb P \\{ X_t = x \\}\n",
    "$$\n",
    "\n",
    "(We are using the [law of total probability](https://en.wikipedia.org/wiki/Law_of_total_probability).)\n",
    "\n",
    "Rewriting this statement in terms of  marginal and conditional probabilities gives\n",
    "\n",
    "$$\n",
    "\\psi_{t+1}(y) = \\sum_{x \\in S} P(x,y) \\psi_t(x)\n",
    "$$\n",
    "\n",
    "There are $ n $ such equations, one for each $ y \\in S $.\n",
    "\n",
    "If we think of $ \\psi_{t+1} $ and $ \\psi_t $ as row vectors, these $ n $ equations are summarized by the matrix expression\n",
    "\n",
    "\n",
    "<a id='equation-fin-mc-fr'></a>\n",
    "$$\n",
    "\\psi_{t+1} = \\psi_t P \\tag{33.4}\n",
    "$$\n",
    "\n",
    "Thus, we postmultiply by $ P $ to move a distribution forward one unit of time.\n",
    "\n",
    "By postmultiplying $ m $ times, we move a distribution forward $ m $ steps into the future.\n",
    "\n",
    "Hence, iterating on [(33.4)](#equation-fin-mc-fr), the expression $ \\psi_{t+m} = \\psi_t P^m $ is also valid — here $ P^m $ is the $ m $-th power of $ P $.\n",
    "\n",
    "As a special case, we see that if $ \\psi_0 $ is the initial distribution from\n",
    "which $ X_0 $ is drawn, then $ \\psi_0 P^m $ is the distribution of\n",
    "$ X_m $.\n",
    "\n",
    "This is very important, so let’s repeat it\n",
    "\n",
    "\n",
    "<a id='equation-mdfmc'></a>\n",
    "$$\n",
    "X_0 \\sim \\psi_0 \\quad \\implies \\quad X_m \\sim \\psi_0 P^m \\tag{33.5}\n",
    "$$\n",
    "\n",
    "The general rule is that postmultiplying a distribution by $ P^m $ shifts it forward $ m $ units of time.\n",
    "\n",
    "Hence the following is also valid.\n",
    "\n",
    "\n",
    "<a id='equation-mdfmc2'></a>\n",
    "$$\n",
    "X_t \\sim \\psi_t \\quad \\implies \\quad X_{t+m} \\sim \\psi_t P^m \\tag{33.6}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='finite-mc-mstp'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6e6e356",
   "metadata": {},
   "source": [
    "### Multiple step transition probabilities\n",
    "\n",
    "We know that the probability of transitioning from $ x $ to $ y $ in\n",
    "one step is $ P(x,y) $.\n",
    "\n",
    "It turns out that the probability of transitioning from $ x $ to $ y $ in\n",
    "$ m $ steps is $ P^m(x,y) $, the $ (x,y) $-th element of the\n",
    "$ m $-th power of $ P $.\n",
    "\n",
    "To see why, consider again [(33.6)](#equation-mdfmc2), but now with a $ \\psi_t $ that puts all probability on state $ x $.\n",
    "\n",
    "Then $ \\psi_t $ is a vector with $ 1 $ in position $ x $ and zero elsewhere.\n",
    "\n",
    "Inserting this into [(33.6)](#equation-mdfmc2), we see that, conditional on $ X_t = x $, the distribution of $ X_{t+m} $ is the $ x $-th row of $ P^m $.\n",
    "\n",
    "In particular\n",
    "\n",
    "$$\n",
    "\\mathbb P \\{X_{t+m} = y \\,|\\, X_t = x \\} = P^m(x, y) = (x, y) \\text{-th element of } P^m\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5624b44e",
   "metadata": {},
   "source": [
    "### Example: probability of recession\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "Recall the stochastic matrix $ P $ for recession and growth [considered above](#mc-eg2).\n",
    "\n",
    "Suppose that the current state is unknown — perhaps statistics are available only at the *end* of the current month.\n",
    "\n",
    "We guess that the probability that the economy is in state $ x $ is $ \\psi_t(x) $ at time t.\n",
    "\n",
    "The probability of being in recession (either mild or severe) in 6 months time is given by\n",
    "\n",
    "$$\n",
    "(\\psi_t P^6)(1) + (\\psi_t P^6)(2)\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='mc-eg1-1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd914ce8",
   "metadata": {},
   "source": [
    "### Example 2: cross-sectional distributions\n",
    "\n",
    "The distributions we have been studying can be viewed either\n",
    "\n",
    "1. as probabilities or  \n",
    "1. as cross-sectional frequencies that the law of large numbers leads us to anticipate for large samples.  \n",
    "\n",
    "\n",
    "To illustrate, recall our model of employment/unemployment dynamics for a given worker [discussed above](#mc-eg1).\n",
    "\n",
    "Consider a large population of workers, each of whose lifetime experience is\n",
    "described by the specified dynamics, with each worker’s outcomes being\n",
    "realizations of processes that are statistically independent of all other\n",
    "workers’ processes.\n",
    "\n",
    "Let $ \\psi_t $ be the current *cross-sectional* distribution over $ \\{ 0, 1 \\} $.\n",
    "\n",
    "The cross-sectional distribution records fractions of workers employed and unemployed at a given moment $ t $.\n",
    "\n",
    "- For example, $ \\psi_t(0) $ is the unemployment rate at time $ t $.  \n",
    "\n",
    "\n",
    "What will the cross-sectional distribution be in 10 periods hence?\n",
    "\n",
    "The answer is $ \\psi_t P^{10} $, where $ P $ is the stochastic matrix in\n",
    "[(33.1)](#equation-p-unempemp).\n",
    "\n",
    "This is because each worker’s state evolves according to $ P $, so\n",
    "$ \\psi_t P^{10} $ is a [marginal distribution](https://en.wikipedia.org/wiki/Marginal_distribution)  for a single randomly selected\n",
    "worker.\n",
    "\n",
    "But when the sample is large, outcomes and probabilities are roughly equal (by an application of the law\n",
    "of large numbers).\n",
    "\n",
    "So for a very large (tending to infinite) population,\n",
    "$ \\psi_t P^{10} $ also represents  fractions of workers in\n",
    "each state.\n",
    "\n",
    "This is exactly the cross-sectional distribution.\n",
    "\n",
    "\n",
    "<a id='stationary'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5e5f467",
   "metadata": {},
   "source": [
    "## Stationary distributions\n",
    "\n",
    "As seen in [(33.4)](#equation-fin-mc-fr), we can shift a distribution forward one\n",
    "unit of time via postmultiplication by $ P $.\n",
    "\n",
    "Some distributions are invariant under this updating process — for example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6db2dfd0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = np.array([[0.4, 0.6],\n",
    "              [0.2, 0.8]])\n",
    "ψ = (0.25, 0.75)\n",
    "ψ @ P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04df1252",
   "metadata": {},
   "source": [
    "Notice that `ψ @ P` is the same as `ψ`.\n",
    "\n",
    "Such distributions are called **stationary** or **invariant**.\n",
    "\n",
    "\n",
    "<a id='mc-stat-dd'></a>\n",
    "Formally, a distribution $ \\psi^* $ on $ S $ is called **stationary** for $ P $ if $ \\psi^* P = \\psi^* $.\n",
    "\n",
    "Notice that, postmultiplying by $ P $, we have $ \\psi^* P^2 = \\psi^* P = \\psi^* $.\n",
    "\n",
    "Continuing in the same way leads to $ \\psi^* = \\psi^* P^t $ for all $ t \\ge 0 $.\n",
    "\n",
    "This tells us an important fact: If the distribution of $ \\psi_0 $ is a stationary distribution, then $ \\psi_t $ will have this same distribution for all $ t \\ge 0 $.\n",
    "\n",
    "The following theorem is proved in Chapter 4 of [[Sargent and Stachurski, 2023](https://intro.quantecon.org/zreferences.html#id24)] and numerous other sources."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f93b5ee2",
   "metadata": {},
   "source": [
    "## \n",
    "\n",
    "Every stochastic matrix $ P $ has at least one stationary distribution.\n",
    "\n",
    "Note that there can be many stationary distributions corresponding to a given\n",
    "stochastic matrix $ P $.\n",
    "\n",
    "- For example, if $ P $ is the identity matrix, then all distributions on $ S $ are stationary.  \n",
    "\n",
    "\n",
    "To get uniqueness, we need the Markov chain to “mix around,” so that the state\n",
    "doesn’t get stuck in some part of the state space.\n",
    "\n",
    "This gives some intuition for the following theorem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c3846a",
   "metadata": {},
   "source": [
    "## \n",
    "\n",
    "If $ P $ is everywhere positive, then $ P $ has exactly one stationary\n",
    "distribution.\n",
    "\n",
    "We will come back to this when we introduce irreducibility in the [next lecture](https://intro.quantecon.org/markov_chains_II.html) on Markov chains."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e68823ea",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Recall our model of the employment/unemployment dynamics of a particular worker [discussed above](#mc-eg1).\n",
    "\n",
    "If $ \\alpha \\in (0,1) $ and $ \\beta \\in (0,1) $, then the transition matrix is everywhere positive.\n",
    "\n",
    "Let $ \\psi^* = (p, 1-p) $ be the stationary distribution, so that $ p $\n",
    "corresponds to unemployment (state 0).\n",
    "\n",
    "Using $ \\psi^* = \\psi^* P $ and a bit of algebra yields\n",
    "\n",
    "$$\n",
    "p = \\frac{\\beta}{\\alpha + \\beta}\n",
    "$$\n",
    "\n",
    "This is, in some sense, a steady state probability of unemployment.\n",
    "\n",
    "Not surprisingly it tends to zero as $ \\beta \\to 0 $, and to one as $ \\alpha \\to 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11ab0244",
   "metadata": {},
   "source": [
    "### Calculating stationary distributions\n",
    "\n",
    "A stable algorithm for computing stationary distributions is implemented in [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
    "\n",
    "Here’s an example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c6ff7c0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[0.4, 0.6],\n",
    "     [0.2, 0.8]]\n",
    "\n",
    "mc = qe.MarkovChain(P)\n",
    "mc.stationary_distributions  # Show all stationary distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca3ac3bd",
   "metadata": {},
   "source": [
    "### Asymptotic stationarity\n",
    "\n",
    "Consider an everywhere positive stochastic matrix with unique stationary distribution $ \\psi^* $.\n",
    "\n",
    "Sometimes the distribution $ \\psi_t = \\psi_0 P^t $ of $ X_t $ converges to $ \\psi^* $ regardless of $ \\psi_0 $.\n",
    "\n",
    "For example, we have the following result\n",
    "\n",
    "\n",
    "<a id='strict-stationary'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0466207",
   "metadata": {},
   "source": [
    "### \n",
    "\n",
    "If there exists an integer $ m $ such that all entries of $ P^m $ are\n",
    "strictly positive, then\n",
    "\n",
    "$$\n",
    "\\psi_0 P^t \\to \\psi^*\n",
    "    \\quad \\text{ as } t \\to \\infty\n",
    "$$\n",
    "\n",
    "where $ \\psi^* $ is the unique stationary distribution.\n",
    "\n",
    "This situation is often referred to as **asymptotic stationarity** or **global stability**.\n",
    "\n",
    "A proof of the theorem can be found in Chapter 4 of [[Sargent and Stachurski, 2023](https://intro.quantecon.org/zreferences.html#id24)], as well as many other sources.\n",
    "\n",
    "\n",
    "<a id='hamilton'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec5cc63",
   "metadata": {},
   "source": [
    "#### Example: Hamilton’s chain\n",
    "\n",
    "Hamilton’s chain satisfies the conditions of the theorem because $ P^2 $ is everywhere positive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88435e8a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = np.array([[0.971, 0.029, 0.000],\n",
    "              [0.145, 0.778, 0.077],\n",
    "              [0.000, 0.508, 0.492]])\n",
    "P @ P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3720c15c",
   "metadata": {},
   "source": [
    "Let’s pick an initial distribution $ \\psi_1, \\psi_2, \\psi_3 $ and trace out the sequence of distributions $ \\psi_i P^t $ for $ t = 0, 1, 2, \\ldots $, for $ i=1, 2, 3 $.\n",
    "\n",
    "First, we write a function to iterate the sequence of distributions for `ts_length` period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc17e48",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def iterate_ψ(ψ_0, P, ts_length):\n",
    "    n = len(P)\n",
    "    ψ_t = np.empty((ts_length, n))\n",
    "    ψ_t[0 ]= ψ_0\n",
    "    for t in range(1, ts_length):\n",
    "        ψ_t[t] = ψ_t[t-1] @ P\n",
    "    return ψ_t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0d21134",
   "metadata": {},
   "source": [
    "Now we plot the sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c60a15bc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ_1 = (0.0, 0.0, 1.0)\n",
    "ψ_2 = (1.0, 0.0, 0.0)\n",
    "ψ_3 = (0.0, 1.0, 0.0)                   # Three initial conditions\n",
    "colors = ['blue','red', 'green']   # Different colors for each initial point\n",
    "\n",
    "# Define the vertices of the unit simplex\n",
    "v = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]])\n",
    "\n",
    "# Define the faces of the unit simplex\n",
    "faces = [\n",
    "    [v[0], v[1], v[2]],\n",
    "    [v[0], v[1], v[3]],\n",
    "    [v[0], v[2], v[3]],\n",
    "    [v[1], v[2], v[3]]\n",
    "]\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "\n",
    "def update(n):    \n",
    "    ax.clear()\n",
    "    ax.set_xlim([0, 1])\n",
    "    ax.set_ylim([0, 1])\n",
    "    ax.set_zlim([0, 1])\n",
    "    ax.view_init(45, 45)\n",
    "    \n",
    "    simplex = Poly3DCollection(faces, alpha=0.03)\n",
    "    ax.add_collection3d(simplex)\n",
    "    \n",
    "    for idx, ψ_0 in enumerate([ψ_1, ψ_2, ψ_3]):\n",
    "        ψ_t = iterate_ψ(ψ_0, P, n+1)\n",
    "        \n",
    "        for i, point in enumerate(ψ_t):\n",
    "            ax.scatter(point[0], point[1], point[2], color=colors[idx], s=60, alpha=(i+1)/len(ψ_t))\n",
    "            \n",
    "    mc = qe.MarkovChain(P)\n",
    "    ψ_star = mc.stationary_distributions[0]\n",
    "    ax.scatter(ψ_star[0], ψ_star[1], ψ_star[2], c='yellow', s=60)\n",
    "    \n",
    "    return fig,\n",
    "\n",
    "anim = FuncAnimation(fig, update, frames=range(20), blit=False, repeat=False)\n",
    "plt.close()\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ce24668",
   "metadata": {},
   "source": [
    "Here\n",
    "\n",
    "- $ P $ is the stochastic matrix for recession and growth [considered above](#mc-eg2).  \n",
    "- The red, blue and green dots are initial marginal probability distributions  $ \\psi_1, \\psi_2, \\psi_3 $, each of which is represented as a vector in $ \\mathbb R^3 $.  \n",
    "- The transparent dots are the marginal distributions $ \\psi_i P^t $ for $ t = 1, 2, \\ldots $, for $ i=1,2,3. $.  \n",
    "- The yellow dot is $ \\psi^* $.  \n",
    "\n",
    "\n",
    "You might like to try experimenting with different initial conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1b10950",
   "metadata": {},
   "source": [
    "#### Example: failure of convergence\n",
    "\n",
    "Consider the periodic chain with stochastic matrix\n",
    "\n",
    "$$\n",
    "P = \n",
    "\\begin{bmatrix}\n",
    "    0 & 1 \\\\\n",
    "    1 & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "This matrix does not satisfy the conditions of\n",
    "[strict_stationary](https://intro.quantecon.org/markov_chains_I.html#strict-stationary) because, as you can readily check,\n",
    "\n",
    "- $ P^m = P $ when $ m $ is odd and  \n",
    "- $ P^m = I $, the identity matrix, when $ m $ is even.  \n",
    "\n",
    "\n",
    "Hence there is no $ m $ such that all elements of $ P^m $ are strictly positive.\n",
    "\n",
    "Moreover, we can see that global stability does not hold.\n",
    "\n",
    "For instance, if we start at $ \\psi_0 = (1,0) $, then $ \\psi_m = \\psi_0 P^m $ is $ (1, 0) $ when $ m $ is even and $ (0,1) $ when $ m $ is odd.\n",
    "\n",
    "We can see similar phenomena in higher dimensions.\n",
    "\n",
    "The next figure illustrates this for a periodic Markov chain with three states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f188466",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ_1 = (0.0, 0.0, 1.0)\n",
    "ψ_2 = (0.5, 0.5, 0.0)\n",
    "ψ_3 = (0.25, 0.25, 0.5)\n",
    "ψ_4 = (1/3, 1/3, 1/3)\n",
    "\n",
    "P = np.array([[0.0, 1.0, 0.0],\n",
    "              [0.0, 0.0, 1.0],\n",
    "              [1.0, 0.0, 0.0]])\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "colors = ['red','yellow', 'green', 'blue']  # Different colors for each initial point\n",
    "\n",
    "# Define the vertices of the unit simplex\n",
    "v = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]])\n",
    "\n",
    "# Define the faces of the unit simplex\n",
    "faces = [\n",
    "    [v[0], v[1], v[2]],\n",
    "    [v[0], v[1], v[3]],\n",
    "    [v[0], v[2], v[3]],\n",
    "    [v[1], v[2], v[3]]\n",
    "]\n",
    "\n",
    "def update(n):\n",
    "    ax.clear()\n",
    "    ax.set_xlim([0, 1])\n",
    "    ax.set_ylim([0, 1])\n",
    "    ax.set_zlim([0, 1])\n",
    "    ax.view_init(45, 45)\n",
    "    \n",
    "    # Plot the 3D unit simplex as planes\n",
    "    simplex = Poly3DCollection(faces,alpha=0.05)\n",
    "    ax.add_collection3d(simplex)\n",
    "    \n",
    "    for idx, ψ_0 in enumerate([ψ_1, ψ_2, ψ_3, ψ_4]):\n",
    "        ψ_t = iterate_ψ(ψ_0, P, n+1)\n",
    "        \n",
    "        point = ψ_t[-1]\n",
    "        ax.scatter(point[0], point[1], point[2], color=colors[idx], s=60)\n",
    "        points = np.array(ψ_t)\n",
    "        ax.plot(points[:, 0], points[:, 1], points[:, 2], color=colors[idx],linewidth=0.75)\n",
    "    \n",
    "    return fig,\n",
    "\n",
    "anim = FuncAnimation(fig, update, frames=range(20), blit=False, repeat=False)\n",
    "plt.close()\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "952f59f2",
   "metadata": {},
   "source": [
    "This animation demonstrates the behavior of an irreducible and periodic stochastic matrix.\n",
    "\n",
    "The red, yellow, and green dots represent different initial probability distributions.\n",
    "\n",
    "The blue dot represents the unique stationary distribution.\n",
    "\n",
    "Unlike Hamilton’s Markov chain, these initial distributions do not converge to the unique stationary distribution.\n",
    "\n",
    "Instead, they cycle periodically around the probability simplex, illustrating that asymptotic stability fails.\n",
    "\n",
    "\n",
    "<a id='finite-mc-expec'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de6cebf6",
   "metadata": {},
   "source": [
    "## Computing expectations\n",
    "\n",
    "\n",
    "<a id='index-3'></a>\n",
    "We sometimes want to  compute mathematical  expectations of functions of $ X_t $ of the form\n",
    "\n",
    "\n",
    "<a id='equation-mc-une'></a>\n",
    "$$\n",
    "\\mathbb E [ h(X_t) ] \\tag{33.7}\n",
    "$$\n",
    "\n",
    "and conditional expectations such as\n",
    "\n",
    "\n",
    "<a id='equation-mc-cce'></a>\n",
    "$$\n",
    "\\mathbb E [ h(X_{t + k})  \\mid X_t = x] \\tag{33.8}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ \\{X_t\\} $ is a Markov chain generated by $ n \\times n $ stochastic matrix $ P $.  \n",
    "- $ h $ is a given function, which, in terms of matrix\n",
    "  algebra, we’ll think of as the column vector  \n",
    "\n",
    "\n",
    "$$\n",
    "h =\n",
    "\\begin{bmatrix}\n",
    "    h(x_1) \\\\\n",
    "    \\vdots \\\\\n",
    "    h(x_n)\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Computing the unconditional expectation [(33.7)](#equation-mc-une) is easy.\n",
    "\n",
    "We just sum over the marginal  distribution  of $ X_t $ to get\n",
    "\n",
    "$$\n",
    "\\mathbb E [ h(X_t) ]\n",
    "= \\sum_{x \\in S} (\\psi P^t)(x) h(x)\n",
    "$$\n",
    "\n",
    "Here $ \\psi $ is the distribution of $ X_0 $.\n",
    "\n",
    "Since $ \\psi $ and hence $ \\psi P^t $ are row vectors, we can also\n",
    "write this as\n",
    "\n",
    "$$\n",
    "\\mathbb E [ h(X_t) ]\n",
    "=  \\psi P^t h\n",
    "$$\n",
    "\n",
    "For the conditional expectation [(33.8)](#equation-mc-cce), we need to sum over\n",
    "the conditional distribution of $ X_{t + k} $ given $ X_t = x $.\n",
    "\n",
    "We already know that this is $ P^k(x, \\cdot) $, so\n",
    "\n",
    "\n",
    "<a id='equation-mc-cce2'></a>\n",
    "$$\n",
    "\\mathbb E [ h(X_{t + k})  \\mid X_t = x]\n",
    "= (P^k h)(x) \\tag{33.9}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56c0d9ab",
   "metadata": {},
   "source": [
    "### Expectations of geometric sums\n",
    "\n",
    "Sometimes we want to compute the mathematical expectation of a geometric sum, such as\n",
    "$ \\sum_t \\beta^t h(X_t) $.\n",
    "\n",
    "In view of the preceding discussion, this is\n",
    "\n",
    "$$\n",
    "\\mathbb{E}\n",
    "    \\left[\n",
    "        \\sum_{j=0}^\\infty \\beta^j h(X_{t+j}) \\mid X_t\n",
    "        = x\n",
    "    \\right]\n",
    "    = x + \\beta (Ph)(x) + \\beta^2 (P^2 h)(x) + \\cdots\n",
    "$$\n",
    "\n",
    "By the [Neumann series lemma](https://intro.quantecon.org/eigen_I.html#la-neumann), this sum can be calculated using\n",
    "\n",
    "$$\n",
    "I + \\beta P + \\beta^2 P^2 + \\cdots = (I - \\beta P)^{-1}\n",
    "$$\n",
    "\n",
    "The vector $ P^k h $ stores the conditional expectation $ \\mathbb E [ h(X_{t + k})  \\mid X_t = x] $ over all $ x $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc18a344",
   "metadata": {},
   "source": [
    "### Exercise 33.1\n",
    "\n",
    "Imam and Temple [[Imam and Temple, 2023](https://intro.quantecon.org/zreferences.html#id289)] used a three-state transition matrix to describe the transition of three states of a regime: growth, stagnation, and collapse\n",
    "\n",
    "$$\n",
    "P :=\n",
    "\\begin{bmatrix}\n",
    "    0.68 & 0.12 & 0.20 \\\\\n",
    "    0.50 & 0.24 & 0.26 \\\\\n",
    "    0.36 & 0.18 & 0.46\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where rows, from top to down, correspond to growth, stagnation, and collapse.\n",
    "\n",
    "In this exercise,\n",
    "\n",
    "1. visualize the transition matrix and show this process is asymptotically stationary  \n",
    "1. calculate the stationary distribution using simulations  \n",
    "1. visualize the dynamics of  $ (\\psi_0 P^t)(i) $ where $ t \\in 0, ..., 25 $ and compare the convergent path with the previous transition matrix  \n",
    "\n",
    "\n",
    "Compare your solution to the paper."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaeaca32",
   "metadata": {},
   "source": [
    "### Solution to[ Exercise 33.1](https://intro.quantecon.org/#mc1_ex_1)\n",
    "\n",
    "Solution 1:\n",
    "\n",
    "![https://intro.quantecon.org/_static/lecture_specific/markov_chains_I/Temple.png](https://intro.quantecon.org/_static/lecture_specific/markov_chains_I/Temple.png)\n",
    "\n",
    "Since the matrix is everywhere positive, there is a unique stationary distribution $ \\psi^* $ such that $ \\psi_t\\to \\psi^* $ as $ t\\to \\infty $.\n",
    "\n",
    "Solution 2:\n",
    "\n",
    "One simple way to calculate the stationary distribution is to take the power of the transition matrix as we have shown before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d26d581c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = np.array([[0.68, 0.12, 0.20],\n",
    "              [0.50, 0.24, 0.26],\n",
    "              [0.36, 0.18, 0.46]])\n",
    "P_power = np.linalg.matrix_power(P, 20)\n",
    "P_power"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb72070",
   "metadata": {},
   "source": [
    "Note that rows of the transition matrix converge to the stationary distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7919fe04",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ_star_p = P_power[0]\n",
    "ψ_star_p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c4daf55",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc = qe.MarkovChain(P)\n",
    "ψ_star = mc.stationary_distributions[0]\n",
    "ψ_star"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9db291a",
   "metadata": {},
   "source": [
    "### Exercise 33.2\n",
    "\n",
    "We discussed the six-state transition matrix estimated by Imam & Temple [[Imam and Temple, 2023](https://intro.quantecon.org/zreferences.html#id289)] [before](#mc-eg3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba784db1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "nodes = ['DG', 'DC', 'NG', 'NC', 'AG', 'AC']\n",
    "P = [[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],\n",
    "     [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],\n",
    "     [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],\n",
    "     [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],\n",
    "     [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],\n",
    "     [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d93c625",
   "metadata": {},
   "source": [
    "In this exercise,\n",
    "\n",
    "1. show this process is asymptotically stationary without simulation  \n",
    "1. simulate and visualize the dynamics starting with a uniform distribution across states (each state will have a probability of 1/6)  \n",
    "1. change the initial distribution to P(DG) = 1, while all other states have a probability of 0  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92418774",
   "metadata": {},
   "source": [
    "### Solution to[ Exercise 33.2](https://intro.quantecon.org/#mc1_ex_2)\n",
    "\n",
    "Solution 1:\n",
    "\n",
    "Although $ P $ is not every positive, $ P^m $ when $ m=3 $ is everywhere positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e7c90e2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = np.array([[0.86, 0.11, 0.03, 0.00, 0.00, 0.00],\n",
    "              [0.52, 0.33, 0.13, 0.02, 0.00, 0.00],\n",
    "              [0.12, 0.03, 0.70, 0.11, 0.03, 0.01],\n",
    "              [0.13, 0.02, 0.35, 0.36, 0.10, 0.04],\n",
    "              [0.00, 0.00, 0.09, 0.11, 0.55, 0.25],\n",
    "              [0.00, 0.00, 0.09, 0.15, 0.26, 0.50]])\n",
    "\n",
    "np.linalg.matrix_power(P,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8355ef8f",
   "metadata": {},
   "source": [
    "So it satisfies the requirement.\n",
    "\n",
    "Solution 2:\n",
    "\n",
    "We find the distribution $ \\psi $ converges to the stationary distribution quickly regardless of the initial distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0414ab9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ts_length = 30\n",
    "num_distributions = 20\n",
    "nodes = ['DG', 'DC', 'NG', 'NC', 'AG', 'AC']\n",
    "\n",
    "# Get parameters of transition matrix\n",
    "n = len(P)\n",
    "mc = qe.MarkovChain(P)\n",
    "ψ_star = mc.stationary_distributions[0]\n",
    "ψ_0 = np.array([[1/6 for i in range(6)],\n",
    "                [0 if i != 0 else 1 for i in range(6)]])\n",
    "## Draw the plot\n",
    "fig, axes = plt.subplots(ncols=2)\n",
    "plt.subplots_adjust(wspace=0.35)\n",
    "for idx in range(2):\n",
    "    ψ_t = iterate_ψ(ψ_0[idx], P, ts_length)\n",
    "    for i in range(n):\n",
    "        axes[idx].plot(ψ_t[:, i] - ψ_star[i], alpha=0.5, label=fr'$\\psi_t({i+1})$')\n",
    "        axes[idx].set_ylim([-0.3, 0.3])\n",
    "        axes[idx].set_xlabel('t')\n",
    "        axes[idx].set_ylabel(fr'$\\psi_t$')\n",
    "        axes[idx].legend()\n",
    "        axes[idx].axhline(0, linestyle='dashed', lw=1, color = 'black')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a22cfd",
   "metadata": {},
   "source": [
    "### Exercise 33.3\n",
    "\n",
    "Prove the following: If $ P $ is a stochastic matrix, then so is the $ k $-th\n",
    "power $ P^k $ for all $ k \\in \\mathbb N $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d578c5c",
   "metadata": {},
   "source": [
    "### Solution to[ Exercise 33.3](https://intro.quantecon.org/#mc1_ex_3)\n",
    "\n",
    "Suppose that $ P $ is stochastic and, moreover, that $ P^k $ is\n",
    "stochastic for some integer $ k $.\n",
    "\n",
    "We will prove that $ P^{k+1} = P P^k $ is also stochastic.\n",
    "\n",
    "(We are doing proof by induction — we assume the claim is true at $ k $ and\n",
    "now prove it is true at $ k+1 $.)\n",
    "\n",
    "To see this, observe that, since $ P^k $ is stochastic and the product of\n",
    "nonnegative matrices is nonnegative, $ P^{k+1} = P P^k $ is nonnegative.\n",
    "\n",
    "Also, if $ \\mathbf 1 $ is a column vector of ones, then, since $ P^k $ is stochastic we\n",
    "have $ P^k \\mathbf 1 = \\mathbf 1 $ (rows sum to one).\n",
    "\n",
    "Therefore $ P^{k+1} \\mathbf 1 = P P^k \\mathbf 1 = P \\mathbf 1 = \\mathbf 1 $\n",
    "\n",
    "The proof is done."
   ]
  }
 ],
 "metadata": {
  "date": 1730689203.6814911,
  "filename": "markov_chains_I.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Markov Chains: Basic Concepts"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}