Example - MAR Hamilton.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:972bef9696f0dad281db904aac4c6759cda5918b1f43d0795b507dac9c927a9c"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import pandas as pd\n",
      "from mar_model import MAR"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Model Setup\n",
      "order = 4\n",
      "nstates = 2\n",
      "\n",
      "switch_ar = False\n",
      "switch_sd = False\n",
      "switch_mean = True"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the dataset from Kim and Nelson's website\n",
      "import requests\n",
      "f = requests.get('http://econ.korea.ac.kr/~cjkim/MARKOV/data/gdp4795.prn')\n",
      "\n",
      "# Hamilton's 1989 GNP dataset: Quarterly, 1947.1 - 1995.3\n",
      "data = pd.DataFrame(\n",
      "    [float(line) for line in f.content.split('\\n')[:-3]],\n",
      "    index=pd.date_range('1947-01-01', '1995-07-01', freq='QS'),\n",
      "    columns=['gnp']\n",
      ")\n",
      "data['dlgnp'] = np.log(data['gnp']).diff()*100\n",
      "data = data['1952-01-01':'1984-10-01']\n",
      "\n",
      "# NBER recessions\n",
      "from pandas.io.data import DataReader\n",
      "from datetime import datetime\n",
      "usrec = DataReader('USREC', 'fred', start=datetime(1952, 1, 1), end=datetime(1984, 12, 1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mod = MAR(data.dlgnp, order, nstates)\n",
      "\n",
      "# Set the initial probability values\n",
      "# see the file's docstring for more information\n",
      "# In this case, we have the following transition matrix:\n",
      "# | p_11 p_21 |\n",
      "# | p_12 p_22 |\n",
      "# where p_ij is the probability of transitioning from\n",
      "# state i to state j\n",
      "# Then the initial transition probabilities should be\n",
      "# *the first row* of that matrix, i.e. it should be\n",
      "# | p_11 p_21 |\n",
      "init_trans_prob = np.r_[0.75, 0.1]\n",
      "# so I am saying that p_11 = 0.75 and p_22 = 0.9\n",
      "# i.e. there is a 75% probability of staying in state 1\n",
      "# if you are already there and a 90% probability of staying\n",
      "# in state 2 if you are already there.\n",
      "\n",
      "# Set the full initial parameters array\n",
      "initial_params = np.array(np.r_[\n",
      "   np.log(init_trans_prob / (1 - init_trans_prob)),  # transition probabilities (transformed)\n",
      "   0, 0, 0, 0, # AR parameters\n",
      "   -np.log(0.8),   # standard deviation (transformed)\n",
      "   -0.2, 1.0   # Means\n",
      "])\n",
      "\n",
      "# Estimate the parameters via MLE\n",
      "res = mod.fit(initial_params, method='bfgs');\n",
      "params = res.params\n",
      "\n",
      "# Transform the unconstrained parameters to interpretable parameters\n",
      "final_params = params.copy()\n",
      "final_params[0:2] = np.exp(final_params[0:2]) / (1 + np.exp(final_params[0:2])) # transition probabilities\n",
      "final_params[6] = np.exp(-final_params[6]) # Standard deviation"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Parameter Estimates')\n",
      "print('-------------------')\n",
      "print('Transition probabilities:')\n",
      "print(mod.transition_matrix(final_params[0:2]))\n",
      "print('AR Parameters (non-switching, here)')\n",
      "print(final_params[2:6])\n",
      "print('Standard Deviation (non-switching, here)')\n",
      "print(final_params[6])\n",
      "print('Means (switching)')\n",
      "print('- in state 1:', final_params[7])\n",
      "print('- in state 2:', final_params[8])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Filter the data\n",
      "# note that this uses the unconstained parameters\n",
      "(\n",
      "    marginal_densities, filtered_joint_probabilities,\n",
      "    filtered_joint_probabilities_t1\n",
      ") = mod.filter(params);\n",
      "\n",
      "transitions = mod.separate_params(params)[0]\n",
      "\n",
      "# Smooth the data\n",
      "filtered_marginal_probabilities = mod.marginalize_probabilities(filtered_joint_probabilities[1:])\n",
      "smoothed_marginal_probabilities = mod.smooth(filtered_joint_probabilities, filtered_joint_probabilities_t1, transitions)\n",
      "\n",
      "# Save the data\n",
      "data['filtered'] = np.r_[\n",
      "    [np.NaN]*order,\n",
      "    filtered_marginal_probabilities[:,0]\n",
      "]\n",
      "data['smoothed'] = np.r_[\n",
      "    [np.NaN]*order,\n",
      "    smoothed_marginal_probabilities[:,0]\n",
      "]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "from matplotlib import dates\n",
      "fig = plt.figure(figsize=(9,9))\n",
      "\n",
      "ax = fig.add_subplot(211)\n",
      "ax.fill_between(usrec.index, 0, usrec.USREC, color='gray', alpha=0.3)\n",
      "ax.plot(data.index, data.filtered, 'k')\n",
      "ax.set(\n",
      "    xlim=('1952-04-01', '1984-12-01'),\n",
      "    ylim=(0,1),\n",
      "    title='Filtered probability of a recession (GDP: 1952:II - 1984:IV)'\n",
      ");\n",
      "\n",
      "ax = fig.add_subplot(212)\n",
      "ax.fill_between(usrec.index, 0, usrec.USREC, color='gray', alpha=0.3)\n",
      "ax.plot(data.index, data.smoothed, 'k')\n",
      "ax.set(\n",
      "    xlim=('1952-04-01', '1984-12-01'),\n",
      "    ylim=(0,1),\n",
      "    title='Smoothed probability of a recession (GDP: 1952:II - 1984:IV)'\n",
      ");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}