{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import sympy as sy\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from numpy.linalg import eig, inv\n",
    "from pandas_datareader.data import DataReader\n",
    "\n",
    "sy.init_printing(use_latex='mathjax')\n",
    "np.set_printoptions(precision=4, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import rcParams\n",
    "\n",
    "# Restore old behavior of rounding default axis ranges\n",
    "rcParams['axes.autolimit_mode'] = 'round_numbers'\n",
    "rcParams['axes.xmargin'] = 0\n",
    "rcParams['axes.ymargin'] = 0\n",
    "\n",
    "# Adjust tick placement\n",
    "rcParams['xtick.direction'] = 'in'\n",
    "rcParams['ytick.direction'] = 'in'\n",
    "rcParams['xtick.top'] = True\n",
    "rcParams['ytick.right'] = True\n",
    "\n",
    "# Disable legend frame\n",
    "rcParams['legend.frameon'] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = DataReader('USREC', 'fred', start='1947-01', end='2020-01')\n",
    "rec = rec.dropna()\n",
    "\n",
    "pop = DataReader(['B230RC0Q173SBEA','CNP16OV'], 'fred', start='1947-01', end='2020-01')\n",
    "pop = pop.dropna()\n",
    "pop = pop.resample('QS').mean()\n",
    "pop = pop.dropna()\n",
    "\n",
    "pop.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fred = DataReader(['PRS85006053','GDP','PCND','PCESV','PCDG',\n",
    "                   'FPI','HOANBS','COMPNFB','GDPDEF','TB3MS'],\n",
    "                  'fred', start='1947-01', end='2020-01')\n",
    "fred = fred.dropna()\n",
    "fred = fred.resample('QS').mean()\n",
    "fred = fred.dropna()\n",
    "fred.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fred.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_hp_cycle, pop_hp_trend = sm.tsa.filters.hpfilter(np.log(pop), lamb=10000)\n",
    "pop_smooth = np.exp(pop_hp_trend)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "pop['CNP16OV'].to_period('D').plot(ax=ax, lw=2)\n",
    "pop_smooth['CNP16OV'].to_period('D').plot(ax=ax, lw=2)\n",
    "plt.legend(['raw', 'filtered'])\n",
    "plt.show()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "pop['CNP16OV'].pct_change().to_period('D').plot(ax=ax, lw=2)\n",
    "pop_smooth['CNP16OV'].pct_change().to_period('D').plot(ax=ax, lw=2)\n",
    "plt.legend(['raw', 'filtered'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta = fred[['GDP']]\n",
    "dta.is_copy = False\n",
    "\n",
    "dta['Output'] = np.log(fred['GDP']*10**9/fred['GDPDEF']*100\n",
    "                       /(pop_smooth['CNP16OV']*10**3))\n",
    "dta['Consumption'] = np.log((fred['PCND']+fred['PCESV'])*10**9/fred['GDPDEF']*100\n",
    "                            /(pop_smooth['CNP16OV']*10**3))\n",
    "dta['Investment'] = np.log((fred['PCDG']+fred['FPI'])*10**9/fred['GDPDEF']*100\n",
    "                           /(pop_smooth['CNP16OV']*10**3))\n",
    "dta['Capital'] = 0*dta['Output']\n",
    "dta['Hours'] = np.log(fred['HOANBS']*100*fred['GDP']/np.mean(fred['GDP']['2010-01':'2010-10'])\n",
    "                      /fred['PRS85006053']/pop_smooth['CNP16OV'])\n",
    "dta['Wages'] = np.log(fred['COMPNFB']/fred['GDPDEF']*100)\n",
    "dta['Interest Rate'] = (1+fred['TB3MS']/100)**(1/4)/(1+fred['GDPDEF'].pct_change())\n",
    "dta['TFP'] = 0*dta['Output']\n",
    "dta['Productivity'] = np.log(fred['PRS85006053']/fred['GDPDEF']*100)-dta['Hours']\n",
    "dta['Price Level'] = np.log(fred['GDPDEF'])\n",
    "\n",
    "dta = dta.drop('GDP', 1)\n",
    "dta = dta.dropna()\n",
    "\n",
    "dta.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate Capital series using PIM\n",
    "temp = fred[['GDP']]\n",
    "temp.is_copy = False\n",
    "\n",
    "temp['Inv'] = (fred['PCDG']+fred['FPI'])/fred['GDPDEF']*100/4\n",
    "temp['LnInv'] = np.log(temp['Inv'])\n",
    "temp['t'] = np.arange(len(temp['Inv']))\n",
    "\n",
    "trend = smf.ols(formula='LnInv ~ t', data=temp).fit()\n",
    "intercept, slope = trend.params\n",
    "\n",
    "delta = 0.025\n",
    "K = np.zeros(len(temp['LnInv']))\n",
    "K_init = np.exp(intercept-slope)/(slope+delta)\n",
    "K[0] = (1-delta)*K_init+temp['Inv'][0]\n",
    "for i in range(1,len(temp['Inv'])):\n",
    "    K[i] = (1-delta)*K[i-1]+temp['Inv'][i]\n",
    "temp['Cap'] = K\n",
    "\n",
    "temp['Cap'].to_period('D').plot(lw=2)\n",
    "plt.show()\n",
    "\n",
    "(temp['Cap']/(fred['GDP']/fred['GDPDEF']*100)).to_period('D').plot(lw=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta['Capital'] = np.log(temp['Cap']*10**9/(pop_smooth['CNP16OV']*10**3))\n",
    "\n",
    "α = 1/3\n",
    "dta['TFP'] = dta['Output']-α*dta['Capital']-(1-α)*dta['Hours']\n",
    "dta['TFP'].to_period('D').plot(lw=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RGDP_pc = (fred['GDP']*10**9/fred['GDPDEF']*100/(pop_smooth['B230RC0Q173SBEA']*10**3)).dropna()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "RGDP_pc.to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "\n",
    "ax.set_ylim(10000, 60000)\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP (2009 dollars)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "\n",
    "ticks = [10000, 20000, 30000, 40000, 50000, 60000]\n",
    "ax.set_yticks(np.log(ticks))\n",
    "ax.set_yticklabels(ticks)\n",
    "\n",
    "ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, log scale')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "(100*RGDP_pc.pct_change(4)).to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "\n",
    "avg = np.mean(100*RGDP_pc.to_period('D').pct_change(4))\n",
    "\n",
    "ax.set_ylim(-6, 12)\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.hlines(avg, dta.index[0], dta.index[-1], color='r', linewidth=2)\n",
    "ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, year-over-year change (%)')\n",
    "plt.show()\n",
    "\n",
    "print('Average growth rate (%) =', avg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdp_lin_cycle, gdp_lin_trend = sm.tsa.filters.hpfilter(np.log(RGDP_pc), lamb=1e9)\n",
    "gdp_hp_cycle, gdp_hp_trend = sm.tsa.filters.hpfilter(np.log(RGDP_pc))\n",
    "gdp_cf_cycle, gdp_cf_trend = sm.tsa.filters.cffilter(np.log(RGDP_pc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "gdp_lin_trend.plot(ax=ax, lw=2, style='r-')\n",
    "\n",
    "ticks = [10000, 20000, 30000, 40000, 50000, 60000]\n",
    "ax.set_yticks(np.log(ticks))\n",
    "ax.set_yticklabels(ticks)\n",
    "\n",
    "ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, exponential trend')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "(100*gdp_lin_cycle).to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, exponential trend residual (%)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "gdp_hp_trend.plot(ax=ax, lw=2, style='r-')\n",
    "\n",
    "ticks = [10000, 20000, 30000, 40000, 50000, 60000]\n",
    "ax.set_yticks(np.log(ticks))\n",
    "ax.set_yticklabels(ticks)\n",
    "\n",
    "ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, Hodrick-Prescott trend')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "(100*gdp_hp_cycle.to_period('D')).plot(ax=ax, lw=2, style='k-')\n",
    "ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)\n",
    "\n",
    "ax.set_ylim(-6, 6)\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, Hodrick-Prescott residual (%)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "(100*gdp_cf_cycle).to_period('D').plot(ax=ax, lw=2, style='k-')\n",
    "ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.title('US real per capita GDP, Christiano-Fitzgerald residual (%)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.corrcoef(gdp_hp_cycle, gdp_cf_cycle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_cycles, hp_trend = sm.tsa.filters.hpfilter((100*dta).dropna())\n",
    "cf_cycles, cf_trend = sm.tsa.filters.cffilter((100*dta).dropna())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Standard Deviations')\n",
    "print(hp_cycles.std())\n",
    "\n",
    "print('')\n",
    "print('Autocorrelations')\n",
    "a = list(dta.columns.values)\n",
    "for i in range(len(a)):\n",
    "    print(dta.columns.values[i], '  \\t\\t', hp_cycles[dta.columns.values[i]].autocorr())\n",
    "\n",
    "print('')\n",
    "print('Correlations')\n",
    "print(hp_cycles.corr(method='pearson'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cf_cycles['Investment / 4'] = cf_cycles['Investment'] / 4\n",
    "\n",
    "fig, ((ax1, ax2, ax3), (ax4, ax5, ax6), \n",
    "      (ax7, ax8, ax9)) = plt.subplots(3, 3, figsize=(20,15), sharex=False, sharey=False)\n",
    "\n",
    "cf_cycles[['Output','Consumption']].to_period('D').plot(ax=ax1, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Investment / 4']].to_period('D').plot(ax=ax2, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Capital']].to_period('D').plot(ax=ax3, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Hours']].to_period('D').plot(ax=ax4, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Wages']].to_period('D').plot(ax=ax5, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Interest Rate']].to_period('D').plot(ax=ax6, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','TFP']].to_period('D').plot(ax=ax7, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Productivity']].to_period('D').plot(ax=ax8, style=['k','r'])\n",
    "\n",
    "cf_cycles[['Output','Price Level']].to_period('D').plot(ax=ax9, style=['k','r'])\n",
    "\n",
    "plt.savefig('US_CF.pdf', transparent=True, bbox_inches='tight', pad_inches=0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "cf_cycles[['Output','Consumption','Investment / 4','Hours']].to_period('D').plot(ax=ax)\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "plt.legend(ncol=2, frameon=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta['t'] = np.arange(len(dta['TFP']))\n",
    "trend_TFP = smf.ols(formula='TFP ~ t', data=dta).fit()\n",
    "dta['TFP_resid'] = trend_TFP.resid\n",
    "intercept_TFP, slope_TFP = trend_TFP.params\n",
    "print(trend_TFP.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta['TFP'].to_period('D').plot(lw=2)\n",
    "(intercept_TFP + slope_TFP*dta['t']).plot(lw=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta['TFP_resid'].to_period('D').plot(lw=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resid_TFP = smf.ols(formula='TFP_resid ~ TFP_resid.shift() -1', data=dta).fit()\n",
    "print(resid_TFP.summary())\n",
    "\n",
    "print('')\n",
    "print('TFP residual autocorrelation    =', resid_TFP.params[0])\n",
    "print('TFP residual standard deviation =', resid_TFP.resid.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Dynare(object):\n",
    "    \n",
    "    def __init__(self, var, varexo, param_values, model, initval):\n",
    "        self.var = var\n",
    "        self.varexo = varexo\n",
    "        self.param_values = param_values\n",
    "        self.model = model\n",
    "        self.initval = initval\n",
    "        \n",
    "        self.TranslateInputs()\n",
    "        \n",
    "        return None\n",
    "    \n",
    "    def TranslateInputs(self):\n",
    "        \"\"\"\n",
    "        Convert input strings into sympy objects\n",
    "        and translate timing\n",
    "        \"\"\"\n",
    "        \n",
    "        # Endogenous variables\n",
    "        self.var_symbols = sy.symbols(self.var)\n",
    "        \n",
    "        self.n = len(self.var_symbols)\n",
    "\n",
    "        self.current = ()\n",
    "        self.future  = ()\n",
    "        self.past    = ()\n",
    "        self.steadys = ()\n",
    "\n",
    "        for i in range(self.n):\n",
    "            self.current += (sy.symbols(str(self.var_symbols[i])+str('_{t}'))),\n",
    "            self.future  += (sy.symbols(str(self.var_symbols[i])+str('_{t+1}'))),\n",
    "            self.past    += (sy.symbols(str(self.var_symbols[i])+str('_{t-1}'))),\n",
    "            self.steadys += (sy.symbols(str(self.var_symbols[i])+str('_{ss}'))),\n",
    "            \n",
    "        \n",
    "        # Exogenous variables\n",
    "        self.varexo_symbols = sy.symbols(self.varexo),\n",
    "        \n",
    "        self.p = len(self.varexo_symbols)\n",
    "        \n",
    "        self.shocks = ()\n",
    "        \n",
    "        for i in range(self.p):\n",
    "            self.shocks += (sy.symbols(str(self.varexo_symbols[i])+str('_{t}'))),\n",
    "            \n",
    "        # Model equations\n",
    "        self.symbol_dict = {'betta':sy.symbols('beta'), 'gama':sy.symbols('gamma')}\n",
    "        self.timing_dict = {}\n",
    "        \n",
    "        for i in range(self.n):\n",
    "            self.timing_dict[sy.sympify(str(self.var_symbols[i])+'(+1)')] = self.future[i]\n",
    "            self.timing_dict[sy.sympify(str(self.var_symbols[i])+'(-1)')] = self.past[i]\n",
    "            self.timing_dict[sy.sympify(str(self.var_symbols[i]))]        = self.current[i]\n",
    "        \n",
    "        for i in range(self.p):\n",
    "            self.timing_dict[sy.sympify(str(self.varexo_symbols[i]))] = self.shocks[i]\n",
    "        \n",
    "        self.model_symbols = sy.sympify(self.model)\n",
    "        self.model_symbols = self.model_symbols.subs(self.timing_dict)\n",
    "        self.model_symbols = self.model_symbols.subs(self.symbol_dict)\n",
    "        \n",
    "        try:\n",
    "            temp = len(self.model_symbols)\n",
    "        except:\n",
    "            self.model_symbols = self.model_symbols,\n",
    "        \n",
    "        self.system = sy.Matrix(self.model_symbols)\n",
    "        \n",
    "        return None\n",
    "    \n",
    "    def SteadySystem(self):\n",
    "        \"\"\"\n",
    "        Removes lead/lag structure from the model to prepare for steady state calculation\n",
    "        \"\"\"\n",
    "        \n",
    "        self.steady_vars = {}\n",
    "        for i in range(self.n):\n",
    "            self.steady_vars[self.current[i]] = self.steadys[i]\n",
    "            self.steady_vars[self.future[i]]  = self.steadys[i]\n",
    "            self.steady_vars[self.past[i]]    = self.steadys[i]\n",
    "        for i in range(self.p):\n",
    "            self.steady_vars[self.shocks[i]]  = 0\n",
    "        \n",
    "        self.steady_system = self.system.subs(self.steady_vars)\n",
    "        \n",
    "        return self.steady_system\n",
    "    \n",
    "    def SteadyValues(self):\n",
    "        \"\"\"\n",
    "        Numerically solves for the steady state of the system given initval\n",
    "        \"\"\"\n",
    "        \n",
    "        self.SteadySystem()\n",
    "        \n",
    "        try:\n",
    "            ss = sy.nsolve(self.steady_system.subs(self.param_values), self.steadys, self.initval)\n",
    "        except:\n",
    "            raise RuntimeError('Adjust initial values')\n",
    "        ss = ss.T.tolist()\n",
    "\n",
    "        self.steady_values = {}\n",
    "        for i in range(self.n):\n",
    "            self.steady_values[self.steadys[i]] = np.float(ss[0][i])\n",
    "                \n",
    "        return self.steady_values\n",
    "    \n",
    "    def steady(self):\n",
    "        \"\"\"\n",
    "        Prints out steady state values for the user\n",
    "        \"\"\"\n",
    "        \n",
    "        self.SteadyValues()\n",
    "        \n",
    "        print('\\n' + 'STEADY-STATE RESULTS' + '\\n')\n",
    "\n",
    "        for i in range(self.n):\n",
    "            print(str(self.var_symbols[i]), '\\t%.4f' % self.steady_values[self.steadys[i]])\n",
    "        \n",
    "        return None\n",
    "    \n",
    "    def resid(self):\n",
    "        \n",
    "        self.SteadyValues()\n",
    "        \n",
    "        temp = self.steady_system.subs(self.param_values).subs(self.steady_values)\n",
    "        \n",
    "        print('\\n' + 'Residuals of the static equations' + '\\n')\n",
    "        \n",
    "        for i in range(self.n):\n",
    "            print('Equation number', i, ': %.4f' % temp[i])\n",
    "        \n",
    "        return None\n",
    "        \n",
    "    \n",
    "    def TimeIteration(self):\n",
    "        \"\"\"\n",
    "        Solves the first-order approximation of the model using time iteration (thanks to Pontus Rendahl)\n",
    "        \"\"\"\n",
    "        \n",
    "        self.SteadyValues()\n",
    "        \n",
    "        self.A_symb = self.system.jacobian(self.past)\n",
    "        self.B_symb = self.system.jacobian(self.current)\n",
    "        self.C_symb = self.system.jacobian(self.future)\n",
    "        self.D_symb = self.system.jacobian(self.shocks)\n",
    "        \n",
    "        self.A = np.array(self.A_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)\n",
    "        self.B = np.array(self.B_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)\n",
    "        self.C = np.array(self.C_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)\n",
    "        self.D = np.array(self.D_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)\n",
    "        \n",
    "        self.metric = 1\n",
    "        self.F = np.zeros((self.n, self.n))\n",
    "        self.S = np.zeros((self.n, self.n))\n",
    "\n",
    "        # Add maxit to while loop?\n",
    "        \n",
    "        while self.metric > 1e-13:\n",
    "            self.F = inv(self.B + self.C @ self.F) @ (-self.A)\n",
    "            self.S = inv(self.B + self.A @ self.S) @ (-self.C)\n",
    "\n",
    "            self.metric = np.max(np.max(np.abs(self.A + self.B @ self.F + self.C @ self.F @ self.F)))\n",
    "\n",
    "        self.Q = -inv(self.B + self.C @ self.F) @ self.D\n",
    "        \n",
    "        # Need formal BK check?\n",
    "\n",
    "        if sum(eig(self.F)[0] > 1) != 0:\n",
    "            raise RuntimeError('Blanchard Kahn conditions are not satisfied: no stable equilibrium')\n",
    "        if sum(eig(self.S)[0] > 1) != 0:\n",
    "            raise RuntimeError('Blanchard Kahn conditions are not satisfied: indeterminacy')\n",
    "        \n",
    "        return None\n",
    "    \n",
    "    def SimulatedMoments(self, hp_filter=None, shocks_stderr=0.01, periods=10000):\n",
    "        \n",
    "        self.TimeIteration()\n",
    "        \n",
    "        x = np.zeros((self.n, periods))\n",
    "        ɛ = np.zeros((self.p, periods))\n",
    "        \n",
    "        for i in range(self.p):\n",
    "            ɛ[i, :] = shocks_stderr * np.random.randn(periods)\n",
    "\n",
    "        for t in range(1, periods):\n",
    "            x[:, t] = self.F @ x[:, t-1] + self.Q @ ɛ[:, t]\n",
    "        \n",
    "        print('SIMULATED MOMENTS')\n",
    "        print('')\n",
    "        print('VARIABLE \\t STD. DEV.')\n",
    "            \n",
    "        if hp_filter == None:\n",
    "            for i in range(self.n):\n",
    "                print(str(self.var_symbols[i]), np.std(x[i, :]))\n",
    "        else:\n",
    "            self.SteadyValues()\n",
    "            try:\n",
    "                for i in range(self.n):\n",
    "                    hp_cycle, hp_trend = sm.tsa.filters.hpfilter(100*np.log(x[i, :]+self.steady_values[self.steadys[i]]), \n",
    "                                                                 lamb=hp_filter)\n",
    "                    print(str(self.var_symbols[i]), '\\t\\t {:.4f}'.format(np.std(hp_cycle)))\n",
    "            except:\n",
    "                print('Error: hp_filter takes only numbers as parameters')\n",
    "            \n",
    "        return None\n",
    "    \n",
    "    def stoch_simul(self, irf=40):\n",
    "        \n",
    "        self.TimeIteration()\n",
    "\n",
    "        FT = self.F.T\n",
    "        QT = self.Q.T\n",
    "\n",
    "        print('\\n'+'POLICY AND TRANSITION FUNCTIONS'+'\\n')\n",
    "\n",
    "        header = '\\t'\n",
    "        for v in self.var_symbols:\n",
    "            header += '\\t' + str(v)\n",
    "        print(header)\n",
    "        \n",
    "        line = ''\n",
    "        for i in range(self.n):\n",
    "            line += '\\t%.4f' % self.steady_values[self.steadys[i]]\n",
    "        print('Constant' + line)\n",
    "\n",
    "        for i in range(self.n):\n",
    "            if (FT[i] != np.zeros((1, self.n))).any():\n",
    "                line = '\\t'\n",
    "                for j in range(self.n):\n",
    "                    line += '\\t%.4f' % FT[i, j]\n",
    "                print(str(self.var_symbols[i])+'(-1)', line)\n",
    "        for i in range(self.p):\n",
    "            line = '\\t'\n",
    "            for j in range(self.n):\n",
    "                line += '\\t%.4f' % QT[i, j]\n",
    "            print(str(self.varexo_symbols[i]), '   ', line)\n",
    "        \n",
    "        # Impulse response functions\n",
    "        if irf > 0:\n",
    "            x = np.zeros((self.n, irf+2))\n",
    "            ɛ = np.zeros((self.p, irf+2))\n",
    "            ɛ[:, 1] = 1\n",
    "\n",
    "            for t in range(1, irf+2):\n",
    "                x[:, t] = self.F @ x[:, t-1] + self.Q @ ɛ[:, t]\n",
    "\n",
    "            for i in range(self.n):\n",
    "                plt.plot(x[i, 1:].T, 'k')\n",
    "                plt.hlines(0, 0, irf, 'r')\n",
    "                plt.title(str(self.var_symbols[i]))\n",
    "                plt.show()\n",
    "            \n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = 'y c i k h w r z yh R'\n",
    "varexo = 'e'\n",
    "\n",
    "param_values = {'alpha':0.33, sy.symbols('beta'):0.99, 'delta':0.025, 'phi':1.75, 'sigma':1, 'rho':0.962159059423142}\n",
    "\n",
    "model = ('-log(z) + rho*log(z(-1)) + e',\n",
    "         '-k + i + (1-delta)*k(-1)',\n",
    "         '-c^(-sigma) + betta*c(+1)**(-sigma)*(1+r(+1))',\n",
    "         '-y + z*k(-1)^alpha*h^(1-alpha)',\n",
    "         '-r + alpha*y/k(-1) - delta',\n",
    "         '-w + (1-alpha)*y/h',\n",
    "         '-phi/(1-h) + w/c',\n",
    "         '-y + c + i',\n",
    "         '-yh + y/h',\n",
    "         '-R + 1+r')\n",
    "\n",
    "initval = (1, 0.8, 0.2, 10, 0.33, 2, 0.01, 1, 3, 1.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbc = Dynare(var, varexo, param_values, model, initval)\n",
    "\n",
    "rbc.system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "rbc.steady()\n",
    "\n",
    "rbc.stoch_simul(irf=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbc.SimulatedMoments(hp_filter=1600, shocks_stderr=0.007828026718977966)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Standard Deviations')\n",
    "print(hp_cycles.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
