multiplicative_functionals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='multiplicative-functionals'></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>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiplicative Functionals\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Multiplicative Functionals](#Multiplicative-Functionals)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [A Log-Likelihood Process](#A-Log-Likelihood-Process)  \n",
    "  - [Benefits from Reduced Aggregate Fluctuations](#Benefits-from-Reduced-Aggregate-Fluctuations)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Co-authors: Chase Coleman and Balint Szoke**\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture is a sequel to the [lecture on additive functionals](https://lectures.quantecon.org/additive_functionals.html).\n",
    "\n",
    "That lecture\n",
    "\n",
    "1. defined a special class of **additive functionals** driven by a first-order vector VAR  \n",
    "1. by taking the exponential of that additive functional, created an associated **multiplicative functional**  \n",
    "\n",
    "\n",
    "This lecture uses this special class to create and analyze two examples\n",
    "\n",
    "- A  **log-likelihood process**, an object at the foundation of both frequentist and Bayesian approaches to statistical inference.  \n",
    "- A version of Robert E. Lucas’s [[LJ03]](https://lectures.quantecon.org/zreferences.html#lucas2003macroeconomic) and Thomas Tallarini’s [[Tal00]](https://lectures.quantecon.org/zreferences.html#tall2000) approaches to measuring the benefits of moderating aggregate fluctuations.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Log-Likelihood Process\n",
    "\n",
    "Consider a vector of additive functionals $ \\{y_t\\}_{t=0}^\\infty $\n",
    "described by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    x_{t+1} & = A x_t + B z_{t+1}\n",
    "    \\\\\n",
    "    y_{t+1} - y_t & = D x_{t} + F z_{t+1},\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $ A $ is a stable matrix, $ \\{z_{t+1}\\}_{t=0}^\\infty $ is\n",
    "an IID sequence of $ {\\cal N}(0,I) $ random vectors, $ F $ is\n",
    "nonsingular, and $ x_0 $ and $ y_0 $ are vectors of known\n",
    "numbers.\n",
    "\n",
    "Evidently,\n",
    "\n",
    "$$\n",
    "x_{t+1} = \\left(A - B F^{-1}D \\right)x_t +\n",
    "B F^{-1} \\left(y_{t+1} - y_t \\right),\n",
    "$$\n",
    "\n",
    "so that $ x_{t+1} $ can be constructed from observations on\n",
    "$ \\{y_{s}\\}_{s=0}^{t+1} $ and $ x_0 $.\n",
    "\n",
    "The distribution of $ y_{t+1} - y_t $ conditional on $ x_t $ is normal with mean $ Dx_t $ and nonsingular covariance matrix $ FF' $.\n",
    "\n",
    "Let $ \\theta $ denote the vector of free parameters of the model.\n",
    "\n",
    "These parameters pin down the elements of $ A, B, D, F $.\n",
    "\n",
    "The **log-likelihood function** of $ \\{y_s\\}_{s=1}^t $ is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\log L_{t}(\\theta)  =\n",
    "    & - {\\frac 1 2} \\sum_{j=1}^{t} (y_{j} - y_{j-1} -\n",
    "         D x_{j-1})'(FF')^{-1}(y_{j} - y_{j-1} - D x_{j-1})\n",
    "    \\\\\n",
    "    & - {\\frac t 2} \\log \\det (FF') - {\\frac {k t} 2} \\log( 2 \\pi)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Let’s consider the case of a scalar process in which $ A, B, D, F $ are scalars and $ z_{t+1} $ is a scalar stochastic process.\n",
    "\n",
    "We let $ \\theta_o $ denote the “true” values of $ \\theta $, meaning the values that generate the data.\n",
    "\n",
    "For the purposes of this exercise,  set $ \\theta_o = (A, B, D, F) = (0.8, 1, 0.5, 0.2) $.\n",
    "\n",
    "Set $ x_0 = y_0 = 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating Sample Paths\n",
    "\n",
    "Let’s write a program to simulate sample paths of $ \\{ x_t, y_{t} \\}_{t=0}^{\\infty} $.\n",
    "\n",
    "We’ll do this by formulating the additive functional as a linear state space model and putting the [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to work.\n",
    "\n",
    "> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "\n",
    "@authors: Chase Coleman, Balint Skoze, Tom Sargent\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.linalg as la\n",
    "import quantecon as qe\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import lognorm\n",
    "\n",
    "\n",
    "class AMF_LSS_VAR:\n",
    "    \"\"\"\n",
    "    This class is written to transform a scalar additive functional\n",
    "    into a linear state space system.\n",
    "    \"\"\"\n",
    "    def __init__(self, A, B, D, F=0.0, ν=0.0):\n",
    "        # Unpack required elements\n",
    "        self.A, self.B, self.D, self.F, self.ν = A, B, D, F, ν\n",
    " \n",
    "        # Create space for additive decomposition\n",
    "        self.add_decomp = None\n",
    "        self.mult_decomp = None\n",
    " \n",
    "        # Construct BIG state space representation\n",
    "        self.lss = self.construct_ss()\n",
    " \n",
    "    def construct_ss(self):\n",
    "        \"\"\"\n",
    "        This creates the state space representation that can be passed\n",
    "        into the quantecon LSS class.\n",
    "        \"\"\"\n",
    "        # Pull out useful info\n",
    "        A, B, D, F, ν = self.A, self.B, self.D, self.F, self.ν\n",
    "        nx, nk, nm = 1, 1, 1\n",
    "        if self.add_decomp:\n",
    "            ν, H, g = self.add_decomp\n",
    "        else:\n",
    "            ν, H, g = self.additive_decomp()\n",
    " \n",
    "        # Build A matrix for LSS\n",
    "        # Order of states is: [1, t, xt, yt, mt]\n",
    "        A1 = np.hstack([1, 0, 0, 0, 0])       # Transition for 1\n",
    "        A2 = np.hstack([1, 1, 0, 0, 0])       # Transition for t\n",
    "        A3 = np.hstack([0, 0, A, 0, 0])       # Transition for x_{t+1}\n",
    "        A4 = np.hstack([ν, 0, D, 1, 0])       # Transition for y_{t+1}\n",
    "        A5 = np.hstack([0, 0, 0, 0, 1])       # Transition for m_{t+1}\n",
    "        Abar = np.vstack([A1, A2, A3, A4, A5])\n",
    " \n",
    "        # Build B matrix for LSS\n",
    "        Bbar = np.vstack([0, 0, B, F, H])\n",
    " \n",
    "        # Build G matrix for LSS\n",
    "        # Order of observation is: [xt, yt, mt, st, tt]\n",
    "        G1 = np.hstack([0, 0, 1, 0, 0])               # Selector for x_{t}\n",
    "        G2 = np.hstack([0, 0, 0, 1, 0])               # Selector for y_{t}\n",
    "        G3 = np.hstack([0, 0, 0, 0, 1])               # Selector for martingale\n",
    "        G4 = np.hstack([0, 0, -g, 0, 0])              # Selector for stationary\n",
    "        G5 = np.hstack([0, ν, 0, 0, 0])               # Selector for trend\n",
    "        Gbar = np.vstack([G1, G2, G3, G4, G5])\n",
    " \n",
    "        # Build H matrix for LSS\n",
    "        Hbar = np.zeros((1, 1))\n",
    " \n",
    "        # Build LSS type\n",
    "        x0 = np.hstack([1, 0, 0, 0, 0])\n",
    "        S0 = np.zeros((5, 5))\n",
    "        lss = qe.lss.LinearStateSpace(Abar, Bbar, Gbar, Hbar, mu_0=x0, Sigma_0=S0)\n",
    " \n",
    "        return lss\n",
    " \n",
    "    def additive_decomp(self):\n",
    "        \"\"\"\n",
    "        Return values for the martingale decomposition (Proposition 4.3.3.)\n",
    "            - ν         : unconditional mean difference in Y\n",
    "            - H         : coefficient for the (linear) martingale component (kappa_a)\n",
    "            - g         : coefficient for the stationary component g(x)\n",
    "            - Y_0       : it should be the function of X_0 (for now set it to 0.0)\n",
    "        \"\"\"\n",
    "        A_res = 1 / (1 - self.A)\n",
    "        g = self.D * A_res\n",
    "        H = self.F + self.D * A_res * self.B\n",
    " \n",
    "        return self.ν, H, g\n",
    " \n",
    "    def multiplicative_decomp(self):\n",
    "        \"\"\"\n",
    "        Return values for the multiplicative decomposition (Example 5.4.4.)\n",
    "            - ν_tilde  : eigenvalue\n",
    "            - H        : vector for the Jensen term\n",
    "        \"\"\"\n",
    "        ν, H, g = self.additive_decomp()\n",
    "        ν_tilde = ν + (.5) * H**2\n",
    " \n",
    "        return ν_tilde, H, g\n",
    " \n",
    "    def loglikelihood_path(self, x, y):\n",
    "        A, B, D, F = self.A, self.B, self.D, self.F\n",
    "        T = y.T.size\n",
    "        FF = F**2\n",
    "        FFinv = 1 / FF\n",
    "        temp = y[1:] - y[:-1] - D * x[:-1]\n",
    "        obs = temp * FFinv * temp\n",
    "        obssum = np.cumsum(obs)\n",
    "        scalar = (np.log(FF) + np.log(2 * np.pi)) * np.arange(1, T)\n",
    " \n",
    "        return (-0.5) * (obssum + scalar)\n",
    " \n",
    "    def loglikelihood(self, x, y):\n",
    "        llh = self.loglikelihood_path(x, y)\n",
    " \n",
    "        return llh[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The heavy lifting is done inside the AMF_LSS_VAR class.\n",
    "\n",
    "The following code adds some simple functions that make it straightforward to generate sample paths from an instance of AMF_LSS_VAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def simulate_xy(amf, T):\n",
    "    \"Simulate individual paths.\"\n",
    "    foo, bar = amf.lss.simulate(T)\n",
    "    x = bar[0, :]\n",
    "    y = bar[1, :]\n",
    "\n",
    "    return x, y\n",
    "\n",
    "def simulate_paths(amf, T=150, I=5000):\n",
    "    \"Simulate multiple independent paths.\"\n",
    "\n",
    "    # Allocate space\n",
    "    storeX = np.empty((I, T))\n",
    "    storeY = np.empty((I, T))\n",
    "\n",
    "    for i in range(I):\n",
    "        # Do specific simulation\n",
    "        x, y = simulate_xy(amf, T)\n",
    "\n",
    "        # Fill in our storage matrices\n",
    "        storeX[i, :] = x\n",
    "        storeY[i, :] = y\n",
    "\n",
    "    return storeX, storeY\n",
    "\n",
    "def population_means(amf, T=150):\n",
    "    # Allocate Space\n",
    "    xmean = np.empty(T)\n",
    "    ymean = np.empty(T)\n",
    "\n",
    "    # Pull out moment generator\n",
    "    moment_generator = amf.lss.moment_sequence()\n",
    "\n",
    "    for tt in range (T):\n",
    "        tmoms = next(moment_generator)\n",
    "        ymeans = tmoms[1]\n",
    "        xmean[tt] = ymeans[0]\n",
    "        ymean[tt] = ymeans[1]\n",
    "\n",
    "    return xmean, ymean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have these functions in our tool kit, let’s apply them to run some\n",
    "simulations.\n",
    "\n",
    "In particular, let’s use our program to generate $ I = 5000 $ sample paths of length $ T = 150 $, labeled $ \\{ x_{t}^i, y_{t}^i \\}_{t=0}^\\infty $ for $ i = 1, ..., I $.\n",
    "\n",
    "Then we compute averages of $ \\frac{1}{I} \\sum_i x_t^i $ and $ \\frac{1}{I} \\sum_i y_t^i $ across the sample paths and compare them with the population means of $ x_t $ and $ y_t $.\n",
    "\n",
    "Here goes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A, B, D, F = [0.8, 1.0, 0.5, 0.2]\n",
    "amf = AMF_LSS_VAR(A, B, D, F=F)\n",
    "\n",
    "T = 150\n",
    "I = 5000\n",
    "\n",
    "# Simulate and compute sample means\n",
    "Xit, Yit = simulate_paths(amf, T, I)\n",
    "Xmean_t = np.mean(Xit, 0)\n",
    "Ymean_t = np.mean(Yit, 0)\n",
    "\n",
    "# Compute population means\n",
    "Xmean_pop, Ymean_pop = population_means(amf, T)\n",
    "\n",
    "# Plot sample means vs population means\n",
    "fig, ax = plt.subplots(2, figsize=(14, 8))\n",
    "\n",
    "ax[0].plot(Xmean_t, label=r'$\\frac{1}{I}\\sum_i x_t^i$', color=\"b\")\n",
    "ax[0].plot(Xmean_pop, label='$\\mathbb{E} x_t$', color=\"k\")\n",
    "ax[0].set_title('$x_t$')\n",
    "ax[0].set_xlim((0, T))\n",
    "ax[0].legend(loc=0)\n",
    "\n",
    "ax[1].plot(Ymean_t, label=r'$\\frac{1}{I}\\sum_i y_t^i$', color=\"b\")\n",
    "ax[1].plot(Ymean_pop, label='$\\mathbb{E} y_t$', color=\"k\")\n",
    "ax[1].set_title('$y_t$')\n",
    "ax[1].set_xlim((0, T))\n",
    "ax[1].legend(loc=0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating Log-likelihoods\n",
    "\n",
    "Our next aim is to write a program to simulate $ \\{\\log L_t \\mid \\theta_o\\}_{t=1}^T $.\n",
    "\n",
    "We want as inputs to this program the *same* sample paths $ \\{x_t^i, y_t^i\\}_{t=0}^T $ that we  have already computed.\n",
    "\n",
    "We now want to simulate $ I = 5000 $ paths of $ \\{\\log L_t^i  \\mid \\theta_o\\}_{t=1}^T $.\n",
    "\n",
    "- For each path, we compute $ \\log L_T^i / T $.  \n",
    "- We also compute $ \\frac{1}{I} \\sum_{i=1}^I \\log L_T^i / T $.  \n",
    "\n",
    "\n",
    "Then we to compare these objects.\n",
    "\n",
    "Below we plot the histogram of $ \\log L_T^i / T $ for realizations $ i = 1, \\ldots, 5000 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def simulate_likelihood(amf, Xit, Yit):\n",
    "    # Get size\n",
    "    I, T = Xit.shape\n",
    "\n",
    "    # Allocate space\n",
    "    LLit = np.empty((I, T-1))\n",
    "\n",
    "    for i in range(I):\n",
    "        LLit[i, :] = amf.loglikelihood_path(Xit[i, :], Yit[i, :])\n",
    "\n",
    "    return LLit\n",
    "\n",
    "# Get likelihood from each path x^{i}, Y^{i}\n",
    "LLit = simulate_likelihood(amf, Xit, Yit)\n",
    "\n",
    "LLT = 1/T * LLit[:, -1]\n",
    "LLmean_t = np.mean(LLT)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.hist(LLT)\n",
    "ax.vlines(LLmean_t, ymin=0, ymax=I//3, color=\"k\", linestyle=\"--\", alpha=0.6)\n",
    "plt.title(r\"Distribution of $\\frac{1}{T} \\log L_{T}  \\mid \\theta_0$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the log-likelihood is almost always nonnegative, implying that $ L_t $ is typically bigger than 1.\n",
    "\n",
    "Recall that the likelihood function is a pdf (probability density function) and **not** a probability measure, so it can take values larger than 1.\n",
    "\n",
    "In the current case, the conditional variance of $ \\Delta y_{t+1} $, which equals  $ FF^T=0.04 $, is so small that the maximum value of the pdf is 2 (see the figure below).\n",
    "\n",
    "This implies that approximately $ 75\\% $ of the time (a bit more than one sigma deviation),  we should expect the **increment** of the log-likelihood to be nonnegative.\n",
    "\n",
    "Let’s see this in a simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "normdist = sp.stats.norm(0, F)\n",
    "mult = 1.175\n",
    "print(f'The pdf at +/- {mult} sigma takes the value:  {normdist.pdf(mult * F)}')\n",
    "print(f'Probability of dL being larger than 1 is approx: {normdist.cdf(mult * F) - normdist.cdf(-mult * F)}')\n",
    "\n",
    "# Compare this to the sample analogue:\n",
    "L_increment = LLit[:, 1:] - LLit[:, :-1]\n",
    "r, c = L_increment.shape\n",
    "frac_nonegative = np.sum(L_increment >= 0) / (c * r)\n",
    "print(f'Fraction of dlogL being nonnegative in the sample is: {frac_nonegative}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s also plot the conditional pdf of $ \\Delta y_{t+1} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "xgrid = np.linspace(-1, 1, 100)\n",
    "plt.plot(xgrid, normdist.pdf(xgrid))\n",
    "plt.title('Conditional pdf $f(\\Delta y_{t+1} \\mid x_t)$')\n",
    "print(f'The pdf at +/- one sigma takes the value: {normdist.pdf(F)}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An Alternative Parameter Vector\n",
    "\n",
    "Now consider alternative parameter vector $ \\theta_1 = [A, B, D, F] = [0.9, 1.0, 0.55, 0.25] $.\n",
    "\n",
    "We want to compute $ \\{\\log L_t \\mid \\theta_1\\}_{t=1}^T $.\n",
    "\n",
    "The $ x_t, y_t $ inputs to this program should be exactly the **same** sample paths $ \\{x_t^i, y_t^i\\}_{t=0}^T $ that we computed above.\n",
    "\n",
    "This is because we want to generate data under the $ \\theta_o $ probability model but evaluate the likelihood under the $ \\theta_1 $ model.\n",
    "\n",
    "So our task is to use our program to simulate $ I = 5000 $ paths of $ \\{\\log L_t^i  \\mid \\theta_1\\}_{t=1}^T $\n",
    "\n",
    "- For each path, compute $ \\frac{1}{T} \\log L_T^i $.  \n",
    "- Then compute $ \\frac{1}{I}\\sum_{i=1}^I \\frac{1}{T} \\log L_T^i $.  \n",
    "\n",
    "\n",
    "We want to compare these objects with each other and with the analogous objects that we computed above.\n",
    "\n",
    "Then we want to interpret outcomes.\n",
    "\n",
    "A function that we constructed can  handle these tasks.\n",
    "\n",
    "The only innovation is that we must create an alternative model to feed in.\n",
    "\n",
    "We will creatively call the new model `amf2`.\n",
    "\n",
    "We make three graphs\n",
    "\n",
    "- the first sets the stage by repeating an earlier graph  \n",
    "- the second contains two histograms of values of  log-likelihoods of the two models  over the period $ T $  \n",
    "- the third compares likelihoods under the true and alternative models  \n",
    "\n",
    "\n",
    "Here’s the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Create the second (wrong) alternative model\n",
    "A2, B2, D2, F2 = [0.9, 1.0, 0.55, 0.25]   #  parameters for θ_1 closer to Θ_0\n",
    "amf2 = AMF_LSS_VAR(A2, B2, D2, F=F2)\n",
    "\n",
    "# Get likelihood from each path x^{i}, y^{i}\n",
    "LLit2 = simulate_likelihood(amf2, Xit, Yit)\n",
    "\n",
    "LLT2 = 1/(T-1) * LLit2[:, -1]\n",
    "LLmean_t2 = np.mean(LLT2)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.hist(LLT2)\n",
    "ax.vlines(LLmean_t2, ymin=0, ymax=1400, color=\"k\", linestyle=\"--\", alpha=0.6)\n",
    "\n",
    "plt.title(r\"Distribution of $\\frac{1}{T} \\log L_{T}  \\mid \\theta_1$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s see a histogram of the log-likelihoods under the true and the alternative model (same sample paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "plt.hist(LLT, bins=50, alpha=0.5, label='True', density=True)\n",
    "plt.hist(LLT2, bins=50, alpha=0.5, label='Alternative', density=True)\n",
    "plt.vlines(np.mean(LLT), 0, 10, color='k', linestyle=\"--\", linewidth= 4)\n",
    "plt.vlines(np.mean(LLT2), 0, 10, color='k', linestyle=\"--\", linewidth= 4)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we’ll plot the histogram of the difference in log-likelihood ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "LLT_diff = LLT - LLT2\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.hist(LLT_diff, bins=50)\n",
    "plt.title(r\"$\\frac{1}{T}\\left[\\log (L_T^i  \\mid \\theta_0) - \\log (L_T^i \\mid \\theta_1)\\right]$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation\n",
    "\n",
    "These histograms of  log-likelihood ratios illustrate  important features of **likelihood ratio tests** as tools for discriminating between statistical models.\n",
    "\n",
    "- The log-likelihood is higher on average under the true model – obviously a very useful property.  \n",
    "- Nevertheless, for a positive fraction of realizations, the log-likelihood is higher for the incorrect than for the true model  \n",
    "\n",
    "\n",
    "> - in these instances, a likelihood ratio test mistakenly selects the wrong model  \n",
    "\n",
    "\n",
    "\n",
    "- These mechanics underlie the statistical theory of **mistake probabilities** associated with model selection tests based on  likelihood ratio.  \n",
    "\n",
    "\n",
    "(In a subsequent lecture, we’ll use some of the code prepared in this lecture to illustrate mistake probabilities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benefits from Reduced Aggregate Fluctuations\n",
    "\n",
    "Now let’s turn to a new example of multiplicative functionals.\n",
    "\n",
    "This example illustrates  ideas in the literature on\n",
    "\n",
    "- **long-run risk** in the consumption based asset pricing literature (e.g., [[BY04]](https://lectures.quantecon.org/zreferences.html#bansal2004risks), [[HHL08]](https://lectures.quantecon.org/zreferences.html#hansen2008consumption), [[Han07]](https://lectures.quantecon.org/zreferences.html#hansen2007beliefs))  \n",
    "- **benefits of eliminating aggregate fluctuations** in representative agent macro models (e.g., [[Tal00]](https://lectures.quantecon.org/zreferences.html#tall2000), [[LJ03]](https://lectures.quantecon.org/zreferences.html#lucas2003macroeconomic))  \n",
    "\n",
    "\n",
    "Let $ c_t $ be consumption at date $ t \\geq 0 $.\n",
    "\n",
    "Suppose that $ \\{\\log c_t \\}_{t=0}^\\infty $ is an additive functional described by\n",
    "\n",
    "$$\n",
    "\\log c_{t+1} - \\log c_t = \\nu + D \\cdot x_t + F \\cdot z_{t+1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "x_{t+1} = A x_t + B z_{t+1}\n",
    "$$\n",
    "\n",
    "Here $ \\{z_{t+1}\\}_{t=0}^\\infty $ is an IID sequence of $ {\\cal N}(0,I) $ random vectors.\n",
    "\n",
    "A representative household ranks consumption processes $ \\{c_t\\}_{t=0}^\\infty $ with a utility functional $ \\{V_t\\}_{t=0}^\\infty $ that satisfies\n",
    "\n",
    "\n",
    "<a id='equation-old1mf'></a>\n",
    "$$\n",
    "\\log V_t - \\log c_t = U \\cdot x_t + {\\sf u} \\tag{1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "U = \\exp(-\\delta) \\left[ I - \\exp(-\\delta) A' \\right]^{-1} D\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "{\\sf u}\n",
    "  = {\\frac {\\exp( -\\delta)}{ 1 - \\exp(-\\delta)}} {\\nu} + \\frac{(1 - \\gamma)}{2} {\\frac {\\exp(-\\delta)}{1 - \\exp(-\\delta)}}\n",
    "\\biggl| D' \\left[ I - \\exp(-\\delta) A \\right]^{-1}B + F \\biggl|^2,\n",
    "$$\n",
    "\n",
    "Here $ \\gamma \\geq 1 $ is a risk-aversion coefficient and $ \\delta > 0 $ is a rate of time preference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Consumption as a Multiplicative Process\n",
    "\n",
    "We begin by showing that consumption is a **multiplicative functional** with representation\n",
    "\n",
    "\n",
    "<a id='equation-old2mf'></a>\n",
    "$$\n",
    "\\frac{c_t}{c_0}\n",
    "= \\exp(\\tilde{\\nu}t )\n",
    "\\left( \\frac{\\tilde{M}_t}{\\tilde{M}_0} \\right)\n",
    "\\left( \\frac{\\tilde{e}(x_0)}{\\tilde{e}(x_t)} \\right) \\tag{2}\n",
    "$$\n",
    "\n",
    "where $ \\left( \\frac{\\tilde{M}_t}{\\tilde{M}_0} \\right) $ is a likelihood ratio process and $ \\tilde M_0 = 1 $.\n",
    "\n",
    "At this point, as an exercise, we ask the reader please to verify the following formulas for $ \\tilde{\\nu} $ and $ \\tilde{e}(x_t) $ as functions of $ A, B, D, F $:\n",
    "\n",
    "$$\n",
    "\\tilde \\nu =  \\nu + \\frac{H \\cdot H}{2}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\tilde e(x) = \\exp[g(x)] = \\exp \\bigl[ D' (I - A)^{-1} x \\bigr]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating a Likelihood Ratio Process Again\n",
    "\n",
    "Next, we want a program to simulate the likelihood ratio process $ \\{ \\tilde{M}_t \\}_{t=0}^\\infty $.\n",
    "\n",
    "In particular, we want to simulate 5000 sample paths of length $ T=1000 $ for the case in which $ x $ is a scalar and $ [A, B, D, F] = [0.8, 0.001, 1.0, 0.01] $ and $ \\nu = 0.005 $.\n",
    "\n",
    "After accomplishing this, we want to display a histogram of $ \\tilde{M}_T^i $ for\n",
    "$ T=1000 $.\n",
    "\n",
    "Here is code that accomplishes these tasks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def simulate_martingale_components(amf, T=1000, I=5000):\n",
    "    # Get the multiplicative decomposition\n",
    "    ν, H, g = amf.multiplicative_decomp()\n",
    "\n",
    "    # Allocate space\n",
    "    add_mart_comp = np.empty((I, T))\n",
    "\n",
    "    # Simulate and pull out additive martingale component\n",
    "    for i in range(I):\n",
    "        foo, bar = amf.lss.simulate(T)\n",
    "\n",
    "        # Martingale component is the third component\n",
    "        add_mart_comp[i, :] = bar[2, :]\n",
    "\n",
    "    mul_mart_comp = np.exp(add_mart_comp - (np.arange(T) * H**2) / 2)\n",
    "\n",
    "    return add_mart_comp, mul_mart_comp\n",
    "\n",
    "\n",
    "# Build model\n",
    "amf_2 = AMF_LSS_VAR(0.8, 0.001, 1.0, 0.01,.005)\n",
    "\n",
    "amc, mmc = simulate_martingale_components(amf_2, 1000, 5000)\n",
    "\n",
    "amcT = amc[:, -1]\n",
    "mmcT = mmc[:, -1]\n",
    "\n",
    "print(\"The (min, mean, max) of additive Martingale component in period T is\")\n",
    "print(f\"\\t ({np.min(amcT)}, {np.mean(amcT)}, {np.max(amcT)})\")\n",
    "\n",
    "print(\"The (min, mean, max) of multiplicative Martingale component in period T is\")\n",
    "print(f\"\\t ({np.min(mmcT)}, {np.mean(mmcT)}, {np.max(mmcT)})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "- The preceding min, mean, and max of the cross-section of the date\n",
    "  $ T $ realizations of the multiplicative martingale component of\n",
    "  $ c_t $ indicate that the sample mean is close to its population\n",
    "  mean of 1.  \n",
    "  \n",
    "  - This outcome prevails for all values of the horizon $ T $.  \n",
    "  \n",
    "- The cross-section distribution of the multiplicative martingale\n",
    "  component of $ c $ at date $ T $ approximates a log-normal\n",
    "  distribution well.  \n",
    "- The histogram of the additive martingale component of\n",
    "  $ \\log c_t $ at date $ T $ approximates a normal distribution\n",
    "  well.  \n",
    "\n",
    "\n",
    "Here’s a histogram of the additive martingale component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.hist(amcT, bins=25, normed=True)\n",
    "plt.title(\"Histogram of Additive Martingale Component\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here’s a histogram of the multiplicative martingale component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.hist(mmcT, bins=25, normed=True)\n",
    "plt.title(\"Histogram of Multiplicative Martingale Component\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Representing the Likelihood Ratio Process\n",
    "\n",
    "The likelihood ratio process $ \\{\\widetilde M_t\\}_{t=0}^\\infty $ can be represented as\n",
    "\n",
    "$$\n",
    "\\widetilde M_t = \\exp \\biggl( \\sum_{j=1}^t \\biggl(H \\cdot z_j -\\frac{ H \\cdot H }{2} \\biggr) \\biggr),  \\quad \\widetilde M_0 =1 ,\n",
    "$$\n",
    "\n",
    "where $ H =  [F + B'(I-A')^{-1} D] $.\n",
    "\n",
    "It follows that $ \\log {\\widetilde M}_t \\sim {\\mathcal N} ( -\\frac{t H \\cdot H}{2}, t H \\cdot H ) $ and that consequently $ {\\widetilde M}_t $ is log-normal.\n",
    "\n",
    "Let’s plot the probability density functions for $ \\log {\\widetilde M}_t $ for\n",
    "$ t=100, 500, 1000, 10000, 100000 $.\n",
    "\n",
    "Then let’s use the plots to  investigate how these densities evolve through time.\n",
    "\n",
    "We will plot the densities of $ \\log {\\widetilde M}_t $ for different values of $ t $.\n",
    "\n",
    "Note: `scipy.stats.lognorm` expects you to pass the standard deviation\n",
    "first $ (tH \\cdot H) $ and then the exponent of the mean as a\n",
    "keyword argument `scale` (`scale=`$ \\exp(-tH \\cdot H/2) $).\n",
    "\n",
    "- See the documentation [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html#scipy.stats.lognorm).  \n",
    "\n",
    "\n",
    "This is peculiar, so make sure you are careful in working with the log-normal distribution.\n",
    "\n",
    "Here is some code that tackles these tasks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def Mtilde_t_density(amf, t, xmin=1e-8, xmax=5.0, npts=5000):\n",
    "\n",
    "    # Pull out the multiplicative decomposition\n",
    "    νtilde, H, g = amf.multiplicative_decomp()\n",
    "    H2 = H * H\n",
    "\n",
    "    # The distribution\n",
    "    mdist = lognorm(np.sqrt(t * H2), scale=np.exp(-t * H2 / 2))\n",
    "    x = np.linspace(xmin, xmax, npts)\n",
    "    pdf = mdist.pdf(x)\n",
    "\n",
    "    return x, pdf\n",
    "\n",
    "\n",
    "def logMtilde_t_density(amf, t, xmin=-15.0, xmax=15.0, npts=5000):\n",
    "\n",
    "    # Pull out the multiplicative decomposition\n",
    "    νtilde, H, g = amf.multiplicative_decomp()\n",
    "    H2 = H * H\n",
    "\n",
    "    # The distribution\n",
    "    lmdist = norm(-t * H2 / 2, np.sqrt(t * H2))\n",
    "    x = np.linspace(xmin, xmax, npts)\n",
    "    pdf = lmdist.pdf(x)\n",
    "\n",
    "    return x, pdf\n",
    "\n",
    "\n",
    "times_to_plot = [10, 100, 500, 1000, 2500, 5000]\n",
    "dens_to_plot = map(lambda t: Mtilde_t_density(amf_2, t, xmin=1e-8, xmax=6.0), times_to_plot)\n",
    "ldens_to_plot = map(lambda t: logMtilde_t_density(amf_2, t, xmin=-10.0, xmax=10.0), times_to_plot)\n",
    "\n",
    "fig, ax = plt.subplots(3, 2, figsize=(8, 14))\n",
    "ax = ax.flatten()\n",
    "\n",
    "fig.suptitle(r\"Densities of $\\tilde{M}_t$\", fontsize=18, y=1.02)\n",
    "for (it, dens_t) in enumerate(dens_to_plot):\n",
    "    x, pdf = dens_t\n",
    "    ax[it].set_title(f\"Density for time {times_to_plot[it]}\")\n",
    "    ax[it].fill_between(x, np.zeros_like(pdf), pdf)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These probability density functions illustrate a **peculiar property** of log-likelihood ratio processes:\n",
    "\n",
    "- With respect to the true model probabilities, they have mathematical expectations equal to $ 1 $ for all $ t \\geq 0 $.  \n",
    "- They almost surely converge to zero.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Welfare Benefits of Reduced Random Aggregate Fluctuations\n",
    "\n",
    "Suppose in the tradition of a strand of macroeconomics (for example Tallarini [[Tal00]](https://lectures.quantecon.org/zreferences.html#tall2000), [[LJ03]](https://lectures.quantecon.org/zreferences.html#lucas2003macroeconomic)) we want to estimate the welfare benefits from removing random fluctuations around trend growth.\n",
    "\n",
    "We shall  compute how much initial consumption $ c_0 $ a representative consumer who ranks consumption streams according to [(1)](#equation-old1mf) would be willing to sacrifice to enjoy the consumption stream\n",
    "\n",
    "$$\n",
    "\\frac{c_t}{c_0} = \\exp (\\tilde{\\nu} t)\n",
    "$$\n",
    "\n",
    "rather than the stream described by equation [(2)](#equation-old2mf).\n",
    "\n",
    "We want to compute the implied percentage reduction in $ c_0 $ that the representative consumer would accept.\n",
    "\n",
    "To accomplish this, we write a function that computes the coefficients $ U $\n",
    "and $ u $ for the original values of $ A, B, D, F, \\nu $, but\n",
    "also for the case that  $ A, B, D, F = [0, 0, 0, 0] $ and\n",
    "$ \\nu = \\tilde{\\nu} $.\n",
    "\n",
    "Here’s our code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def Uu(amf, δ, γ):\n",
    "    A, B, D, F, ν = amf.A, amf.B, amf.D, amf.F, amf.ν\n",
    "    ν_tilde, H, g = amf.multiplicative_decomp()\n",
    "\n",
    "    resolv = 1 / (1 - np.exp(-δ) * A)\n",
    "    vect = F + D * resolv * B\n",
    "\n",
    "    U_risky = np.exp(-δ) * resolv * D\n",
    "    u_risky = (np.exp(-δ) / (1 - np.exp(-δ))) * (ν + (.5) * (1 - γ) * (vect**2))\n",
    "\n",
    "    U_det = 0\n",
    "    u_det = (np.exp(-δ) / (1 - np.exp(-δ))) * ν_tilde\n",
    "\n",
    "    return U_risky, u_risky, U_det, u_det\n",
    "\n",
    "# Set remaining parameters\n",
    "δ = 0.02\n",
    "γ = 2.0\n",
    "\n",
    "# Get coefficients\n",
    "U_r, u_r, U_d, u_d = Uu(amf_2, δ, γ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values of the two processes are\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\log V^r_0 &= \\log c^r_0 + U^r x_0 + u^r\n",
    "     \\\\\n",
    "    \\log V^d_0 &= \\log c^d_0 + U^d x_0 + u^d\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We look for the ratio $ \\frac{c^r_0-c^d_0}{c^r_0} $ that makes\n",
    "$ \\log V^r_0 - \\log V^d_0 = 0 $\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\underbrace{ \\log V^r_0 - \\log V^d_0}_{=0} + \\log c^d_0 - \\log c^r_0\n",
    "      &= (U^r-U^d) x_0 + u^r - u^d\n",
    "    \\\\\n",
    " \\frac{c^d_0}{ c^r_0}\n",
    "     &= \\exp\\left((U^r-U^d) x_0 + u^r - u^d\\right)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hence, the implied percentage reduction in $ c_0 $ that the\n",
    "representative consumer would accept is given by\n",
    "\n",
    "$$\n",
    "\\frac{c^r_0-c^d_0}{c^r_0} = 1 - \\exp\\left((U^r-U^d) x_0 + u^r - u^d\\right)\n",
    "$$\n",
    "\n",
    "Let’s compute this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x0 = 0.0  # initial conditions\n",
    "logVC_r = U_r * x0 + u_r\n",
    "logVC_d = U_d * x0 + u_d\n",
    "\n",
    "perc_reduct = 100 * (1 - np.exp(logVC_r - logVC_d))\n",
    "perc_reduct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find that the consumer would be willing to take a percentage reduction of initial consumption equal to around 1.081."
   ]
  }
 ],
 "metadata": {
  "filename": "multiplicative_functionals.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Multiplicative Functionals"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}