ar1_processes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='ar1'></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": [
    "# AR1 Processes\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [AR1 Processes](#AR1-Processes)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [The AR(1) Model](#The-AR%281%29-Model)  \n",
    "  - [Stationarity and Asymptotic Stability](#Stationarity-and-Asymptotic-Stability)  \n",
    "  - [Ergodicity](#Ergodicity)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture we are going to study a very simple class of stochastic\n",
    "models called AR(1) processes.\n",
    "\n",
    "These simple models are used again and again in economic research to represent the dynamics of series such as\n",
    "\n",
    "- labor income  \n",
    "- dividends  \n",
    "- productivity, etc.  \n",
    "\n",
    "\n",
    "AR(1) processes can take negative values but are easily converted into positive processes when necessary by a transformation such as exponentiation.\n",
    "\n",
    "We are going to study AR(1) processes partly because they are useful and\n",
    "partly because they help us understand important concepts.\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The AR(1) Model\n",
    "\n",
    "The **AR(1) model** (autoregressive model of order 1) takes the form\n",
    "\n",
    "\n",
    "<a id='equation-can-ar1'></a>\n",
    "$$\n",
    "X_{t+1} = a X_t + b + c W_{t+1} \\tag{1}\n",
    "$$\n",
    "\n",
    "where $ a, b, c $ are scalar-valued parameters.\n",
    "\n",
    "This law of motion generates a time series $ \\{ X_t\\} $ as soon as we\n",
    "specify an initial condition $ X_0 $.\n",
    "\n",
    "This is called the **state process** and the state space is $ \\mathbb R $.\n",
    "\n",
    "To make things even simpler, we will assume that\n",
    "\n",
    "- the process $ \\{ W_t \\} $ is IID and standard normal,  \n",
    "- the initial condition $ X_0 $ is drawn from the normal distribution $ N(\\mu_0, v_0) $ and  \n",
    "- the initial condition $ X_0 $ is independent of $ \\{ W_t \\} $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Moving Average Representation\n",
    "\n",
    "Iterating backwards from time $ t $, we obtain\n",
    "\n",
    "$$\n",
    "X_t = a X_{t-1} + b +  c W_t\n",
    "        = a^2 X_{t-2} + a b + a c W_{t-1} + b + c W_t\n",
    "        = \\cdots\n",
    "$$\n",
    "\n",
    "If we work all the way back to time zero, we get\n",
    "\n",
    "\n",
    "<a id='equation-ar1-ma'></a>\n",
    "$$\n",
    "X_t = a^t X_0 + b \\sum_{j=0}^{t-1} a^j +\n",
    "        c \\sum_{j=0}^{t-1} a^j  W_{t-j} \\tag{2}\n",
    "$$\n",
    "\n",
    "Equation [(2)](#equation-ar1-ma) shows that $ X_t $ is a well defined random variable, the value of which depends on\n",
    "\n",
    "- the parameters,  \n",
    "- the initial condition $ X_0 $ and  \n",
    "- the shocks $ W_1, \\ldots W_t $ from time $ t=1 $ to the present.  \n",
    "\n",
    "\n",
    "Throughout, the symbol $ \\psi_t $ will be used to refer to the\n",
    "density of this random variable $ X_t $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution Dynamics\n",
    "\n",
    "One of the nice things about this model is that it’s so easy to trace out the sequence of distributions $ \\{ \\psi_t \\} $ corresponding to the time\n",
    "series $ \\{ X_t\\} $.\n",
    "\n",
    "To see this, we first note that $ X_t $ is normally distributed for each $ t $.\n",
    "\n",
    "This is immediate form [(2)](#equation-ar1-ma), since linear combinations of independent\n",
    "normal random variables are normal.\n",
    "\n",
    "Given that $ X_t $ is normally distributed, we will know the full distribution\n",
    "$ \\psi_t $ if we can pin down its first two moments.\n",
    "\n",
    "Let $ \\mu_t $ and $ v_t $ denote the mean and variance\n",
    "of $ X_t $ respectively.\n",
    "\n",
    "We can pin down these values from [(2)](#equation-ar1-ma) or we can use the following\n",
    "recursive expressions:\n",
    "\n",
    "\n",
    "<a id='equation-dyn-tm'></a>\n",
    "$$\n",
    "\\mu_{t+1} = a \\mu_t + b\n",
    "\\quad \\text{and} \\quad\n",
    "v_{t+1} = a^2 v_t + c^2 \\tag{3}\n",
    "$$\n",
    "\n",
    "These expressions are obtained from [(1)](#equation-can-ar1) by taking, respectively, the expectation and variance of both sides of the equality.\n",
    "\n",
    "In calculating the second expression, we are using the fact that $ X_t $\n",
    "and $ W_{t+1} $ are independent.\n",
    "\n",
    "(This follows from our assumptions and [(2)](#equation-ar1-ma).)\n",
    "\n",
    "Given the dynamics in [(2)](#equation-ar1-ma) and initial conditions $ \\mu_0,\n",
    "v_0 $, we obtain $ \\mu_t, v_t $ and hence\n",
    "\n",
    "$$\n",
    "\\psi_t = N(\\mu_t, v_t)\n",
    "$$\n",
    "\n",
    "The following code uses these facts to track the sequence of marginal\n",
    "distributions $ \\{ \\psi_t \\} $.\n",
    "\n",
    "The parameters are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "a, b, c = 0.9, 0.1, 0.5\n",
    "\n",
    "mu, v = -3.0, 0.6  # initial conditions mu_0, v_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here’s the sequence of distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "sim_length = 10\n",
    "grid = np.linspace(-5, 7, 120)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for t in range(sim_length):\n",
    "    mu = a * mu + b\n",
    "    v = a**2 * v + c**2\n",
    "    ax.plot(grid, norm.pdf(grid, loc=mu, scale=np.sqrt(v)),\n",
    "            label=f\"$\\psi_{t}$\",\n",
    "            alpha=0.7)\n",
    "\n",
    "ax.legend(bbox_to_anchor=[1.05,1],loc=2,borderaxespad=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stationarity and Asymptotic Stability\n",
    "\n",
    "Notice that, in the figure above, the sequence $ \\{ \\psi_t \\} $ seems to be converging to a limiting distribution.\n",
    "\n",
    "This is even clearer if we project forward further into the future:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_density_seq(ax, mu_0=-3.0, v_0=0.6, sim_length=60):\n",
    "    mu, v = mu_0, v_0\n",
    "    for t in range(sim_length):\n",
    "        mu = a * mu + b\n",
    "        v = a**2 * v + c**2\n",
    "        ax.plot(grid,\n",
    "                norm.pdf(grid, loc=mu, scale=np.sqrt(v)),\n",
    "                alpha=0.5)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "plot_density_seq(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, the limit does not depend on the initial condition.\n",
    "\n",
    "For example, this alternative density sequence also converges to the same limit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plot_density_seq(ax, mu_0=3.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact it’s easy to show that such convergence will occur, regardless of the initial condition, whenever $ |a| < 1 $.\n",
    "\n",
    "To see this, we just have to look at the dynamics of the first two moments, as\n",
    "given in [(3)](#equation-dyn-tm).\n",
    "\n",
    "When $ |a| < 1 $, these sequence converge to the respective limits\n",
    "\n",
    "\n",
    "<a id='equation-mu-sig-star'></a>\n",
    "$$\n",
    "\\mu^* := \\frac{b}{1-a}\n",
    "\\quad \\text{and} \\quad\n",
    "v^* = \\frac{c^2}{1 - a^2} \\tag{4}\n",
    "$$\n",
    "\n",
    "(See our [lecture on one dimensional dynamics](https://python-programming.quantecon.org/scalar_dynam.html) for background on deterministic convergence.)\n",
    "\n",
    "Hence\n",
    "\n",
    "\n",
    "<a id='equation-ar1-psi-star'></a>\n",
    "$$\n",
    "\\psi_t \\to \\psi^* = N(\\mu^*, v^*)\n",
    "\\quad \\text{as }\n",
    "t \\to \\infty \\tag{5}\n",
    "$$\n",
    "\n",
    "We can confirm this is valid for the sequence above using the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plot_density_seq(ax, mu_0=3.0)\n",
    "\n",
    "mu_star = b / (1 - a)\n",
    "std_star = np.sqrt(c**2 / (1 - a**2))  # square root of v_star\n",
    "psi_star = norm.pdf(grid, loc=mu_star, scale=std_star)\n",
    "ax.plot(grid, psi_star, 'k-', lw=2, label=\"$\\psi^*$\")\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As claimed, the sequence $ \\{ \\psi_t \\} $ converges to $ \\psi^* $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary Distributions\n",
    "\n",
    "A stationary distribution is a distribution that is a fixed\n",
    "point of the update rule for distributions.\n",
    "\n",
    "In other words, if $ \\psi_t $ is stationary, then $ \\psi_{t+j} =\n",
    "\\psi_t $ for all $ j $ in $ \\mathbb N $.\n",
    "\n",
    "A different way to put this, specialized to the current setting, is as follows: a\n",
    "density $ \\psi $ on $ \\mathbb R $ is **stationary** for the AR(1) process if\n",
    "\n",
    "$$\n",
    "X_t \\sim \\psi\n",
    "\\quad \\implies \\quad\n",
    "a X_t + b + c W_{t+1} \\sim \\psi\n",
    "$$\n",
    "\n",
    "The distribution $ \\psi^* $ in [(5)](#equation-ar1-psi-star) has this property —\n",
    "checking this is an exercise.\n",
    "\n",
    "(Of course, we are assuming that $ |a| < 1 $ so that $ \\psi^* $ is\n",
    "well defined.)\n",
    "\n",
    "In fact, it can be shown that no other distribution on $ \\mathbb R $ has this property.\n",
    "\n",
    "Thus, when $ |a| < 1 $, the AR(1) model has exactly one stationary density and that density is given by $ \\psi^* $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ergodicity\n",
    "\n",
    "The concept of ergodicity is used in different ways by different authors.\n",
    "\n",
    "One way to understand it in the present setting is that a version of the Law\n",
    "of Large Numbers is valid for $ \\{X_t\\} $, even though it is not IID.\n",
    "\n",
    "In particular, averages over time series converge to expectations under the\n",
    "stationary distribution.\n",
    "\n",
    "Indeed, it can be proved that, whenever $ |a| < 1 $, we have\n",
    "\n",
    "\n",
    "<a id='equation-ar1-ergo'></a>\n",
    "$$\n",
    "\\frac{1}{m} \\sum_{t = 1}^m h(X_t)  \\to\n",
    "\\int h(x) \\psi^*(x) dx\n",
    "    \\quad \\text{as } m \\to \\infty \\tag{6}\n",
    "$$\n",
    "\n",
    "whenever the integral on the right hand side is finite and well defined.\n",
    "\n",
    "Notes:\n",
    "\n",
    "- In [(6)](#equation-ar1-ergo), convergence holds with probability one.  \n",
    "- The textbook by [[MT09]](https://python-programming.quantecon.org/zreferences.html#meyntweedie2009) is a classic reference on ergodicity.  \n",
    "\n",
    "\n",
    "For example, if we consider the identity function $ h(x) = x $, we get\n",
    "\n",
    "$$\n",
    "\\frac{1}{m} \\sum_{t = 1}^m X_t  \\to\n",
    "\\int x \\psi^*(x) dx\n",
    "    \\quad \\text{as } m \\to \\infty\n",
    "$$\n",
    "\n",
    "In other words, the time series sample mean converges to the mean of the\n",
    "stationary distribution.\n",
    "\n",
    "As will become clear over the next few lectures, ergodicity is a very\n",
    "important concept for statistics and simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Let $ k $ be a natural number.\n",
    "\n",
    "The $ k $-th central moment of a  random variable is defined as\n",
    "\n",
    "$$\n",
    "M_k := \\mathbb E [ (X - \\mathbb E X )^k ]\n",
    "$$\n",
    "\n",
    "When that random variable is $ N(\\mu, \\sigma^2) $, it is known that\n",
    "\n",
    "$$\n",
    "M_k =\n",
    "\\begin{cases}\n",
    "    0 & \\text{ if } k \\text{ is odd} \\\\\n",
    "    \\sigma^k (k-1)!! & \\text{ if } k \\text{ is even}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Here $ n!! $ is the double factorial.\n",
    "\n",
    "According to [(6)](#equation-ar1-ergo), we should have, for any $ k \\in \\mathbb N $,\n",
    "\n",
    "$$\n",
    "\\frac{1}{m} \\sum_{t = 1}^m\n",
    "    (X_t - \\mu^* )^k\n",
    "    \\approx M_k\n",
    "$$\n",
    "\n",
    "when $ m $ is large.\n",
    "\n",
    "Confirm this by simulation at a range of $ k $ using the default parameters from the lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Write your own version of a one dimensional [kernel density\n",
    "estimator](https://en.wikipedia.org/wiki/Kernel_density_estimation),\n",
    "which estimates a density from a sample.\n",
    "\n",
    "Write it as a class that takes the data $ X $ and bandwidth\n",
    "$ h $ when initialized and provides a method $ f $ such that\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{hn} \\sum_{i=1}^n\n",
    "K \\left( \\frac{x-X_i}{h} \\right)\n",
    "$$\n",
    "\n",
    "For $ K $ use the Gaussian kernel ($ K $ is the standard normal\n",
    "density).\n",
    "\n",
    "Write the class so that the bandwidth defaults to Silverman’s rule (see\n",
    "the “rule of thumb” discussion on [this\n",
    "page](https://en.wikipedia.org/wiki/Kernel_density_estimation)). Test\n",
    "the class you have written by going through the steps\n",
    "\n",
    "1. simulate data $ X_1, \\ldots, X_n $ from distribution $ \\phi $  \n",
    "1. plot the kernel density estimate over a suitable range  \n",
    "1. plot the density of $ \\phi $ on the same figure  \n",
    "\n",
    "\n",
    "for distributions $ \\phi $ of the following types\n",
    "\n",
    "- [beta\n",
    "  distribution](https://en.wikipedia.org/wiki/Beta_distribution)\n",
    "  with $ \\alpha = \\beta = 2 $  \n",
    "- [beta\n",
    "  distribution](https://en.wikipedia.org/wiki/Beta_distribution)\n",
    "  with $ \\alpha = 2 $ and $ \\beta = 5 $  \n",
    "- [beta\n",
    "  distribution](https://en.wikipedia.org/wiki/Beta_distribution)\n",
    "  with $ \\alpha = \\beta = 0.5 $  \n",
    "\n",
    "\n",
    "Use $ n=500 $.\n",
    "\n",
    "Make a comment on your results. (Do you think this is a good estimator\n",
    "of these distributions?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    "In the lecture we discussed the following fact: for the $ AR(1) $ process\n",
    "\n",
    "$$\n",
    "X_{t+1} = a X_t + b + c W_{t+1}\n",
    "$$\n",
    "\n",
    "with $ \\{ W_t \\} $ iid and standard normal,\n",
    "\n",
    "$$\n",
    "\\psi_t = N(\\mu, s^2) \\implies \\psi_{t+1}\n",
    "= N(a \\mu + b, a^2 s^2 + c^2)\n",
    "$$\n",
    "\n",
    "Confirm this, at least approximately, by simulation. Let\n",
    "\n",
    "- $ a = 0.9 $  \n",
    "- $ b = 0.0 $  \n",
    "- $ c = 0.1 $  \n",
    "- $ \\mu = -3 $  \n",
    "- $ s = 0.2 $  \n",
    "\n",
    "\n",
    "First, plot $ \\psi_t $ and $ \\psi_{t+1} $ using the true\n",
    "distributions described above.\n",
    "\n",
    "Second, plot $ \\psi_{t+1} $ on the same figure (in a different\n",
    "color) as follows:\n",
    "\n",
    "1. Generate $ n $ draws of $ X_t $ from the $ N(\\mu, s^2) $\n",
    "  distribution  \n",
    "1. Update them all using the rule\n",
    "  $ X_{t+1} = a X_t + b + c W_{t+1} $  \n",
    "1. Use the resulting sample of $ X_{t+1} $ values to produce a\n",
    "  density estimate via kernel density estimation.  \n",
    "\n",
    "\n",
    "Try this for $ n=2000 $ and confirm that the\n",
    "simulation based estimate of $ \\psi_{t+1} $ does converge to the\n",
    "theoretical distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from numba import njit\n",
    "from scipy.special import factorial2\n",
    "\n",
    "@njit\n",
    "def sample_moments_ar1(k, m=100_000, mu_0=0.0, sigma_0=1.0, seed=1234):\n",
    "    np.random.seed(seed)\n",
    "    sample_sum = 0.0\n",
    "    x = mu_0 + sigma_0 * np.random.randn()\n",
    "    for t in range(m):\n",
    "        sample_sum += (x - mu_star)**k\n",
    "        x = a * x + b + c * np.random.randn()\n",
    "    return sample_sum / m\n",
    "\n",
    "def true_moments_ar1(k):\n",
    "    if k % 2 == 0:\n",
    "        return std_star**k * factorial2(k - 1)\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "k_vals = np.arange(6) + 1\n",
    "sample_moments = np.empty_like(k_vals)\n",
    "true_moments = np.empty_like(k_vals)\n",
    "\n",
    "for k_idx, k in enumerate(k_vals):\n",
    "    sample_moments[k_idx] = sample_moments_ar1(k)\n",
    "    true_moments[k_idx] = true_moments_ar1(k)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(k_vals, true_moments, label=\"true moments\")\n",
    "ax.plot(k_vals, sample_moments, label=\"sample moments\")\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Here is one solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "K = norm.pdf\n",
    "\n",
    "class KDE:\n",
    "\n",
    "    def __init__(self, x_data, h=None):\n",
    "\n",
    "        if h is None:\n",
    "            c = x_data.std()\n",
    "            n = len(x_data)\n",
    "            h = 1.06 * c * n**(-1/5)\n",
    "        self.h = h\n",
    "        self.x_data = x_data\n",
    "\n",
    "    def f(self, x):\n",
    "        if np.isscalar(x):\n",
    "            return K((x - self.x_data) / self.h).mean() * (1/self.h)\n",
    "        else:\n",
    "            y = np.empty_like(x)\n",
    "            for i, x_val in enumerate(x):\n",
    "                y[i] = K((x_val - self.x_data) / self.h).mean() * (1/self.h)\n",
    "            return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_kde(ϕ, x_min=-0.2, x_max=1.2):\n",
    "    x_data = ϕ.rvs(n)\n",
    "    kde = KDE(x_data)\n",
    "\n",
    "    x_grid = np.linspace(-0.2, 1.2, 100)\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(x_grid, kde.f(x_grid), label=\"estimate\")\n",
    "    ax.plot(x_grid, ϕ.pdf(x_grid), label=\"true density\")\n",
    "    ax.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import beta\n",
    "\n",
    "n = 500\n",
    "parameter_pairs= (2, 2), (2, 5), (0.5, 0.5)\n",
    "for α, β in parameter_pairs:\n",
    "    plot_kde(beta(α, β))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the kernel density estimator is effective when the underlying\n",
    "distribution is smooth but less so otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    "Here is our solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "a = 0.9\n",
    "b = 0.0\n",
    "c = 0.1\n",
    "μ = -3\n",
    "s = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "μ_next = a * μ + b\n",
    "s_next = np.sqrt(a**2 * s**2 + c**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ = lambda x: K((x - μ) / s)\n",
    "ψ_next = lambda x: K((x - μ_next) / s_next)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ψ = norm(μ, s)\n",
    "ψ_next = norm(μ_next, s_next)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 2000\n",
    "x_draws = ψ.rvs(n)\n",
    "x_draws_next = a * x_draws + b + c * np.random.randn(n)\n",
    "kde = KDE(x_draws_next)\n",
    "\n",
    "x_grid = np.linspace(μ - 1, μ + 1, 100)\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x_grid, ψ.pdf(x_grid), label=\"$\\psi_t$\")\n",
    "ax.plot(x_grid, ψ_next.pdf(x_grid), label=\"$\\psi_{t+1}$\")\n",
    "ax.plot(x_grid, kde.f(x_grid), label=\"estimate of $\\psi_{t+1}$\")\n",
    "\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulated distribution approximately coincides with the theoretical\n",
    "distribution, as predicted."
   ]
  }
 ],
 "metadata": {
  "date": 1591697892.0366428,
  "filename": "ar1_processes.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "AR1 Processes"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}