ar1_turningpts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d8c01aa5",
   "metadata": {},
   "source": [
    "# Forecasting  an AR(1) Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6db45500",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install arviz pymc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf6c7336",
   "metadata": {},
   "source": [
    "This lecture describes methods for forecasting statistics that are functions of future values of a univariate autogressive process.\n",
    "\n",
    "The methods are designed to take into account two possible sources of uncertainty about these statistics:\n",
    "\n",
    "- random shocks that impinge of the transition law  \n",
    "- uncertainty about the parameter values of the AR(1) process  \n",
    "\n",
    "\n",
    "We consider two sorts of statistics:\n",
    "\n",
    "- prospective values $ y_{t+j} $ of a random process $ \\{y_t\\} $ that is governed by the AR(1) process  \n",
    "- sample path properties that are defined as non-linear functions of future values $ \\{y_{t+j}\\}_{j \\geq 1} $ at time $ t $  \n",
    "\n",
    "\n",
    "**Sample path properties** are things like “time to next turning point” or “time to next recession”.\n",
    "\n",
    "To investigate sample path properties we’ll use a simulation procedure recommended by Wecker [[Wecker, 1979](https://python.quantecon.org/zreferences.html#id8)].\n",
    "\n",
    "To acknowledge uncertainty about parameters, we’ll deploy `pymc` to construct a Bayesian joint posterior distribution for unknown parameters.\n",
    "\n",
    "Let’s start with some imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d4efd90",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import arviz as az\n",
    "import pymc as pmc\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set_style('white')\n",
    "colors = sns.color_palette()\n",
    "\n",
    "import logging\n",
    "logging.basicConfig()\n",
    "logger = logging.getLogger('pymc')\n",
    "logger.setLevel(logging.CRITICAL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77d162c5",
   "metadata": {},
   "source": [
    "## A Univariate First-Order Autoregressive Process\n",
    "\n",
    "Consider the univariate AR(1) model:\n",
    "\n",
    "\n",
    "<a id='equation-ar1-tp-eq1'></a>\n",
    "$$\n",
    "y_{t+1} = \\rho y_t + \\sigma \\epsilon_{t+1}, \\quad t \\geq 0 \\tag{47.1}\n",
    "$$\n",
    "\n",
    "where the scalars $ \\rho $ and $ \\sigma $ satisfy $ |\\rho| < 1 $ and $ \\sigma > 0 $;\n",
    "$ \\{\\epsilon_{t+1}\\} $ is a sequence of i.i.d. normal random variables with mean $ 0 $ and variance $ 1 $.\n",
    "\n",
    "The initial condition $ y_{0} $ is a known number.\n",
    "\n",
    "Equation [(47.1)](#equation-ar1-tp-eq1) implies that for $ t \\geq 0 $, the conditional density of $ y_{t+1} $ is\n",
    "\n",
    "\n",
    "<a id='equation-ar1-tp-eq2'></a>\n",
    "$$\n",
    "f(y_{t+1} | y_{t}; \\rho, \\sigma) \\sim {\\mathcal N}(\\rho y_{t}, \\sigma^2) \\ \\tag{47.2}\n",
    "$$\n",
    "\n",
    "Further, equation [(47.1)](#equation-ar1-tp-eq1) also implies that for $ t \\geq 0 $, the conditional density of $ y_{t+j} $ for $ j \\geq 1 $ is\n",
    "\n",
    "\n",
    "<a id='equation-ar1-tp-eq3'></a>\n",
    "$$\n",
    "f(y_{t+j} | y_{t}; \\rho, \\sigma) \\sim {\\mathcal N}\\left(\\rho^j y_{t}, \\sigma^2 \\frac{1 - \\rho^{2j}}{1 - \\rho^2} \\right) \\tag{47.3}\n",
    "$$\n",
    "\n",
    "The predictive distribution [(47.3)](#equation-ar1-tp-eq3) that assumes that the parameters $ \\rho, \\sigma $ are known, which we express\n",
    "by conditioning on them.\n",
    "\n",
    "We also want to compute a  predictive distribution that does not condition on $ \\rho,\\sigma $ but instead takes account of our uncertainty about them.\n",
    "\n",
    "We form this predictive distribution by integrating [(47.3)](#equation-ar1-tp-eq3) with respect to a joint posterior distribution $ \\pi_t(\\rho,\\sigma | y^t) $\n",
    "that conditions on an observed history $ y^t = \\{y_s\\}_{s=0}^t $:\n",
    "\n",
    "\n",
    "<a id='equation-ar1-tp-eq4'></a>\n",
    "$$\n",
    "f(y_{t+j} | y^t)  = \\int f(y_{t+j} | y_{t}; \\rho, \\sigma) \\pi_t(\\rho,\\sigma | y^t ) d \\rho d \\sigma \\tag{47.4}\n",
    "$$\n",
    "\n",
    "Predictive distribution [(47.3)](#equation-ar1-tp-eq3) assumes that parameters $ (\\rho,\\sigma) $ are known.\n",
    "\n",
    "Predictive distribution [(47.4)](#equation-ar1-tp-eq4) assumes that parameters $ (\\rho,\\sigma) $ are uncertain, but have known probability distribution $ \\pi_t(\\rho,\\sigma | y^t ) $.\n",
    "\n",
    "We also want to compute some  predictive distributions of “sample path statistics” that might  include, for example\n",
    "\n",
    "- the time until the next “recession”,  \n",
    "- the minimum value of $ Y $ over the next 8 periods,  \n",
    "- “severe recession”, and  \n",
    "- the time until the next turning point (positive or negative).  \n",
    "\n",
    "\n",
    "To accomplish that for situations in which we are uncertain about parameter values, we shall extend Wecker’s [[Wecker, 1979](https://python.quantecon.org/zreferences.html#id8)] approach in the following way.\n",
    "\n",
    "- first simulate an initial path of length $ T_0 $;  \n",
    "- for a given prior, draw a sample of size $ N $ from the posterior joint distribution of parameters $ \\left(\\rho,\\sigma\\right) $ after observing the initial path;  \n",
    "- for each draw $ n=0,1,...,N $, simulate a “future path” of length $ T_1 $ with parameters $ \\left(\\rho_n,\\sigma_n\\right) $ and compute our three “sample path statistics”;  \n",
    "- finally, plot the desired statistics from the $ N $ samples as an empirical distribution.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9048cfb6",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "First, we’ll simulate a  sample path from which to launch our forecasts.\n",
    "\n",
    "In addition to plotting the sample path, under the assumption that the true parameter values are known,\n",
    "we’ll plot $ .9 $ and $ .95 $ coverage intervals using conditional distribution\n",
    "[(47.3)](#equation-ar1-tp-eq3) described above.\n",
    "\n",
    "We’ll also plot a bunch of samples of sequences of future values and watch where they fall relative to the coverage interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb39946b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def AR1_simulate(rho, sigma, y0, T):\n",
    "\n",
    "    # Allocate space and draw epsilons\n",
    "    y = np.empty(T)\n",
    "    eps = np.random.normal(0, sigma, T)\n",
    "\n",
    "    # Initial condition and step forward\n",
    "    y[0] = y0\n",
    "    for t in range(1, T):\n",
    "        y[t] = rho * y[t-1] + eps[t]\n",
    "        \n",
    "    return y\n",
    "\n",
    "\n",
    "def plot_initial_path(initial_path):\n",
    "    \"\"\"\n",
    "    Plot the initial path and the preceding predictive densities\n",
    "    \"\"\"\n",
    "    # Compute .9 confidence interval]\n",
    "    y0 = initial_path[-1]\n",
    "    center = np.array([rho**j * y0 for j in range(T1)])\n",
    "    vars = np.array([sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)])\n",
    "    y_bounds1_c95, y_bounds2_c95 = center + 1.96 * np.sqrt(vars), center - 1.96 * np.sqrt(vars)\n",
    "    y_bounds1_c90, y_bounds2_c90 = center + 1.65 * np.sqrt(vars), center - 1.65 * np.sqrt(vars)\n",
    "\n",
    "    # Plot\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(12, 6))\n",
    "    ax.set_title(\"Initial Path and Predictive Densities\", fontsize=15)\n",
    "    ax.plot(np.arange(-T0 + 1, 1), initial_path)\n",
    "    ax.set_xlim([-T0, T1])\n",
    "    ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1)\n",
    "\n",
    "    # Simulate future paths\n",
    "    for i in range(10):\n",
    "        y_future = AR1_simulate(rho, sigma, y0, T1)\n",
    "        ax.plot(np.arange(T1), y_future, color='grey', alpha=.5)\n",
    "    \n",
    "    # Plot 90% CI\n",
    "    ax.fill_between(np.arange(T1), y_bounds1_c95, y_bounds2_c95, alpha=.3, label='95% CI')\n",
    "    ax.fill_between(np.arange(T1), y_bounds1_c90, y_bounds2_c90, alpha=.35, label='90% CI')\n",
    "    ax.plot(np.arange(T1), center, color='red', alpha=.7, label='expected mean')\n",
    "    ax.legend(fontsize=12)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "sigma =  1\n",
    "rho = 0.9\n",
    "T0, T1 = 100, 100\n",
    "y0 = 10\n",
    "\n",
    "# Simulate\n",
    "np.random.seed(145)\n",
    "initial_path = AR1_simulate(rho, sigma, y0, T0)\n",
    "\n",
    "# Plot\n",
    "plot_initial_path(initial_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7e7fa9b",
   "metadata": {},
   "source": [
    "As functions of forecast horizon, the coverage intervals have shapes like those described in\n",
    "[https://python.quantecon.org/perm_income_cons.html](https://python.quantecon.org/perm_income_cons.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11545ac4",
   "metadata": {},
   "source": [
    "## Predictive Distributions of Path Properties\n",
    "\n",
    "Wecker [[Wecker, 1979](https://python.quantecon.org/zreferences.html#id8)] proposed using simulation techniques to characterize  predictive distribution of some statistics that are  non-linear functions of $ y $.\n",
    "\n",
    "He called these functions “path properties” to contrast them with properties of single data points.\n",
    "\n",
    "He studied two special  prospective path properties of a given series $ \\{y_t\\} $.\n",
    "\n",
    "The first was **time until the next turning point**.\n",
    "\n",
    "- he defined a **“turning point”** to be the date of the  second of two successive declines in  $ y $.  \n",
    "\n",
    "\n",
    "To examine this statistic, let $ Z $ be an indicator process\n",
    "\n",
    "$$\n",
    "Z_t(Y(\\omega)) :=  \n",
    "\\begin{cases} \n",
    "\\ 1 & \\text{if } Y_t(\\omega)< Y_{t-1}(\\omega)< Y_{t-2}(\\omega) \\geq Y_{t-3}(\\omega) \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Then the random variable **time until the next turning point**  is defined as the following **stopping time** with respect to $ Z $:\n",
    "\n",
    "$$\n",
    "W_t(\\omega):= \\inf \\{ k\\geq 1 \\mid Z_{t+k}(\\omega) = 1\\}\n",
    "$$\n",
    "\n",
    "Wecker  [[Wecker, 1979](https://python.quantecon.org/zreferences.html#id8)] also studied **the minimum value of $ Y $ over the next 8 quarters**\n",
    "which can be defined as the random variable.\n",
    "\n",
    "$$\n",
    "M_t(\\omega) := \\min \\{ Y_{t+1}(\\omega); Y_{t+2}(\\omega); \\dots; Y_{t+8}(\\omega)\\}\n",
    "$$\n",
    "\n",
    "It is interesting to study yet another possible concept of a **turning point**.\n",
    "\n",
    "Thus, let\n",
    "\n",
    "$$\n",
    "T_t(Y(\\omega)) := \n",
    "\\begin{cases}\n",
    "\\ 1 & \\text{if } Y_{t-2}(\\omega)> Y_{t-1}(\\omega) > Y_{t}(\\omega) \\ \\text{and } \\ Y_{t}(\\omega) < Y_{t+1}(\\omega) < Y_{t+2}(\\omega) \\\\\n",
    "\\ -1 & \\text{if } Y_{t-2}(\\omega)< Y_{t-1}(\\omega) < Y_{t}(\\omega) \\ \\text{and } \\ Y_{t}(\\omega) > Y_{t+1}(\\omega) > Y_{t+2}(\\omega) \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Define a **positive turning point today or tomorrow** statistic as\n",
    "\n",
    "$$\n",
    "P_t(\\omega) := \n",
    "\\begin{cases}\n",
    "\\ 1 & \\text{if } T_t(\\omega)=1 \\ \\text{or} \\ T_{t+1}(\\omega)=1 \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "This is designed to express the event\n",
    "\n",
    "- ``after one or two decrease(s), $ Y $ will grow for two consecutive quarters’’  \n",
    "\n",
    "\n",
    "Following [[Wecker, 1979](https://python.quantecon.org/zreferences.html#id8)], we can use simulations to calculate  probabilities of $ P_t $ and $ N_t $ for each period $ t $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b41bbf35",
   "metadata": {},
   "source": [
    "## A Wecker-Like Algorithm\n",
    "\n",
    "The procedure consists of the following steps:\n",
    "\n",
    "- index a sample path by $ \\omega_i $  \n",
    "- for a given date $ t $, simulate $ I $ sample paths of length $ N $  \n",
    "\n",
    "\n",
    "$$\n",
    "Y(\\omega_i) = \\left\\{ Y_{t+1}(\\omega_i), Y_{t+2}(\\omega_i), \\dots, Y_{t+N}(\\omega_i)\\right\\}_{i=1}^I\n",
    "$$\n",
    "\n",
    "- for each path $ \\omega_i $, compute the associated value of $ W_t(\\omega_i), W_{t+1}(\\omega_i), \\dots $  \n",
    "- consider the sets $ \\{W_t(\\omega_i)\\}^{T}_{i=1}, \\ \\{W_{t+1}(\\omega_i)\\}^{T}_{i=1}, \\ \\dots, \\ \\{W_{t+N}(\\omega_i)\\}^{T}_{i=1} $ as samples from the predictive distributions $ f(W_{t+1} \\mid \\mathcal y_t, \\dots) $, $ f(W_{t+2} \\mid y_t, y_{t-1}, \\dots) $, $ \\dots $, $ f(W_{t+N} \\mid y_t, y_{t-1}, \\dots) $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad93d46c",
   "metadata": {},
   "source": [
    "## Using Simulations to Approximate a Posterior Distribution\n",
    "\n",
    "The next code cells use `pymc` to compute the time $ t $ posterior distribution of $ \\rho, \\sigma $.\n",
    "\n",
    "Note that in defining the likelihood function, we choose to condition on the initial value $ y_0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d6c85b0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def draw_from_posterior(sample):\n",
    "    \"\"\"\n",
    "    Draw a sample of size N from the posterior distribution.\n",
    "    \"\"\"\n",
    "\n",
    "    AR1_model = pmc.Model()\n",
    "\n",
    "    with AR1_model:\n",
    "        \n",
    "        # Start with priors\n",
    "        rho = pmc.Uniform('rho',lower=-1.,upper=1.)  # Assume stable rho\n",
    "        sigma = pmc.HalfNormal('sigma', sigma = np.sqrt(10))\n",
    "\n",
    "        # Expected value of y at the next period (rho * y)\n",
    "        yhat = rho * sample[:-1]\n",
    "\n",
    "        # Likelihood of the actual realization.\n",
    "        y_like = pmc.Normal('y_obs', mu=yhat, sigma=sigma, observed=sample[1:])\n",
    "\n",
    "    with AR1_model:\n",
    "        trace = pmc.sample(10000, tune=5000)\n",
    "\n",
    "    # check condition\n",
    "    with AR1_model:\n",
    "        az.plot_trace(trace, figsize=(17, 6))\n",
    "    \n",
    "    rhos = trace.posterior.rho.values.flatten()\n",
    "    sigmas = trace.posterior.sigma.values.flatten()\n",
    "\n",
    "    post_sample = {\n",
    "        'rho': rhos,\n",
    "        'sigma': sigmas\n",
    "    }\n",
    "    \n",
    "    return post_sample\n",
    "\n",
    "post_samples = draw_from_posterior(initial_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8524d904",
   "metadata": {},
   "source": [
    "The graphs on the left portray posterior marginal distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4118787a",
   "metadata": {},
   "source": [
    "## Calculating Sample Path Statistics\n",
    "\n",
    "Our next step is to prepare Python code to compute our sample path statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25f96b79",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# define statistics\n",
    "def next_recession(omega):\n",
    "    n = omega.shape[0] - 3\n",
    "    z = np.zeros(n, dtype=int)\n",
    "    \n",
    "    for i in range(n):\n",
    "        z[i] = int(omega[i] <= omega[i+1] and omega[i+1] > omega[i+2] and omega[i+2] > omega[i+3])\n",
    "\n",
    "    if np.any(z) == False:\n",
    "        return 500\n",
    "    else:\n",
    "        return np.where(z==1)[0][0] + 1\n",
    "\n",
    "def minimum_value(omega):\n",
    "    return min(omega[:8])\n",
    "\n",
    "def severe_recession(omega):\n",
    "    z = np.diff(omega)\n",
    "    n = z.shape[0]\n",
    "    \n",
    "    sr = (z < -.02).astype(int)\n",
    "    indices = np.where(sr == 1)[0]\n",
    "\n",
    "    if len(indices) == 0:\n",
    "        return T1\n",
    "    else:\n",
    "        return indices[0] + 1\n",
    "\n",
    "def next_turning_point(omega):\n",
    "    \"\"\"\n",
    "    Suppose that omega is of length 6\n",
    "    \n",
    "        y_{t-2}, y_{t-1}, y_{t}, y_{t+1}, y_{t+2}, y_{t+3}\n",
    "    \n",
    "    that is sufficient for determining the value of P/N    \n",
    "    \"\"\"\n",
    "    \n",
    "    n = np.asarray(omega).shape[0] - 4\n",
    "    T = np.zeros(n, dtype=int)\n",
    "    \n",
    "    for i in range(n):\n",
    "        if ((omega[i] > omega[i+1]) and (omega[i+1] > omega[i+2]) and \n",
    "            (omega[i+2] < omega[i+3]) and (omega[i+3] < omega[i+4])):\n",
    "            T[i] = 1\n",
    "        elif ((omega[i] < omega[i+1]) and (omega[i+1] < omega[i+2]) and \n",
    "            (omega[i+2] > omega[i+3]) and (omega[i+3] > omega[i+4])):\n",
    "            T[i] = -1\n",
    "    \n",
    "    up_turn = np.where(T == 1)[0][0] + 1 if (1 in T) == True else T1\n",
    "    down_turn = np.where(T == -1)[0][0] + 1 if (-1 in T) == True else T1\n",
    "\n",
    "    return up_turn, down_turn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "890610d6",
   "metadata": {},
   "source": [
    "## Original Wecker Method\n",
    "\n",
    "Now we  apply Wecker’s original  method by simulating future paths and compute predictive distributions, conditioning\n",
    "on the true  parameters associated with the data-generating model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ab93341",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_Wecker(initial_path, N, ax):\n",
    "    \"\"\"\n",
    "    Plot the predictive distributions from \"pure\" Wecker's method.\n",
    "    \"\"\"\n",
    "    # Store outcomes\n",
    "    next_reces = np.zeros(N)\n",
    "    severe_rec = np.zeros(N)\n",
    "    min_vals = np.zeros(N)\n",
    "    next_up_turn, next_down_turn = np.zeros(N), np.zeros(N)\n",
    "    \n",
    "    # Compute .9 confidence interval]\n",
    "    y0 = initial_path[-1]\n",
    "    center = np.array([rho**j * y0 for j in range(T1)])\n",
    "    vars = np.array([sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)])\n",
    "    y_bounds1_c95, y_bounds2_c95 = center + 1.96 * np.sqrt(vars), center - 1.96 * np.sqrt(vars)\n",
    "    y_bounds1_c90, y_bounds2_c90 = center + 1.65 * np.sqrt(vars), center - 1.65 * np.sqrt(vars)\n",
    "\n",
    "    # Plot\n",
    "    ax[0, 0].set_title(\"Initial path and predictive densities\", fontsize=15)\n",
    "    ax[0, 0].plot(np.arange(-T0 + 1, 1), initial_path)\n",
    "    ax[0, 0].set_xlim([-T0, T1])\n",
    "    ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1)\n",
    "\n",
    "    # Plot 90% CI\n",
    "    ax[0, 0].fill_between(np.arange(T1), y_bounds1_c95, y_bounds2_c95, alpha=.3)\n",
    "    ax[0, 0].fill_between(np.arange(T1), y_bounds1_c90, y_bounds2_c90, alpha=.35)\n",
    "    ax[0, 0].plot(np.arange(T1), center, color='red', alpha=.7)\n",
    "\n",
    "    # Simulate future paths\n",
    "    for n in range(N):\n",
    "        sim_path = AR1_simulate(rho, sigma, initial_path[-1], T1)\n",
    "        next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path]))\n",
    "        severe_rec[n] = severe_recession(sim_path)\n",
    "        min_vals[n] = minimum_value(sim_path)\n",
    "        next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path)\n",
    "\n",
    "        if n%(N/10) == 0:\n",
    "            ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1)\n",
    "            \n",
    "    # Return next_up_turn, next_down_turn\n",
    "    sns.histplot(next_reces, kde=True, stat='density', ax=ax[0, 1], alpha=.8, label='True parameters')\n",
    "    ax[0, 1].set_title(\"Predictive distribution of time until the next recession\", fontsize=13)\n",
    "\n",
    "    sns.histplot(severe_rec, kde=False, stat='density', ax=ax[1, 0], binwidth=0.9, alpha=.7, label='True parameters')\n",
    "    ax[1, 0].set_title(r\"Predictive distribution of stopping time of growth$<-2\\%$\", fontsize=13)\n",
    "\n",
    "    sns.histplot(min_vals, kde=True, stat='density', ax=ax[1, 1], alpha=.8, label='True parameters')\n",
    "    ax[1, 1].set_title(\"Predictive distribution of minimum value in the next 8 periods\", fontsize=13)\n",
    "\n",
    "    sns.histplot(next_up_turn, kde=True, stat='density', ax=ax[2, 0], alpha=.8, label='True parameters')\n",
    "    ax[2, 0].set_title(\"Predictive distribution of time until the next positive turn\", fontsize=13)\n",
    "\n",
    "    sns.histplot(next_down_turn, kde=True, stat='density', ax=ax[2, 1], alpha=.8, label='True parameters')\n",
    "    ax[2, 1].set_title(\"Predictive distribution of time until the next negative turn\", fontsize=13)\n",
    "\n",
    "fig, ax = plt.subplots(3, 2, figsize=(15,12))\n",
    "plot_Wecker(initial_path, 1000, ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1902848e",
   "metadata": {},
   "source": [
    "## Extended Wecker Method\n",
    "\n",
    "Now we apply we apply our  “extended” Wecker method based on  predictive densities of $ y $ defined by\n",
    "[(47.4)](#equation-ar1-tp-eq4) that acknowledge posterior uncertainty in the parameters $ \\rho, \\sigma $.\n",
    "\n",
    "To approximate  the intergration on the right side of [(47.4)](#equation-ar1-tp-eq4), we  repeatedly draw parameters from the joint posterior distribution each time we simulate a sequence of future values from model [(47.1)](#equation-ar1-tp-eq1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "976f0ddc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_extended_Wecker(post_samples, initial_path, N, ax):\n",
    "    \"\"\"\n",
    "    Plot the extended Wecker's predictive distribution\n",
    "    \"\"\"\n",
    "    # Select a sample\n",
    "    index = np.random.choice(np.arange(len(post_samples['rho'])), N + 1, replace=False)\n",
    "    rho_sample = post_samples['rho'][index]\n",
    "    sigma_sample = post_samples['sigma'][index]\n",
    "\n",
    "    # Store outcomes\n",
    "    next_reces = np.zeros(N)\n",
    "    severe_rec = np.zeros(N)\n",
    "    min_vals = np.zeros(N)\n",
    "    next_up_turn, next_down_turn = np.zeros(N), np.zeros(N)\n",
    "\n",
    "    # Plot\n",
    "    ax[0, 0].set_title(\"Initial path and future paths simulated from posterior draws\", fontsize=15)\n",
    "    ax[0, 0].plot(np.arange(-T0 + 1, 1), initial_path)\n",
    "    ax[0, 0].set_xlim([-T0, T1])\n",
    "    ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1)\n",
    "\n",
    "    # Simulate future paths\n",
    "    for n in range(N):\n",
    "        sim_path = AR1_simulate(rho_sample[n], sigma_sample[n], initial_path[-1], T1)\n",
    "        next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path]))\n",
    "        severe_rec[n] = severe_recession(sim_path)\n",
    "        min_vals[n] = minimum_value(sim_path)\n",
    "        next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path)\n",
    "\n",
    "        if n % (N / 10) == 0:\n",
    "            ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1)\n",
    "        \n",
    "    # Return next_up_turn, next_down_turn\n",
    "    sns.histplot(next_reces, kde=True, stat='density', ax=ax[0, 1], alpha=.6, color=colors[1], label='Sampling from posterior')\n",
    "    ax[0, 1].set_title(\"Predictive distribution of time until the next recession\", fontsize=13)\n",
    "\n",
    "    sns.histplot(severe_rec, kde=False, stat='density', ax=ax[1, 0], binwidth=.9, alpha=.6, color=colors[1], label='Sampling from posterior')\n",
    "    ax[1, 0].set_title(r\"Predictive distribution of stopping time of growth$<-2\\%$\", fontsize=13)\n",
    "\n",
    "    sns.histplot(min_vals, kde=True, stat='density', ax=ax[1, 1], alpha=.6, color=colors[1], label='Sampling from posterior')\n",
    "    ax[1, 1].set_title(\"Predictive distribution of minimum value in the next 8 periods\", fontsize=13)\n",
    "\n",
    "    sns.histplot(next_up_turn, kde=True, stat='density', ax=ax[2, 0], alpha=.6, color=colors[1], label='Sampling from posterior')\n",
    "    ax[2, 0].set_title(\"Predictive distribution of time until the next positive turn\", fontsize=13)\n",
    "\n",
    "    sns.histplot(next_down_turn, kde=True, stat='density', ax=ax[2, 1], alpha=.6, color=colors[1], label='Sampling from posterior')\n",
    "    ax[2, 1].set_title(\"Predictive distribution of time until the next negative turn\", fontsize=13)\n",
    "\n",
    "fig, ax = plt.subplots(3, 2, figsize=(15, 12))\n",
    "plot_extended_Wecker(post_samples, initial_path, 1000, ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b125efea",
   "metadata": {},
   "source": [
    "## Comparison\n",
    "\n",
    "Finally, we plot both the original Wecker method and the extended method with parameter values drawn from the posterior together to compare the differences that emerge from pretending to know parameter values when they are actually uncertain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "613aedd5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(3, 2, figsize=(15,12))\n",
    "plot_Wecker(initial_path, 1000, ax)\n",
    "ax[0, 0].clear()\n",
    "plot_extended_Wecker(post_samples, initial_path, 1000, ax)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "date": 1731921898.1311622,
  "filename": "ar1_turningpts.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Forecasting  an AR(1) Process"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}