muth_kalman.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='muth-kalman'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reverse Engineering a la Muth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Reverse Engineering a la Muth](#Reverse-Engineering-a-la-Muth)  \n",
    "  - [Friedman (1956) and Muth (1960)](#Friedman-%281956%29-and-Muth-%281960%29)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Co-author: Chase Coleman**\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture uses the quantecon library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll also need the following imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "\n",
    "from quantecon import Kalman\n",
    "from quantecon import LinearStateSpace\n",
    "from scipy.stats import norm\n",
    "np.set_printoptions(linewidth=120, precision=4, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lecture uses the Kalman filter to reformulate John F. Muth’s first\n",
    "paper [[Mut60]](https://lectures.quantecon.org/zreferences.html#muth1960) about rational expectations.\n",
    "\n",
    "Muth used *classical* prediction methods to reverse engineer a\n",
    "stochastic process that renders optimal Milton Friedman’s [[Fri56]](https://lectures.quantecon.org/zreferences.html#friedman1956) “adaptive\n",
    "expectations” scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Friedman (1956) and Muth (1960)\n",
    "\n",
    "Milton Friedman [[Fri56]](https://lectures.quantecon.org/zreferences.html#friedman1956) (1956) posited that\n",
    "consumer’s forecast their future disposable income with the adaptive\n",
    "expectations scheme\n",
    "\n",
    "\n",
    "<a id='equation-expectations'></a>\n",
    "$$\n",
    "y_{t+i,t}^* = K \\sum_{j=0}^\\infty (1 - K)^j y_{t-j} \\tag{1}\n",
    "$$\n",
    "\n",
    "where $ K \\in (0,1) $ and $ y_{t+i,t}^* $ is a forecast of\n",
    "future $ y $ over horizon $ i $.\n",
    "\n",
    "Milton Friedman justified the **exponential smoothing** forecasting\n",
    "scheme [(1)](#equation-expectations) informally, noting  that it seemed a plausible way to use\n",
    "past income to forecast future income.\n",
    "\n",
    "In his first paper about rational expectations, John F. Muth [[Mut60]](https://lectures.quantecon.org/zreferences.html#muth1960)\n",
    "reverse-engineered a univariate stochastic process\n",
    "$ \\{y_t\\}_{t=- \\infty}^\\infty $ for which Milton Friedman’s adaptive\n",
    "expectations scheme gives linear least forecasts of $ y_{t+j} $ for\n",
    "any horizon $ i $.\n",
    "\n",
    "Muth sought a setting and a sense in which Friedman’s forecasting scheme\n",
    "is optimal.\n",
    "\n",
    "That is, Muth asked for what optimal forecasting **question** is Milton\n",
    "Friedman’s adaptive expectation scheme the **answer**.\n",
    "\n",
    "Muth (1960) used classical prediction methods based on lag-operators and\n",
    "$ z $-transforms to find the answer to his question.\n",
    "\n",
    "Please see lectures [Classical Control with Linear Algebra](https://lectures.quantecon.org/lu_tricks.html) and\n",
    "[Classical Filtering and Prediction with Linear Algebra](https://lectures.quantecon.org/classical_filtering.html) for an introduction to the classical\n",
    "tools that Muth used.\n",
    "\n",
    "Rather than using those classical tools, in this lecture we apply the\n",
    "Kalman filter to express the heart of Muth’s analysis concisely.\n",
    "\n",
    "The lecture [First Look at Kalman Filter](https://lectures.quantecon.org/kalman.html) describes the Kalman filter.\n",
    "\n",
    "We’ll use limiting versions of the Kalman filter corresponding to what are called **stationary values** in that lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Process for Which Adaptive Expectations are Optimal\n",
    "\n",
    "Suppose that an observable $ y_t $ is the sum of an unobserved\n",
    "random walk $ x_t $ and an IID shock $ \\epsilon_{2,t} $:\n",
    "\n",
    "\n",
    "<a id='equation-state-space'></a>\n",
    "$$\n",
    "\\begin{aligned} x_{t+1} & = x_t + \\sigma_x \\epsilon_{1,t+1} \\cr\n",
    "          y_t & = x_t + \\sigma_y \\epsilon_{2,t} \\end{aligned} \\tag{2}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} \\epsilon_{1,t+1} \\cr \\epsilon_{2,t} \\end{bmatrix} \\sim {\\mathcal N} (0, I)\n",
    "$$\n",
    "\n",
    "is an IID process.\n",
    "\n",
    "**Note:** A property of the state-space representation [(2)](#equation-state-space) is that in\n",
    "general neither $ \\epsilon_{1,t} $ nor $ \\epsilon_{2,t} $ is in\n",
    "the space spanned by square-summable linear combinations of\n",
    "$ y_t, y_{t-1}, \\ldots $.\n",
    "\n",
    "In general\n",
    "$ \\begin{bmatrix} \\epsilon_{1,t} \\cr \\epsilon_{2t} \\end{bmatrix} $\n",
    "has more information about future $ y_{t+j} $’s than is contained\n",
    "in $ y_t, y_{t-1}, \\ldots $.\n",
    "\n",
    "We can use the asymptotic or stationary values of the Kalman gain and\n",
    "the one-step-ahead conditional state covariance matrix to compute a\n",
    "time-invariant *innovations representation*\n",
    "\n",
    "\n",
    "<a id='equation-innovations'></a>\n",
    "$$\n",
    "\\begin{aligned} \\hat x_{t+1} & = \\hat x_t + K a_t  \\cr\n",
    "           y_t & = \\hat x_t + a_t \\end{aligned} \\tag{3}\n",
    "$$\n",
    "\n",
    "where $ \\hat x_t = E [x_t | y_{t-1}, y_{t-2}, \\ldots ] $ and\n",
    "$ a_t = y_t - E[y_t |y_{t-1}, y_{t-2}, \\ldots ] $.\n",
    "\n",
    "**Note:** A key property about an *innovations representation* is that\n",
    "$ a_t $ is in the space spanned by square summable linear\n",
    "combinations of $ y_t, y_{t-1}, \\ldots $.\n",
    "\n",
    "For more ramifications of this property, see the lectures  [Shock Non-Invertibility](https://lectures.quantecon.org/hs_invertibility_example.html)  and\n",
    "[Recursive Models of Dynamic Linear Economies](https://lectures.quantecon.org/hs_recursive_models.html).\n",
    "\n",
    "Later we’ll stack these state-space systems [(2)](#equation-state-space) and [(3)](#equation-innovations) to display some\n",
    "classic findings of Muth.\n",
    "\n",
    "But first, let’s create an instance of the state-space system [(2)](#equation-state-space) then\n",
    "apply the quantecon `Kalman` class, then uses it to construct the associated “innovations representation”"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Make some parameter choices\n",
    "# sigx/sigy are state noise std err and measurement noise std err\n",
    "μ_0, σ_x, σ_y = 10, 1, 5\n",
    "\n",
    "# Create a LinearStateSpace object\n",
    "A, C, G, H = 1, σ_x, 1, σ_y\n",
    "ss = LinearStateSpace(A, C, G, H, mu_0=μ_0)\n",
    "\n",
    "# Set prior and initialize the Kalman type\n",
    "x_hat_0, Σ_0 = 10, 1\n",
    "kmuth = Kalman(ss, x_hat_0, Σ_0)\n",
    "\n",
    "# Computes stationary values which we need for the innovation\n",
    "# representation\n",
    "S1, K1 = kmuth.stationary_values()\n",
    "\n",
    "# Form innovation representation state-space\n",
    "Ak, Ck, Gk, Hk = A, K1, G, 1\n",
    "\n",
    "ssk = LinearStateSpace(Ak, Ck, Gk, Hk, mu_0=x_hat_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some Useful State-Space Math\n",
    "\n",
    "Now we want to map the time-invariant innovations representation [(3)](#equation-innovations) and\n",
    "the original state-space system [(2)](#equation-state-space) into a convenient form for deducing\n",
    "the impulse responses from the original shocks to the $ x_t $ and\n",
    "$ \\hat x_t $.\n",
    "\n",
    "Putting both of these representations into a single state-space system\n",
    "is yet another application of the insight that “finding the state is an\n",
    "art”.\n",
    "\n",
    "We’ll define a state vector and appropriate state-space matrices that\n",
    "allow us to represent both systems in one fell swoop.\n",
    "\n",
    "Note that\n",
    "\n",
    "$$\n",
    "a_t = x_t + \\sigma_y \\epsilon_{2,t} - \\hat x_t\n",
    "$$\n",
    "\n",
    "so that\n",
    "\n",
    "$$\n",
    "\\begin{aligned} \\hat x_{t+1} & = \\hat x_t + K (x_t + \\sigma_y \\epsilon_{2,t} - \\hat x_t) \\cr\n",
    "       & = (1-K) \\hat x_t + K x_t + K \\sigma_y \\epsilon_{2,t} \\end{aligned}\n",
    "$$\n",
    "\n",
    "The stacked system\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} x_{t+1} \\cr \\hat x_{t+1} \\cr \\epsilon_{2,t+1} \\end{bmatrix} =\n",
    "\\begin{bmatrix} 1 & 0 & 0 \\cr K & (1-K) & K \\sigma_y \\cr 0 & 0 & 0 \\end{bmatrix}\n",
    "\\begin{bmatrix} x_{t} \\cr \\hat x_t \\cr \\epsilon_{2,t} \\end{bmatrix}+\n",
    "\\begin{bmatrix} \\sigma_x & 0 \\cr 0 & 0 \\cr 0 & 1 \\end{bmatrix}\n",
    "\\begin{bmatrix} \\epsilon_{1,t+1} \\cr \\epsilon_{2,t+1} \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} y_t \\cr a_t \\end{bmatrix} = \\begin{bmatrix} 1 & 0 & \\sigma_y \\cr\n",
    "                                      1 & -1 & \\sigma_y \\end{bmatrix}  \\begin{bmatrix} x_{t} \\cr \\hat x_t \\cr \\epsilon_{2,t} \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "is a state-space system that tells us how the shocks\n",
    "$ \\begin{bmatrix} \\epsilon_{1,t+1} \\cr \\epsilon_{2,t+1} \\end{bmatrix} $\n",
    "affect states $ \\hat x_{t+1}, x_t $, the observable $ y_t $, and\n",
    "the innovation $ a_t $.\n",
    "\n",
    "With this tool at our disposal, let’s form the composite system and\n",
    "simulate it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Create grand state-space for y_t, a_t as observed vars -- Use\n",
    "# stacking trick above\n",
    "Af = np.array([[ 1,      0,        0],\n",
    "               [K1, 1 - K1, K1 * σ_y],\n",
    "               [ 0,      0,        0]])\n",
    "Cf = np.array([[σ_x,        0],\n",
    "               [  0, K1 * σ_y],\n",
    "               [  0,        1]])\n",
    "Gf = np.array([[1,  0, σ_y],\n",
    "               [1, -1, σ_y]])\n",
    "\n",
    "μ_true, μ_prior = 10, 10\n",
    "μ_f = np.array([μ_true, μ_prior, 0]).reshape(3, 1)\n",
    "\n",
    "# Create the state-space\n",
    "ssf = LinearStateSpace(Af, Cf, Gf, mu_0=μ_f)\n",
    "\n",
    "# Draw observations of y from the state-space model\n",
    "N = 50\n",
    "xf, yf = ssf.simulate(N)\n",
    "\n",
    "print(f\"Kalman gain = {K1}\")\n",
    "print(f\"Conditional variance = {S1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have simulated our joint system, we have $ x_t $,\n",
    "$ \\hat{x_t} $, and $ y_t $.\n",
    "\n",
    "We can now investigate how these\n",
    "variables are related by plotting some key objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimates of Unobservables\n",
    "\n",
    "First, let’s plot the hidden state $ x_t $ and the filtered version\n",
    "$ \\hat x_t $ that is linear-least squares projection of $ x_t $\n",
    "on the history $ y_{t-1}, y_{t-2}, \\ldots $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(xf[0, :], label=\"$x_t$\")\n",
    "ax.plot(xf[1, :], label=\"Filtered $x_t$\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_title(r\"$x$ vs $\\hat{x}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how $ x_t $ and $ \\hat{x_t} $ differ.\n",
    "\n",
    "For Friedman, $ \\hat x_t $ and not $ x_t $ is the consumer’s\n",
    "idea about her/his *permanent income*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relation between Unobservable and Observable\n",
    "\n",
    "Now let’s plot $ x_t $ and $ y_t $.\n",
    "\n",
    "Recall that $ y_t $ is just $ x_t $ plus white noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(yf[0, :], label=\"y\")\n",
    "ax.plot(xf[0, :], label=\"x\")\n",
    "ax.legend()\n",
    "ax.set_title(r\"$x$ and $y$\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see above that $ y $ seems to look like white noise around the\n",
    "values of $ x $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Innovations\n",
    "\n",
    "Recall that we wrote down the innovation representation that depended on\n",
    "$ a_t $. We now plot the innovations $ \\{a_t\\} $:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(yf[1, :], label=\"a\")\n",
    "ax.legend()\n",
    "ax.set_title(r\"Innovation $a_t$\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MA and AR Representations\n",
    "\n",
    "Now we shall extract from the `Kalman` instance `kmuth` coefficients of\n",
    "\n",
    "- a fundamental moving average representation that represents\n",
    "  $ y_t $ as a one-sided moving sum of current and past\n",
    "  $ a_t $s that are square summable linear combinations of $ y_t, y_{t-1}, \\ldots $.  \n",
    "- a univariate autoregression representation that depicts the\n",
    "  coefficients in a linear least square projection of $ y_t $ on\n",
    "  the semi-infinite history $ y_{t-1}, y_{t-2}, \\ldots $.  \n",
    "\n",
    "\n",
    "Then we’ll plot each of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Kalman Methods for MA and VAR\n",
    "coefs_ma = kmuth.stationary_coefficients(5, \"ma\")\n",
    "coefs_var = kmuth.stationary_coefficients(5, \"var\")\n",
    "\n",
    "# Coefficients come in a list of arrays, but we\n",
    "# want to plot them and so need to stack into an array\n",
    "coefs_ma_array = np.vstack(coefs_ma)\n",
    "coefs_var_array = np.vstack(coefs_var)\n",
    "\n",
    "fig, ax = plt.subplots(2)\n",
    "ax[0].plot(coefs_ma_array, label=\"MA\")\n",
    "ax[0].legend()\n",
    "ax[1].plot(coefs_var_array, label=\"VAR\")\n",
    "ax[1].legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **moving average** coefficients in the top panel show tell-tale\n",
    "signs of $ y_t $ being a process whose first difference is a first-order\n",
    "autoregression.\n",
    "\n",
    "The **autoregressive coefficients** decline geometrically with decay\n",
    "rate $ (1-K) $.\n",
    "\n",
    "These are exactly the target outcomes that Muth (1960) aimed to reverse\n",
    "engineer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "print(f'decay parameter 1 - K1 = {1 - K1}')"
   ]
  }
 ],
 "metadata": {
  "filename": "muth_kalman.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Reverse Engineering a la Muth"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}