finite_markov.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b2fd630a",
   "metadata": {},
   "source": [
    "\n",
    "<a id='mc'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed3894d7",
   "metadata": {},
   "source": [
    "# Finite Markov Chains\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72c4f3f9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68e50683",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Markov chains are one of the most useful classes of stochastic processes, being\n",
    "\n",
    "- simple, flexible and supported by many elegant theoretical results  \n",
    "- valuable for building intuition about random dynamic models  \n",
    "- central to quantitative modeling in their own right  \n",
    "\n",
    "\n",
    "You will find them in many of the workhorse models of economics and finance.\n",
    "\n",
    "In this lecture, we review some of the theory of Markov chains.\n",
    "\n",
    "We will also introduce some of the high-quality routines for working with Markov chains available in [QuantEcon.py](https://quantecon.org/quantecon-py/).\n",
    "\n",
    "Prerequisite knowledge is basic probability and linear algebra.\n",
    "\n",
    "Let’s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4e18268",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)  #set default figure size\n",
    "import quantecon as qe\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4991731a",
   "metadata": {},
   "source": [
    "## Definitions\n",
    "\n",
    "The following concepts are fundamental.\n",
    "\n",
    "\n",
    "<a id='finite-dp-stoch-mat'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4467a6e1",
   "metadata": {},
   "source": [
    "### Stochastic Matrices\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "A **stochastic matrix** (or **Markov matrix**)  is an $ n \\times n $ square matrix $ P $\n",
    "such that\n",
    "\n",
    "1. each element of $ P $ is nonnegative, and  \n",
    "1. each row of $ P $ sums to one  \n",
    "\n",
    "\n",
    "Each row of $ P $ can be regarded as a probability mass function over $ n $ possible outcomes.\n",
    "\n",
    "It is too not difficult to check <sup><a href=#pm id=pm-link>[1]</a></sup> that if $ P $ is a stochastic matrix, then so is the $ k $-th power $ P^k $ for all $ k \\in \\mathbb N $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f47083a",
   "metadata": {},
   "source": [
    "### Markov Chains\n",
    "\n",
    "\n",
    "<a id='index-4'></a>\n",
    "There is a close connection between stochastic matrices and Markov chains.\n",
    "\n",
    "To begin, let $ S $ be a finite set with $ n $ elements $ \\{x_1, \\ldots, x_n\\} $.\n",
    "\n",
    "The set $ S $ is called the **state space** and $ x_1, \\ldots, x_n $ are the **state values**.\n",
    "\n",
    "A **Markov chain** $ \\{X_t\\} $ on $ S $ is a sequence of random variables on $ S $ 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{13.1}\n",
    "$$\n",
    "\n",
    "In other words, knowing the current state is enough to know probabilities for future states.\n",
    "\n",
    "In particular, the dynamics of a Markov chain are fully determined by the set of values\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{13.2}\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 marginal distribution $ \\psi $  \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 [(13.2)](#equation-mpp).\n",
    "\n",
    "\n",
    "<a id='mc-eg1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0f92235",
   "metadata": {},
   "source": [
    "### Example 1\n",
    "\n",
    "Consider a worker who, at any given time $ t $, is either unemployed (state 0) or employed (state 1).\n",
    "\n",
    "Suppose that, over a one month period,\n",
    "\n",
    "1. An unemployed worker finds a job with probability $ \\alpha \\in (0, 1) $.  \n",
    "1. An employed worker loses her job and becomes unemployed with probability $ \\beta \\in (0, 1) $.  \n",
    "\n",
    "\n",
    "In terms of a Markov model, we have\n",
    "\n",
    "- $ S = \\{ 0, 1\\} $  \n",
    "- $ P(0, 1) = \\alpha $ and $ P(1, 0) = \\beta $  \n",
    "\n",
    "\n",
    "We can write out the transition probabilities in matrix form as\n",
    "\n",
    "\n",
    "<a id='equation-p-unempemp'></a>\n",
    "$$\n",
    "P\n",
    "= \\left(\n",
    "\\begin{array}{cc}\n",
    "    1 - \\alpha & \\alpha \\\\\n",
    "    \\beta & 1 - \\beta\n",
    "\\end{array}\n",
    "  \\right) \\tag{13.3}\n",
    "$$\n",
    "\n",
    "Once we have the values $ \\alpha $ and $ \\beta $, 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 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 such applications below.\n",
    "\n",
    "\n",
    "<a id='mc-eg2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1189c768",
   "metadata": {},
   "source": [
    "### Example 2\n",
    "\n",
    "From  US unemployment data, Hamilton [[Hamilton, 2005](https://tools-techniques.quantecon.org/zreferences.html#id155)] estimated the stochastic matrix\n",
    "\n",
    "$$\n",
    "P =\n",
    "\\left(\n",
    "  \\begin{array}{ccc}\n",
    "     0.971 & 0.029 & 0 \\\\\n",
    "     0.145 & 0.778 & 0.077 \\\\\n",
    "     0 & 0.508 & 0.492\n",
    "  \\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- the frequency is monthly  \n",
    "- the first state represents “normal growth”  \n",
    "- the second state represents “mild recession”  \n",
    "- the third state represents “severe recession”  \n",
    "\n",
    "\n",
    "For example, the matrix tells us that when the state is normal growth, the state will again be normal growth next month with probability 0.97.\n",
    "\n",
    "In general, large values on the main diagonal indicate persistence in the process $ \\{ X_t \\} $.\n",
    "\n",
    "This Markov process can also be represented as a directed graph, with edges labeled by transition probabilities\n",
    "\n",
    "![https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/hamilton_graph.png](https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/hamilton_graph.png)\n",
    "\n",
    "  \n",
    "Here “ng” is normal growth, “mr” is mild recession, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91157896",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "\n",
    "<a id='index-5'></a>\n",
    "One natural way to answer questions about Markov chains is to simulate them.\n",
    "\n",
    "(To approximate the probability of event $ E $, we can simulate many times and count the fraction of times that $ E $ occurs).\n",
    "\n",
    "Nice functionality for simulating Markov chains exists in [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
    "\n",
    "- Efficient, bundled with lots of other useful routines for handling Markov chains.  \n",
    "\n",
    "\n",
    "However, it’s also a good exercise to roll our own routines — let’s do that first and then come back to the methods in [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
    "\n",
    "In these exercises, we’ll take the state space to be $ S = 0,\\ldots, n-1 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3784a73b",
   "metadata": {},
   "source": [
    "### Rolling Our Own\n",
    "\n",
    "To simulate a Markov chain, we need its stochastic matrix $ P $ and a marginal probability distribution $ \\psi $  from which to  draw a realization of $ X_0 $.\n",
    "\n",
    "The Markov chain is then constructed as discussed above.  To repeat:\n",
    "\n",
    "1. At time $ t=0 $, draw a realization of  $ X_0 $  from $ \\psi $.  \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",
    "To implement this simulation procedure, we need a method for generating draws from a discrete distribution.\n",
    "\n",
    "For this task, we’ll use `random.draw` from [QuantEcon](http://quantecon.org/quantecon-py), which works as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f2147b0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ = (0.3, 0.7)           # probabilities over {0, 1}\n",
    "cdf = np.cumsum(ψ)       # convert into cummulative distribution\n",
    "qe.random.draw(cdf, 5)   # generate 5 independent draws from ψ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a54219f3",
   "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 state `init`  \n",
    "- A positive integer `sample_size` representing the length of the time series the function should return  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88fffeb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def mc_sample_path(P, ψ_0=None, sample_size=1_000):\n",
    "\n",
    "    # set up\n",
    "    P = np.asarray(P)\n",
    "    X = np.empty(sample_size, dtype=int)\n",
    "\n",
    "    # Convert each row of P into a cdf\n",
    "    n = len(P)\n",
    "    P_dist = [np.cumsum(P[i, :]) for i in range(n)]\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(sample_size - 1):\n",
    "        X[t+1] = qe.random.draw(P_dist[X[t]])\n",
    "\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "133a9310",
   "metadata": {},
   "source": [
    "Let’s see how it works using the small matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fb412fb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[0.4, 0.6],\n",
    "     [0.2, 0.8]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f94df54a",
   "metadata": {},
   "source": [
    "As we’ll see later, for a long series drawn from `P`, the fraction of the sample that takes value 0 will be about 0.25.\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": "644dee7a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = mc_sample_path(P, ψ_0=[0.1, 0.9], sample_size=100_000)\n",
    "np.mean(X == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74313c19",
   "metadata": {},
   "source": [
    "You can try changing the initial distribution to confirm that the output is\n",
    "always close to 0.25, at least for the `P` matrix above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7798cf3b",
   "metadata": {},
   "source": [
    "### Using QuantEcon’s Routines\n",
    "\n",
    "As discussed above, [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": "ac719a1b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from quantecon import MarkovChain\n",
    "\n",
    "mc = qe.MarkovChain(P)\n",
    "X = mc.simulate(ts_length=1_000_000)\n",
    "np.mean(X == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed63e044",
   "metadata": {},
   "source": [
    "The [QuantEcon.py](http://quantecon.org/quantecon-py) routine is [JIT compiled](https://python-programming.quantecon.org/numba.html#numba-link) and much faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "798a6259",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%time mc_sample_path(P, sample_size=1_000_000) # Our homemade code version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bafeb80d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%time mc.simulate(ts_length=1_000_000) # qe code version"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b48abc3",
   "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": "78456e0c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc = qe.MarkovChain(P, state_values=('unemployed', 'employed'))\n",
    "mc.simulate(ts_length=4, init='employed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4da76e7a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.simulate(ts_length=4, init='unemployed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97b20dbf",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.simulate(ts_length=4)  # Start at randomly chosen initial state"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c357d858",
   "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": "3980c205",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.simulate_indices(ts_length=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c5fb5b3",
   "metadata": {},
   "source": [
    "\n",
    "<a id='mc-md'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeff24b8",
   "metadata": {},
   "source": [
    "## Marginal Distributions\n",
    "\n",
    "\n",
    "<a id='index-7'></a>\n",
    "Suppose that\n",
    "\n",
    "1. $ \\{X_t\\} $ is a Markov chain with stochastic matrix $ P $  \n",
    "1. the marginal distribution of $ X_t $ is known to be $ \\psi_t $  \n",
    "\n",
    "\n",
    "What then is the marginal distribution of $ X_{t+1} $, or, more generally, of $ X_{t+m} $?\n",
    "\n",
    "To answer this, we let $ \\psi_t $ be the marginal 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",
    "Using the [law of total probability](https://en.wikipedia.org/wiki/Law_of_total_probability), we can decompose the probability that $ X_{t+1} = y $ as follows:\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",
    "In words, to get the probability of being at $ y $ tomorrow, we account for\n",
    "all  ways this can happen and sum their probabilities.\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{13.4}\n",
    "$$\n",
    "\n",
    "Thus, to move a marginal distribution forward one unit of time, we postmultiply by $ P $.\n",
    "\n",
    "By postmultiplying $ m $ times, we move a marginal distribution forward $ m $ steps into the future.\n",
    "\n",
    "Hence, iterating on [(13.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{13.5}\n",
    "$$\n",
    "\n",
    "and, more generally,\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{13.6}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='finite-mc-mstp'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58292b5d",
   "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 [(13.6)](#equation-mdfmc2), but now with a $ \\psi_t $ that puts all probability on state $ x $ so that the transition probabilities are\n",
    "\n",
    "- 1 in the $ x $-th position and zero elsewhere  \n",
    "\n",
    "\n",
    "Inserting this into [(13.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": "d455a90a",
   "metadata": {},
   "source": [
    "### Example: Probability of Recession\n",
    "\n",
    "\n",
    "<a id='index-8'></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(x) $.\n",
    "\n",
    "The probability of being in recession (either mild or severe) in 6 months time is given by the inner product\n",
    "\n",
    "$$\n",
    "\\psi P^6\n",
    "\\cdot\n",
    "\\left(\n",
    "  \\begin{array}{c}\n",
    "     0 \\\\\n",
    "     1 \\\\\n",
    "     1\n",
    "  \\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='mc-eg1-1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "006a80d5",
   "metadata": {},
   "source": [
    "### Example 2: Cross-Sectional Distributions\n",
    "\n",
    "\n",
    "<a id='index-9'></a>\n",
    "The marginal distributions we have been studying can be viewed either as\n",
    "probabilities or as cross-sectional frequencies that a Law of Large Numbers leads us to anticipate for  large samples.\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 described by the specified dynamics, with each worker’s\n",
    "outcomes being realizations of processes that are statistically independent of all other workers’ processes.\n",
    "\n",
    "Let $ \\psi $ 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.\n",
    "\n",
    "- For example, $ \\psi(0) $ is the unemployment rate.  \n",
    "\n",
    "\n",
    "What will the cross-sectional distribution be in 10 periods hence?\n",
    "\n",
    "The answer is $ \\psi P^{10} $, where $ P $ is the stochastic matrix in\n",
    "[(13.3)](#equation-p-unempemp).\n",
    "\n",
    "This is because each worker’s state evolves according to $ P $, so\n",
    "$ \\psi P^{10} $ is a marginal distibution  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 P^{10} $ also represents  fractions of workers in\n",
    "each state.\n",
    "\n",
    "This is exactly the cross-sectional distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d11ae1b",
   "metadata": {},
   "source": [
    "## Irreducibility and Aperiodicity\n",
    "\n",
    "\n",
    "<a id='index-11'></a>\n",
    "Irreducibility and aperiodicity are central concepts of modern Markov chain theory.\n",
    "\n",
    "Let’s see what they’re about."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1d7cbb9",
   "metadata": {},
   "source": [
    "### Irreducibility\n",
    "\n",
    "Let $ P $ be a fixed stochastic matrix.\n",
    "\n",
    "Two states $ x $ and $ y $ are said to **communicate** with each other if\n",
    "there exist positive integers $ j $ and $ k $ such that\n",
    "\n",
    "$$\n",
    "P^j(x, y) > 0\n",
    "\\quad \\text{and} \\quad\n",
    "P^k(y, x) > 0\n",
    "$$\n",
    "\n",
    "In view of our discussion [above](#finite-mc-mstp), this means precisely\n",
    "that\n",
    "\n",
    "- state $ x $ can eventually be reached  from state $ y $, and  \n",
    "- state $ y $ can eventually  be reached from state $ x $  \n",
    "\n",
    "\n",
    "The stochastic matrix $ P $ is called **irreducible** if all states\n",
    "communicate; that is, if $ x $ and $ y $ communicate for all\n",
    "$ (x, y) $ in $ S \\times S $.\n",
    "\n",
    "For example, consider the following transition probabilities for wealth of a fictitious set of\n",
    "households\n",
    "\n",
    "![https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_irreducibility1.png](https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_irreducibility1.png)\n",
    "\n",
    "  \n",
    "We can translate this into a stochastic matrix, putting zeros where\n",
    "there’s no edge between nodes\n",
    "\n",
    "$$\n",
    "P :=\n",
    "\\left(\n",
    "  \\begin{array}{ccc}\n",
    "     0.9 & 0.1 & 0 \\\\\n",
    "     0.4 & 0.4 & 0.2 \\\\\n",
    "     0.1 & 0.1 & 0.8\n",
    "  \\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "It’s clear from the graph that this stochastic matrix is irreducible: we can  eventually\n",
    "reach any state from any other state.\n",
    "\n",
    "We can also test this using [QuantEcon.py](http://quantecon.org/quantecon-py)’s MarkovChain class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0534469",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[0.9, 0.1, 0.0],\n",
    "     [0.4, 0.4, 0.2],\n",
    "     [0.1, 0.1, 0.8]]\n",
    "\n",
    "mc = qe.MarkovChain(P, ('poor', 'middle', 'rich'))\n",
    "mc.is_irreducible"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be027f68",
   "metadata": {},
   "source": [
    "Here’s a more pessimistic scenario in which  poor people remain poor forever\n",
    "\n",
    "![https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_irreducibility2.png](https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_irreducibility2.png)\n",
    "\n",
    "  \n",
    "This stochastic matrix is not irreducible, since, for example, rich is not accessible from poor.\n",
    "\n",
    "Let’s confirm this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38af2da6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[1.0, 0.0, 0.0],\n",
    "     [0.1, 0.8, 0.1],\n",
    "     [0.0, 0.2, 0.8]]\n",
    "\n",
    "mc = qe.MarkovChain(P, ('poor', 'middle', 'rich'))\n",
    "mc.is_irreducible"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac63d147",
   "metadata": {},
   "source": [
    "We can also determine the “communication classes”"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15a6864f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.communication_classes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43c6dba8",
   "metadata": {},
   "source": [
    "It might be clear to you already that irreducibility is going to be important in terms of long run outcomes.\n",
    "\n",
    "For example, poverty is a life sentence in the second graph but not the first.\n",
    "\n",
    "We’ll come back to this a bit later."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe702e8",
   "metadata": {},
   "source": [
    "### Aperiodicity\n",
    "\n",
    "Loosely speaking, a Markov chain is called **periodic** if it cycles in a predictable way, and **aperiodic** otherwise.\n",
    "\n",
    "Here’s a trivial example with three states\n",
    "\n",
    "![https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_aperiodicity1.png](https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_aperiodicity1.png)\n",
    "\n",
    "  \n",
    "The chain cycles with period 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1efed039",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[0, 1, 0],\n",
    "     [0, 0, 1],\n",
    "     [1, 0, 0]]\n",
    "\n",
    "mc = qe.MarkovChain(P)\n",
    "mc.period"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee5e2682",
   "metadata": {},
   "source": [
    "More formally, the **period** of a state $ x $ is the largest common divisor\n",
    "of a set of integers\n",
    "\n",
    "$$\n",
    "D(x) := \\{j \\geq 1 : P^j(x, x) > 0\\}\n",
    "$$\n",
    "\n",
    "In the last example, $ D(x) = \\{3, 6, 9, \\ldots\\} $ for every state $ x $, so the period is 3.\n",
    "\n",
    "A stochastic matrix is called **aperiodic** if the period of every state is 1, and **periodic** otherwise.\n",
    "\n",
    "For example, the stochastic matrix associated with the transition probabilities below is periodic because, for example, state $ a $ has period 2\n",
    "\n",
    "![https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_aperiodicity2.png](https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/mc_aperiodicity2.png)\n",
    "\n",
    "  \n",
    "We can confirm that the stochastic matrix is periodic with the following code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "922b6fce",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = [[0.0, 1.0, 0.0, 0.0],\n",
    "     [0.5, 0.0, 0.5, 0.0],\n",
    "     [0.0, 0.5, 0.0, 0.5],\n",
    "     [0.0, 0.0, 1.0, 0.0]]\n",
    "\n",
    "mc = qe.MarkovChain(P)\n",
    "mc.period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c9b44f2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mc.is_aperiodic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f65a239",
   "metadata": {},
   "source": [
    "## Stationary Distributions\n",
    "\n",
    "\n",
    "<a id='index-13'></a>\n",
    "As seen in [(13.4)](#equation-fin-mc-fr), we can shift a marginal distribution forward one 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": "33345e6b",
   "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": "d1b1509b",
   "metadata": {},
   "source": [
    "Such distributions are called **stationary** or **invariant**.\n",
    "\n",
    "\n",
    "<a id='mc-stat-dd'></a>\n",
    "Formally, a marginal distribution $ \\psi^* $ on $ S $ is called **stationary** for $ P $ if $ \\psi^* = \\psi^* P $.\n",
    "\n",
    "(This is the same notion of stationarity that we learned about in the\n",
    "[lecture on AR(1) processes](https://intro.quantecon.org/ar1_processes.html) applied to a different setting.)\n",
    "\n",
    "From this equality, we immediately get $ \\psi^* = \\psi^* P^t $ for all $ t $.\n",
    "\n",
    "This tells us an important fact: If the distribution of $ X_0 $ is a stationary distribution, then $ X_t $ will have this same distribution for all $ t $.\n",
    "\n",
    "Hence stationary distributions have a natural interpretation as **stochastic steady states** — we’ll discuss this more soon.\n",
    "\n",
    "Mathematically, a stationary distribution is a fixed point of $ P $ when $ P $ is thought of as the map $ \\psi \\mapsto \\psi P $ from (row) vectors to (row) vectors.\n",
    "\n",
    "**Theorem.** Every stochastic matrix $ P $ has at least one stationary distribution.\n",
    "\n",
    "(We are assuming here that the state space $ S $ is finite; if not more assumptions are required)\n",
    "\n",
    "For proof of this result, you can apply [Brouwer’s fixed point theorem](https://en.wikipedia.org/wiki/Brouwer_fixed-point_theorem), or see [EDTC](https://johnstachurski.net/edtc.html), theorem 4.3.5.\n",
    "\n",
    "There can be many stationary distributions corresponding to a given stochastic matrix $ P $.\n",
    "\n",
    "- For example, if $ P $ is the identity matrix, then all marginal distributions are stationary.  \n",
    "\n",
    "\n",
    "To get uniqueness an invariant distribution, the transition matrix $ P $ must have the property that no nontrivial subsets of\n",
    "the state space are **infinitely persistent**.\n",
    "\n",
    "A subset of the state space is infinitely persistent if other parts of the\n",
    "state space cannot be accessed from it.\n",
    "\n",
    "Thus, infinite persistence of a non-trivial subset is the opposite of irreducibility.\n",
    "\n",
    "This gives some intuition for the following fundamental theorem.\n",
    "\n",
    "\n",
    "<a id='mc-conv-thm'></a>\n",
    "**Theorem.** If $ P $ is both aperiodic and irreducible, then\n",
    "\n",
    "1. $ P $ has exactly one stationary distribution $ \\psi^* $.  \n",
    "1. For any initial marginal distribution $ \\psi_0 $, we have $ \\| \\psi_0 P^t - \\psi^* \\| \\to 0 $ as $ t \\to \\infty $.  \n",
    "\n",
    "\n",
    "For a proof, see, for example, theorem 5.2 of [[Häggström, 2002](https://tools-techniques.quantecon.org/zreferences.html#id123)].\n",
    "\n",
    "(Note that part 1 of the theorem only requires  irreducibility, whereas part 2\n",
    "requires both irreducibility and aperiodicity)\n",
    "\n",
    "A stochastic matrix that satisfies the conditions of the theorem is sometimes called **uniformly ergodic**.\n",
    "\n",
    "A sufficient condition for aperiodicity and irreducibility is that every element of $ P $ is strictly positive.\n",
    "\n",
    "- Try to convince yourself of this.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e20b0b3",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Recall our model of the employment/unemployment dynamics of a particular worker [discussed above](#mc-eg1).\n",
    "\n",
    "Assuming $ \\alpha \\in (0,1) $ and $ \\beta \\in (0,1) $, the uniform ergodicity condition is satisfied.\n",
    "\n",
    "Let $ \\psi^* = (p, 1-p) $ be the stationary distribution, so that $ p $ 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 — more about the  interpretation of this below.\n",
    "\n",
    "Not surprisingly it tends to zero as $ \\beta \\to 0 $, and to one as $ \\alpha \\to 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42760c33",
   "metadata": {},
   "source": [
    "### Calculating Stationary Distributions\n",
    "\n",
    "\n",
    "<a id='index-14'></a>\n",
    "As discussed above, a particular Markov matrix $ P $ can have many stationary distributions.\n",
    "\n",
    "That is, there can be many row vectors $ \\psi $ such that $ \\psi = \\psi P $.\n",
    "\n",
    "In fact if $ P $ has two distinct stationary distributions $ \\psi_1,\n",
    "\\psi_2 $ then it has infinitely many, since in this case, as you can verify,  for any $ \\lambda \\in [0, 1] $\n",
    "\n",
    "$$\n",
    "\\psi_3 := \\lambda \\psi_1 + (1 - \\lambda) \\psi_2\n",
    "$$\n",
    "\n",
    "is a stationary distribution for $ P $.\n",
    "\n",
    "If we restrict attention to the case in which only one stationary distribution exists, one way to  finding it is to solve the system\n",
    "\n",
    "\n",
    "<a id='equation-eq-eqpsifixed'></a>\n",
    "$$\n",
    "\\psi (I_n - P) = 0 \\tag{13.7}\n",
    "$$\n",
    "\n",
    "for $ \\psi $, where $ I_n $ is the $ n \\times n $ identity.\n",
    "\n",
    "But the zero vector solves system [(13.7)](#equation-eq-eqpsifixed),  so we must proceed cautiously.\n",
    "\n",
    "We want to impose the restriction that $ \\psi $ is  a probability distribution.\n",
    "\n",
    "There are various ways to do this.\n",
    "\n",
    "One option is to regard solving system [(13.7)](#equation-eq-eqpsifixed)  as an eigenvector problem: a vector\n",
    "$ \\psi $ such that $ \\psi = \\psi P $ is a left eigenvector associated\n",
    "with the unit eigenvalue $ \\lambda = 1 $.\n",
    "\n",
    "A stable and sophisticated algorithm specialized for stochastic matrices is implemented in [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
    "\n",
    "This is the one we recommend:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d274f00",
   "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": "7a192e2d",
   "metadata": {},
   "source": [
    "### Convergence to Stationarity\n",
    "\n",
    "\n",
    "<a id='index-15'></a>\n",
    "Part 2 of the Markov chain convergence theorem [stated above](#mc-conv-thm) tells us that the marginal distribution of $ X_t $ converges to the stationary distribution regardless of where we begin.\n",
    "\n",
    "This adds considerable authority to our interpretation of $ \\psi^* $ as a stochastic steady state.\n",
    "\n",
    "The convergence in the theorem is illustrated in the next figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb2242bb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P = ((0.971, 0.029, 0.000),\n",
    "     (0.145, 0.778, 0.077),\n",
    "     (0.000, 0.508, 0.492))\n",
    "P = np.array(P)\n",
    "\n",
    "ψ = (0.0, 0.2, 0.8)        # Initial condition\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "ax.set(xlim=(0, 1), ylim=(0, 1), zlim=(0, 1),\n",
    "       xticks=(0.25, 0.5, 0.75),\n",
    "       yticks=(0.25, 0.5, 0.75),\n",
    "       zticks=(0.25, 0.5, 0.75))\n",
    "\n",
    "x_vals, y_vals, z_vals = [], [], []\n",
    "for t in range(20):\n",
    "    x_vals.append(ψ[0])\n",
    "    y_vals.append(ψ[1])\n",
    "    z_vals.append(ψ[2])\n",
    "    ψ = ψ @ P\n",
    "\n",
    "ax.scatter(x_vals, y_vals, z_vals, c='r', s=60)\n",
    "ax.view_init(30, 210)\n",
    "\n",
    "mc = qe.MarkovChain(P)\n",
    "ψ_star = mc.stationary_distributions[0]\n",
    "ax.scatter(ψ_star[0], ψ_star[1], ψ_star[2], c='k', s=60)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf15c95b",
   "metadata": {},
   "source": [
    "Here\n",
    "\n",
    "- $ P $ is the stochastic matrix for recession and growth [considered above](#mc-eg2).  \n",
    "- The highest red dot is an arbitrarily chosen initial marginal probability distribution  $ \\psi $, represented as a vector in $ \\mathbb R^3 $.  \n",
    "- The other red dots are the marginal distributions $ \\psi P^t $ for $ t = 1, 2, \\ldots $.  \n",
    "- The black dot is $ \\psi^* $.  \n",
    "\n",
    "\n",
    "You might like to try experimenting with different initial conditions.\n",
    "\n",
    "\n",
    "<a id='ergodicity'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d8f5c99",
   "metadata": {},
   "source": [
    "## Ergodicity\n",
    "\n",
    "\n",
    "<a id='index-17'></a>\n",
    "Under irreducibility, yet another important result obtains: for all $ x \\in S $,\n",
    "\n",
    "\n",
    "<a id='equation-llnfmc0'></a>\n",
    "$$\n",
    "\\frac{1}{m} \\sum_{t = 1}^m \\mathbf{1}\\{X_t = x\\}  \\to \\psi^*(x)\n",
    "    \\quad \\text{as } m \\to \\infty \\tag{13.8}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ \\mathbf{1}\\{X_t = x\\} = 1 $ if $ X_t = x $ and zero otherwise  \n",
    "- convergence is with probability one  \n",
    "- the result does not depend on the marginal distribution  of $ X_0 $  \n",
    "\n",
    "\n",
    "The result tells us that the fraction of time the chain spends at state $ x $ converges to $ \\psi^*(x) $ as time goes to infinity.\n",
    "\n",
    "\n",
    "<a id='new-interp-sd'></a>\n",
    "This gives us another way to interpret the stationary distribution — provided that the convergence result in [(13.8)](#equation-llnfmc0) is valid.\n",
    "\n",
    "The convergence asserted in [(13.8)](#equation-llnfmc0) is a special case of a law of large numbers result for Markov chains — see [EDTC](http://johnstachurski.net/edtc.html), section 4.3.4 for some additional information.\n",
    "\n",
    "\n",
    "<a id='mc-eg1-2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abd31922",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Recall our cross-sectional interpretation of the employment/unemployment model [discussed above](#mc-eg1-1).\n",
    "\n",
    "Assume that $ \\alpha \\in (0,1) $ and $ \\beta \\in (0,1) $, so that irreducibility and aperiodicity both hold.\n",
    "\n",
    "We saw that the stationary distribution is $ (p, 1-p) $, where\n",
    "\n",
    "$$\n",
    "p = \\frac{\\beta}{\\alpha + \\beta}\n",
    "$$\n",
    "\n",
    "In the cross-sectional interpretation, this is the fraction of people unemployed.\n",
    "\n",
    "In view of our latest (ergodicity) result, it is also the fraction of time that a single worker can expect to spend unemployed.\n",
    "\n",
    "Thus, in the long-run, cross-sectional averages for a population and time-series averages for a given person coincide.\n",
    "\n",
    "This is one aspect of the concept  of ergodicity.\n",
    "\n",
    "\n",
    "<a id='finite-mc-expec'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71da9d56",
   "metadata": {},
   "source": [
    "## Computing Expectations\n",
    "\n",
    "\n",
    "<a id='index-18'></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{13.9}\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{13.10}\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",
    "= \\left(\n",
    "\\begin{array}{c}\n",
    "    h(x_1) \\\\\n",
    "    \\vdots \\\\\n",
    "    h(x_n)\n",
    "\\end{array}\n",
    "  \\right)\n",
    "$$\n",
    "\n",
    "Computing the unconditional expectation [(13.9)](#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 [(13.10)](#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{13.11}\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": "8f860187",
   "metadata": {},
   "source": [
    "### Iterated Expectations\n",
    "\n",
    "The **law of iterated expectations** states that\n",
    "\n",
    "$$\n",
    "\\mathbb E \\left[ \\mathbb E [ h(X_{t + k})  \\mid X_t = x] \\right] = \\mathbb E [  h(X_{t + k}) ]\n",
    "$$\n",
    "\n",
    "where the outer $ \\mathbb E $ on the left side is an unconditional distribution taken with respect to the marginal distribution  $ \\psi_t $ of $ X_t $\n",
    "(again see equation [(13.6)](#equation-mdfmc2)).\n",
    "\n",
    "To verify the law of iterated expectations, use  equation [(13.11)](#equation-mc-cce2) to substitute $ (P^k h)(x) $ for $ E [ h(X_{t + k})  \\mid X_t = x] $, write\n",
    "\n",
    "$$\n",
    "\\mathbb E \\left[ \\mathbb E [ h(X_{t + k})  \\mid X_t = x] \\right] = \\psi_t P^k h,\n",
    "$$\n",
    "\n",
    "and note $ \\psi_t P^k h = \\psi_{t+k} h = \\mathbb E [  h(X_{t + k}) ] $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02e19f98",
   "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",
    "        \\sum_{j=0}^\\infty \\beta^j h(X_{t+j}) \\mid X_t = x\n",
    "    \\Bigr]\n",
    "= [(I - \\beta P)^{-1} h](x)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "(I - \\beta P)^{-1}  = I + \\beta P + \\beta^2 P^2 + \\cdots\n",
    "$$\n",
    "\n",
    "Premultiplication by $ (I - \\beta P)^{-1} $ amounts to “applying the **resolvent operator**”."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "369ef2c8",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "428492af",
   "metadata": {},
   "source": [
    "## Exercise 13.1\n",
    "\n",
    "According to the discussion [above](#mc-eg1-2), if a worker’s employment dynamics obey the stochastic matrix\n",
    "\n",
    "$$\n",
    "P\n",
    "= \\left(\n",
    "\\begin{array}{cc}\n",
    "    1 - \\alpha & \\alpha \\\\\n",
    "    \\beta & 1 - \\beta\n",
    "\\end{array}\n",
    "  \\right)\n",
    "$$\n",
    "\n",
    "with $ \\alpha \\in (0,1) $ and $ \\beta \\in (0,1) $, then, in the long-run, the fraction\n",
    "of time spent unemployed will be\n",
    "\n",
    "$$\n",
    "p := \\frac{\\beta}{\\alpha + \\beta}\n",
    "$$\n",
    "\n",
    "In other words, if $ \\{X_t\\} $ represents the Markov chain for\n",
    "employment, then $ \\bar X_m \\to p $ as $ m \\to \\infty $, where\n",
    "\n",
    "$$\n",
    "\\bar X_m := \\frac{1}{m} \\sum_{t = 1}^m \\mathbf{1}\\{X_t = 0\\}\n",
    "$$\n",
    "\n",
    "This exercise asks you to illustrate convergence by computing\n",
    "$ \\bar X_m $ for large $ m $ and checking that\n",
    "it is close to $ p $.\n",
    "\n",
    "You will see that this statement is true regardless of the choice of initial\n",
    "condition or the values of $ \\alpha, \\beta $, provided both lie in\n",
    "$ (0, 1) $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df0f5391",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 13.1](https://tools-techniques.quantecon.org/#fm_ex1)\n",
    "\n",
    "We will address this exercise graphically.\n",
    "\n",
    "The plots show the time series of $ \\bar X_m - p $ for two initial\n",
    "conditions.\n",
    "\n",
    "As $ m $ gets large, both series converge to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6dc749",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "α = β = 0.1\n",
    "N = 10000\n",
    "p = β / (α + β)\n",
    "\n",
    "P = ((1 - α,       α),               # Careful: P and p are distinct\n",
    "     (    β,   1 - β))\n",
    "mc = MarkovChain(P)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 6))\n",
    "ax.set_ylim(-0.25, 0.25)\n",
    "ax.grid()\n",
    "ax.hlines(0, 0, N, lw=2, alpha=0.6)   # Horizonal line at zero\n",
    "\n",
    "for x0, col in ((0, 'blue'), (1, 'green')):\n",
    "    # Generate time series for worker that starts at x0\n",
    "    X = mc.simulate(N, init=x0)\n",
    "    # Compute fraction of time spent unemployed, for each n\n",
    "    X_bar = (X == 0).cumsum() / (1 + np.arange(N, dtype=float))\n",
    "    # Plot\n",
    "    ax.fill_between(range(N), np.zeros(N), X_bar - p, color=col, alpha=0.1)\n",
    "    ax.plot(X_bar - p, color=col, label=f'$X_0 = \\, {x0} $')\n",
    "    # Overlay in black--make lines clearer\n",
    "    ax.plot(X_bar - p, 'k-', alpha=0.6)\n",
    "\n",
    "ax.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6af54ef",
   "metadata": {},
   "source": [
    "## Exercise 13.2\n",
    "\n",
    "A topic of interest for economics and many other disciplines is *ranking*.\n",
    "\n",
    "Let’s now consider one of the most practical and important ranking problems\n",
    "— the rank assigned to web pages by search engines.\n",
    "\n",
    "(Although the problem is motivated from outside of economics, there is in fact a deep connection between search ranking systems and prices in certain competitive equilibria — see [[Du *et al.*, 2013](https://tools-techniques.quantecon.org/zreferences.html#id147)].)\n",
    "\n",
    "To understand the issue, consider the set of results returned by a query to a web search engine.\n",
    "\n",
    "For the user, it is desirable to\n",
    "\n",
    "1. receive a large set of accurate matches  \n",
    "1. have the matches returned in order, where the order corresponds to some measure of “importance”  \n",
    "\n",
    "\n",
    "Ranking according to a measure of importance is the problem we now consider.\n",
    "\n",
    "The methodology developed to solve this problem by Google founders Larry Page and Sergey Brin\n",
    "is known as [PageRank](https://en.wikipedia.org/wiki/PageRank).\n",
    "\n",
    "To illustrate the idea, consider the following diagram\n",
    "\n",
    "![https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/web_graph.png](https://tools-techniques.quantecon.org/_static/lecture_specific/finite_markov/web_graph.png)\n",
    "\n",
    "  \n",
    "Imagine that this is a miniature version of the WWW, with\n",
    "\n",
    "- each node representing a web page  \n",
    "- each arrow representing the existence of a link from one page to another  \n",
    "\n",
    "\n",
    "Now let’s think about which pages are likely to be important, in the sense of being valuable to a search engine user.\n",
    "\n",
    "One possible criterion for the importance of a page is the number of inbound links — an indication of popularity.\n",
    "\n",
    "By this measure, `m` and `j` are the most important pages, with 5 inbound links each.\n",
    "\n",
    "However, what if the pages linking to `m`, say, are not themselves important?\n",
    "\n",
    "Thinking this way, it seems appropriate to weight the inbound nodes by relative importance.\n",
    "\n",
    "The PageRank algorithm does precisely this.\n",
    "\n",
    "A slightly simplified presentation that captures the basic idea is as follows.\n",
    "\n",
    "Letting $ j $ be (the integer index of) a typical page and $ r_j $ be its ranking, we set\n",
    "\n",
    "$$\n",
    "r_j = \\sum_{i \\in L_j} \\frac{r_i}{\\ell_i}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ \\ell_i $ is the total number of outbound links from $ i $  \n",
    "- $ L_j $ is the set of all pages $ i $ such that $ i $ has a link to $ j $  \n",
    "\n",
    "\n",
    "This is a measure of the number of inbound links, weighted by their own ranking (and normalized by $ 1 / \\ell_i $).\n",
    "\n",
    "There is, however, another interpretation, and it brings us back to Markov chains.\n",
    "\n",
    "Let $ P $ be the matrix given by $ P(i, j) = \\mathbf 1\\{i \\to j\\} / \\ell_i $ where $ \\mathbf 1\\{i \\to j\\} = 1 $ if $ i $ has a link to $ j $ and zero otherwise.\n",
    "\n",
    "The matrix $ P $ is a stochastic matrix provided that each page has at least one link.\n",
    "\n",
    "With this definition of $ P $ we have\n",
    "\n",
    "$$\n",
    "r_j\n",
    "= \\sum_{i \\in L_j} \\frac{r_i}{\\ell_i}\n",
    "= \\sum_{\\text{all } i} \\mathbf 1\\{i \\to j\\} \\frac{r_i}{\\ell_i}\n",
    "= \\sum_{\\text{all } i} P(i, j) r_i\n",
    "$$\n",
    "\n",
    "Writing $ r $ for the row vector of rankings, this becomes $ r = r P $.\n",
    "\n",
    "Hence $ r $ is the stationary distribution of the stochastic matrix $ P $.\n",
    "\n",
    "Let’s think of $ P(i, j) $ as the probability of “moving” from page $ i $ to page $ j $.\n",
    "\n",
    "The value $ P(i, j) $ has the interpretation\n",
    "\n",
    "- $ P(i, j) = 1/k $ if $ i $ has $ k $ outbound links and $ j $ is one of them  \n",
    "- $ P(i, j) = 0 $ if $ i $ has no direct link to $ j $  \n",
    "\n",
    "\n",
    "Thus, motion from page to page is that of a web surfer who moves from one page to another by randomly clicking on one of the links on that page.\n",
    "\n",
    "Here “random” means that each link is selected with equal probability.\n",
    "\n",
    "Since $ r $ is the stationary distribution of $ P $, assuming that the uniform ergodicity condition is valid, we [can interpret](#new-interp-sd) $ r_j $ as the fraction of time that a (very persistent) random surfer spends at page $ j $.\n",
    "\n",
    "Your exercise is to apply this ranking algorithm to the graph pictured above\n",
    "and return the list of pages ordered by rank.\n",
    "\n",
    "There is a total of 14 nodes (i.e., web pages), the first named `a` and the last named `n`.\n",
    "\n",
    "A typical line from the file has the form"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a611ac0",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```text\n",
    "d -> h;\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b38b99d",
   "metadata": {},
   "source": [
    "This should be interpreted as meaning that there exists a link from `d` to `h`.\n",
    "\n",
    "The data for this graph is shown below, and read into a file called `web_graph_data.txt` when the cell is executed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c4b1a20",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%file web_graph_data.txt\n",
    "a -> d;\n",
    "a -> f;\n",
    "b -> j;\n",
    "b -> k;\n",
    "b -> m;\n",
    "c -> c;\n",
    "c -> g;\n",
    "c -> j;\n",
    "c -> m;\n",
    "d -> f;\n",
    "d -> h;\n",
    "d -> k;\n",
    "e -> d;\n",
    "e -> h;\n",
    "e -> l;\n",
    "f -> a;\n",
    "f -> b;\n",
    "f -> j;\n",
    "f -> l;\n",
    "g -> b;\n",
    "g -> j;\n",
    "h -> d;\n",
    "h -> g;\n",
    "h -> l;\n",
    "h -> m;\n",
    "i -> g;\n",
    "i -> h;\n",
    "i -> n;\n",
    "j -> e;\n",
    "j -> i;\n",
    "j -> k;\n",
    "k -> n;\n",
    "l -> m;\n",
    "m -> g;\n",
    "n -> c;\n",
    "n -> j;\n",
    "n -> m;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea1f53bb",
   "metadata": {},
   "source": [
    "To parse this file and extract the relevant information, you can use [regular expressions](https://docs.python.org/3/library/re.html).\n",
    "\n",
    "The following code snippet provides a hint as to how you can go about this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd65aeac",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import re\n",
    "re.findall('\\w', 'x +++ y ****** z')  # \\w matches alphanumerics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b95e2af",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "re.findall('\\w', 'a ^^ b &&& $$ c')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b083e5fb",
   "metadata": {},
   "source": [
    "When you solve for the ranking, you will find that the highest ranked node is in fact `g`, while the lowest is `a`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec4b6ef",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 13.2](https://tools-techniques.quantecon.org/#fm_ex2)\n",
    "\n",
    "Here is one solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6681cb5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Return list of pages, ordered by rank\n",
    "\"\"\"\n",
    "import re\n",
    "from operator import itemgetter\n",
    "\n",
    "infile = 'web_graph_data.txt'\n",
    "alphabet = 'abcdefghijklmnopqrstuvwxyz'\n",
    "\n",
    "n = 14 # Total number of web pages (nodes)\n",
    "\n",
    "# Create a matrix Q indicating existence of links\n",
    "#  * Q[i, j] = 1 if there is a link from i to j\n",
    "#  * Q[i, j] = 0 otherwise\n",
    "Q = np.zeros((n, n), dtype=int)\n",
    "with open(infile) as f:\n",
    "    edges = f.readlines()\n",
    "for edge in edges:\n",
    "    from_node, to_node = re.findall('\\w', edge)\n",
    "    i, j = alphabet.index(from_node), alphabet.index(to_node)\n",
    "    Q[i, j] = 1\n",
    "# Create the corresponding Markov matrix P\n",
    "P = np.empty((n, n))\n",
    "for i in range(n):\n",
    "    P[i, :] = Q[i, :] / Q[i, :].sum()\n",
    "mc = MarkovChain(P)\n",
    "# Compute the stationary distribution r\n",
    "r = mc.stationary_distributions[0]\n",
    "ranked_pages = {alphabet[i] : r[i] for i in range(n)}\n",
    "# Print solution, sorted from highest to lowest rank\n",
    "print('Rankings\\n ***')\n",
    "for name, rank in sorted(ranked_pages.items(), key=itemgetter(1), reverse=1):\n",
    "    print(f'{name}: {rank:.4}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbb2239f",
   "metadata": {},
   "source": [
    "## Exercise 13.3\n",
    "\n",
    "In numerical work, it is sometimes convenient to replace a continuous model with a discrete one.\n",
    "\n",
    "In particular, Markov chains are routinely generated as discrete approximations to AR(1) processes of the form\n",
    "\n",
    "$$\n",
    "y_{t+1} = \\rho y_t + u_{t+1}\n",
    "$$\n",
    "\n",
    "Here $ {u_t} $ is assumed to be IID and $ N(0, \\sigma_u^2) $.\n",
    "\n",
    "The variance of the stationary probability distribution of $ \\{ y_t \\} $ is\n",
    "\n",
    "$$\n",
    "\\sigma_y^2 := \\frac{\\sigma_u^2}{1-\\rho^2}\n",
    "$$\n",
    "\n",
    "Tauchen’s method [[Tauchen, 1986](https://tools-techniques.quantecon.org/zreferences.html#id212)] is the most common method for approximating this continuous state process with a finite state Markov chain.\n",
    "\n",
    "A routine for this already exists in [QuantEcon.py](http://quantecon.org/quantecon-py) but let’s write our own version as an exercise.\n",
    "\n",
    "As a first step, we choose\n",
    "\n",
    "- $ n $, the number of states for the discrete approximation  \n",
    "- $ m $, an integer that parameterizes the width of the state space  \n",
    "\n",
    "\n",
    "Next, we create a state space $ \\{x_0, \\ldots, x_{n-1}\\} \\subset \\mathbb R $\n",
    "and a stochastic $ n \\times n $ matrix $ P $ such that\n",
    "\n",
    "- $ x_0 = - m \\, \\sigma_y $  \n",
    "- $ x_{n-1} = m \\, \\sigma_y $  \n",
    "- $ x_{i+1} = x_i + s $ where $ s = (x_{n-1} - x_0) / (n - 1) $  \n",
    "\n",
    "\n",
    "Let $ F $ be the cumulative distribution function of the normal distribution $ N(0, \\sigma_u^2) $.\n",
    "\n",
    "The values $ P(x_i, x_j) $ are computed to approximate the AR(1) process — omitting the derivation, the rules are as follows:\n",
    "\n",
    "1. If $ j = 0 $, then set  \n",
    "  $$\n",
    "  P(x_i, x_j) = P(x_i, x_0) = F(x_0-\\rho x_i + s/2)\n",
    "  $$\n",
    "1. If $ j = n-1 $, then set  \n",
    "  $$\n",
    "  P(x_i, x_j) = P(x_i, x_{n-1}) = 1 - F(x_{n-1} - \\rho x_i - s/2)\n",
    "  $$\n",
    "1. Otherwise, set  \n",
    "  $$\n",
    "  P(x_i, x_j) = F(x_j - \\rho x_i + s/2) - F(x_j - \\rho x_i - s/2)\n",
    "  $$\n",
    "\n",
    "\n",
    "The exercise is to write a function `approx_markov(rho, sigma_u, m=3, n=7)` that returns\n",
    "$ \\{x_0, \\ldots, x_{n-1}\\} \\subset \\mathbb R $ and $ n \\times n $ matrix\n",
    "$ P $ as described above.\n",
    "\n",
    "- Even better, write a function that returns an instance of [QuantEcon.py’s](http://quantecon.org/quantecon-py) MarkovChain class.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc0734c8",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 13.3](https://tools-techniques.quantecon.org/#fm_ex3)\n",
    "\n",
    "A solution from the [QuantEcon.py](http://quantecon.org/quantecon-py) library\n",
    "can be found [here](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/markov/approximation.py).\n",
    "\n",
    "<p><a id=pm href=#pm-link><strong>[1]</strong></a> Hint: First show that if $ P $ and $ Q $ are stochastic matrices then so is their product — to check the row sums, try post multiplying by a column vector of ones.  Finally, argue that $ P^n $ is a stochastic matrix using induction."
   ]
  }
 ],
 "metadata": {
  "date": 1718249790.8235233,
  "filename": "finite_markov.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Finite Markov Chains"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}