{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov switching autoregression models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook provides an example of the use of Markov switching models in Statsmodels to replicate a number of results presented in Kim and Nelson (1999). It applies the Hamilton (1989) filter the Kim (1994) smoother.\n",
    "\n",
    "This is tested against the Markov-switching models from E-views 8, which can be found at http://www.eviews.com/EViews8/ev8ecswitch_n.html#MarkovAR or the Markov-switching models of Stata 14 which can be found at http://www.stata.com/manuals14/tsmswitch.pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "import requests\n",
    "from io import BytesIO\n",
    "\n",
    "# NBER recessions\n",
    "from pandas_datareader.data import DataReader\n",
    "from datetime import datetime\n",
    "usrec = DataReader('USREC', 'fred', start=datetime(1947, 1, 1), end=datetime(2013, 4, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamilton (1989) switching model of GNP\n",
    "\n",
    "This replicates Hamilton's (1989) seminal paper introducing Markov-switching models. The model is an autoregressive model of order 4 in which the mean of the process switches between two regimes. It can be written:\n",
    "\n",
    "$$\n",
    "y_t = \\mu_{S_t} + \\phi_1 (y_{t-1} - \\mu_{S_{t-1}}) + \\phi_2 (y_{t-2} - \\mu_{S_{t-2}}) + \\phi_3 (y_{t-3} - \\mu_{S_{t-3}}) + \\phi_4 (y_{t-4} - \\mu_{S_{t-4}}) + \\varepsilon_t\n",
    "$$\n",
    "\n",
    "Each period, the regime transitions according to the following matrix of transition probabilities:\n",
    "\n",
    "$$ P(S_t = s_t | S_{t-1} = s_{t-1}) =\n",
    "\\begin{bmatrix}\n",
    "p_{00} & p_{10} \\\\\n",
    "p_{01} & p_{11}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $p_{ij}$ is the probability of transitioning *from* regime $i$, *to* regime $j$.\n",
    "\n",
    "The model class is `MarkovAutoregression` in the time-series part of `Statsmodels`. In order to create the model, we must specify the number of regimes with `k_regimes=2`, and the order of the autoregression with `order=4`. The default model also includes switching autoregressive coefficients, so here we also need to specify `switching_ar=False` to avoid that.\n",
    "\n",
    "After creation, the model is `fit` via maximum likelihood estimation. Under the hood, good starting parameters are found using a number of steps of the expectation maximization (EM) algorithm, and a quasi-Newton (BFGS) algorithm is applied to quickly find the maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the RGNP data to replicate Hamilton\n",
    "dta = pd.read_stata('https://www.stata-press.com/data/r14/rgnp.dta').iloc[1:]\n",
    "dta.index = pd.DatetimeIndex(dta.date, freq='QS')\n",
    "dta_hamilton = dta.rgnp\n",
    "\n",
    "# Plot the data\n",
    "dta_hamilton.plot(title='Growth rate of Real GNP', figsize=(12,3))\n",
    "\n",
    "# Fit the model\n",
    "mod_hamilton = sm.tsa.MarkovAutoregression(dta_hamilton, k_regimes=2, order=4, switching_ar=False)\n",
    "res_hamilton = mod_hamilton.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_hamilton.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the filtered and smoothed probabilities of a recession. Filtered refers to an estimate of the probability at time $t$ based on data up to and including time $t$ (but excluding time $t+1, ..., T$). Smoothed refers to an estimate of the probability at time $t$ using all the data in the sample.\n",
    "\n",
    "For reference, the shaded periods represent the NBER recessions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, figsize=(7,7))\n",
    "ax = axes[0]\n",
    "ax.plot(res_hamilton.filtered_marginal_probabilities[0])\n",
    "ax.fill_between(usrec.index, 0, 1, where=usrec['USREC'].values, color='k', alpha=0.1)\n",
    "ax.set_xlim(dta_hamilton.index[4], dta_hamilton.index[-1])\n",
    "ax.set(title='Filtered probability of recession')\n",
    "\n",
    "ax = axes[1]\n",
    "ax.plot(res_hamilton.smoothed_marginal_probabilities[0])\n",
    "ax.fill_between(usrec.index, 0, 1, where=usrec['USREC'].values, color='k', alpha=0.1)\n",
    "ax.set_xlim(dta_hamilton.index[4], dta_hamilton.index[-1])\n",
    "ax.set(title='Smoothed probability of recession')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the estimated transition matrix we can calculate the expected duration of a recession versus an expansion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_hamilton.expected_durations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, it is expected that a recession will last about one year (4 quarters) and an expansion about two and a half years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kim, Nelson, and Startz (1998) Three-state Variance Switching\n",
    "\n",
    "This model demonstrates estimation with regime heteroskedasticity (switching of variances) and no mean effect. The dataset can be reached at http://econ.korea.ac.kr/~cjkim/MARKOV/data/ew_excs.prn.\n",
    "\n",
    "The model in question is:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y_t & = \\varepsilon_t \\\\\n",
    "\\varepsilon_t & \\sim N(0, \\sigma_{S_t}^2)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Since there is no autoregressive component, this model can be fit using the `MarkovRegression` class. Since there is no mean effect, we specify `trend='nc'`. There are hypotheized to be three regimes for the switching variances, so we specify `k_regimes=3` and `switching_variance=True` (by default, the variance is assumed to be the same across regimes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the dataset\n",
    "ew_excs = requests.get('http://econ.korea.ac.kr/~cjkim/MARKOV/data/ew_excs.prn').content\n",
    "raw = pd.read_table(BytesIO(ew_excs), header=None, skipfooter=1, engine='python')\n",
    "raw.index = pd.date_range('1926-01-01', '1995-12-01', freq='MS')\n",
    "\n",
    "dta_kns = raw.loc[:'1986'] - raw.loc[:'1986'].mean()\n",
    "\n",
    "# Plot the dataset\n",
    "dta_kns[0].plot(title='Excess returns', figsize=(12, 3))\n",
    "\n",
    "# Fit the model\n",
    "mod_kns = sm.tsa.MarkovRegression(dta_kns, k_regimes=3, trend='nc', switching_variance=True)\n",
    "res_kns = mod_kns.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_kns.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot the probabilities of being in each of the regimes; only in a few periods is a high-variance regime probable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, figsize=(10,7))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.plot(res_kns.smoothed_marginal_probabilities[0])\n",
    "ax.set(title='Smoothed probability of a low-variance regime for stock returns')\n",
    "\n",
    "ax = axes[1]\n",
    "ax.plot(res_kns.smoothed_marginal_probabilities[1])\n",
    "ax.set(title='Smoothed probability of a medium-variance regime for stock returns')\n",
    "\n",
    "ax = axes[2]\n",
    "ax.plot(res_kns.smoothed_marginal_probabilities[2])\n",
    "ax.set(title='Smoothed probability of a high-variance regime for stock returns')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filardo (1994) Time-Varying Transition Probabilities\n",
    "\n",
    "This model demonstrates estimation with time-varying transition probabilities. The dataset can be reached at http://econ.korea.ac.kr/~cjkim/MARKOV/data/filardo.prn.\n",
    "\n",
    "In the above models we have assumed that the transition probabilities are constant across time. Here we allow the probabilities to change with the state of the economy. Otherwise, the model is the same Markov autoregression of Hamilton (1989).\n",
    "\n",
    "Each period, the regime now transitions according to the following matrix of time-varying transition probabilities:\n",
    "\n",
    "$$ P(S_t = s_t | S_{t-1} = s_{t-1}) =\n",
    "\\begin{bmatrix}\n",
    "p_{00,t} & p_{10,t} \\\\\n",
    "p_{01,t} & p_{11,t}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $p_{ij,t}$ is the probability of transitioning *from* regime $i$, *to* regime $j$ in period $t$, and is defined to be:\n",
    "\n",
    "$$\n",
    "p_{ij,t} = \\frac{\\exp\\{ x_{t-1}' \\beta_{ij} \\}}{1 + \\exp\\{ x_{t-1}' \\beta_{ij} \\}}\n",
    "$$\n",
    "\n",
    "Instead of estimating the transition probabilities as part of maximum likelihood, the regression coefficients $\\beta_{ij}$ are estimated. These coefficients relate the transition probabilities to a vector of pre-determined or exogenous regressors $x_{t-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the dataset\n",
    "filardo = requests.get('http://econ.korea.ac.kr/~cjkim/MARKOV/data/filardo.prn').content\n",
    "dta_filardo = pd.read_table(BytesIO(filardo), sep=' +', header=None, skipfooter=1, engine='python')\n",
    "dta_filardo.columns = ['month', 'ip', 'leading']\n",
    "dta_filardo.index = pd.date_range('1948-01-01', '1991-04-01', freq='MS')\n",
    "\n",
    "dta_filardo['dlip'] = np.log(dta_filardo['ip']).diff()*100\n",
    "# Deflated pre-1960 observations by ratio of std. devs.\n",
    "# See hmt_tvp.opt or Filardo (1994) p. 302\n",
    "std_ratio = dta_filardo['dlip']['1960-01-01':].std() / dta_filardo['dlip'][:'1959-12-01'].std()\n",
    "dta_filardo['dlip'][:'1959-12-01'] = dta_filardo['dlip'][:'1959-12-01'] * std_ratio\n",
    "\n",
    "dta_filardo['dlleading'] = np.log(dta_filardo['leading']).diff()*100\n",
    "dta_filardo['dmdlleading'] = dta_filardo['dlleading'] - dta_filardo['dlleading'].mean()\n",
    "\n",
    "# Plot the data\n",
    "dta_filardo['dlip'].plot(title='Standardized growth rate of industrial production', figsize=(13,3))\n",
    "plt.figure()\n",
    "dta_filardo['dmdlleading'].plot(title='Leading indicator', figsize=(13,3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time-varying transition probabilities are specified by the `exog_tvtp` parameter.\n",
    "\n",
    "Here we demonstrate another feature of model fitting - the use of a random search for MLE starting parameters. Because Markov switching models are often characterized by many local maxima of the likelihood function, performing an initial optimization step can be helpful to find the best parameters.\n",
    "\n",
    "Below, we specify that 20 random perturbations from the starting parameter vector are examined and the best one used as the actual starting parameters. Because of the random nature of the search, we seed the random number generator beforehand to allow replication of the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod_filardo = sm.tsa.MarkovAutoregression(\n",
    "    dta_filardo.iloc[2:]['dlip'], k_regimes=2, order=4, switching_ar=False,\n",
    "    exog_tvtp=sm.add_constant(dta_filardo.iloc[1:-1]['dmdlleading']))\n",
    "\n",
    "np.random.seed(12345)\n",
    "res_filardo = mod_filardo.fit(search_reps=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_filardo.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot the smoothed probability of the economy operating in a low-production state, and again include the NBER recessions for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,3))\n",
    "\n",
    "ax.plot(res_filardo.smoothed_marginal_probabilities[0])\n",
    "ax.fill_between(usrec.index, 0, 1, where=usrec['USREC'].values, color='gray', alpha=0.2)\n",
    "ax.set_xlim(dta_filardo.index[6], dta_filardo.index[-1])\n",
    "ax.set(title='Smoothed probability of a low-production state');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the time-varying transition probabilities, we can see how the expected duration of a low-production state changes over time:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_filardo.expected_durations[0].plot(\n",
    "    title='Expected duration of a low-production state', figsize=(12,3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During recessions, the expected duration of a low-production state is much higher than in an expansion."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}