Class_10_BusinessCycles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Business Cycle Facts and Real Business Cycles Model\n",
    "\n",
    "We will take a look at the data on business cycle fluctuations in the United States and then see if we can explain their behavior with the help of the basic RBC model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "from pandas_datareader.data import DataReader"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Below we are downloading data from Federal Reserve Bank of St. Louis Economic Database (FRED) https://fred.stlouisfed.org/ and organizing them into three datasets.\n",
    "\n",
    "The first dataset, `rec`, contains monthly data on whether in a given month the US economy was in a recession state (`1`) or not (`0`) - `USREC`\n",
    "\n",
    "The second dataset, `pop`, containts quarterly data on:\n",
    "- total population in the US - `B230RC0Q173SBEA`\n",
    "- civilian noninstitutional population of 16 years of age and older - `CNP16OV`\n",
    "\n",
    "The third dataset, `fred`, contains quarterly data on:\n",
    "\n",
    "- nominal Gross Domestic Product in billions of current dollars - `GDP`\n",
    "- nominal Personal Consumption Expenditures on Nondurable Goods - `PCND`\n",
    "- nominal Personal Consumption Expenditures on Services - `PCESV`\n",
    "- nominal Personal Consumption Expenditures on Durable Goods - `PCDG`\n",
    "- nominal Fixed Private Investment - `FPI`\n",
    "- Implicit Price Deflator for Gross Domestic Product - `GDPDEF`\n",
    "- Nonfarm Business Sector: Current Dollar Output - `PRS85006053`\n",
    "- total hours worked in the nonfarm business sector (index) - `HOANBS`\n",
    "- nominal compensation per hour in the nonfarm business sector (index) - `COMPNFB`\n",
    "- interest rate on 3-month treasury bill - `TB3MS`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = '1945-01'\n",
    "end   = '2019-12'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = DataReader('USREC', 'fred', start=start, end=end)\n",
    "rec = rec.dropna()\n",
    "\n",
    "pop = DataReader(['B230RC0Q173SBEA', 'CNP16OV'], 'fred', start=start, end=end)\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(['GDP', 'PCND', 'PCESV', 'PCDG', 'FPI', 'GDPDEF',\n",
    "                   'PRS85006053', 'HOANBS', 'COMPNFB', 'TB3MS'],\n",
    "                  'fred', start=start, end=end)\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": "markdown",
   "metadata": {},
   "source": [
    "Since the population data are based on the population census, they exhibit periodic \"jumps\" that we will filter out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_hp_cycle = pd.DataFrame()\n",
    "pop_hp_trend = pd.DataFrame()\n",
    "\n",
    "for col in pop.columns:\n",
    "    pop_hp_cycle[col], pop_hp_trend[col] = sm.tsa.filters.hpfilter(np.log(pop[col]), lamb=10000)\n",
    "    \n",
    "pop_smooth = np.exp(pop_hp_trend)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "(100*pop['CNP16OV'].pct_change()).to_period('D').plot(ax=ax, lw=2)\n",
    "(100*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": "markdown",
   "metadata": {},
   "source": [
    "Below we are constructing the series that will correspond to model objects. All variables will be expressed in logarithms:\n",
    "- `Output` is real GDP per person (aged 16 and over)\n",
    "- `Consumption` is the real sum of expenditures on nondurable goods and services per person\n",
    "- `Investment` is the real sum of expenditures on durable goods and fixed private investment per person\n",
    "- `Capital` will contain real capital stock per person (to be constructed later)\n",
    "- `Hours` is hours worked per person, corrected for the difference between the entire economy and the nonfarm business sector\n",
    "- `Wages` is real compensation per hour worked\n",
    "- `Interest Rate` is the real quarterly return on the 3 month treasury bill\n",
    "- `TFP` will contain Total Factor Productivity series (to be constructed later)\n",
    "- `Productivity` is the real output per hour worked in the nonfarm business sector\n",
    "- `Price Level` is the GDP implicit price deflator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta = pd.DataFrame()\n",
    "\n",
    "dta['Output'] = np.log(fred['GDP']*10**9/fred['GDPDEF']*100\n",
    "                       /(pop_smooth['CNP16OV']*10**3))\n",
    "\n",
    "dta['Consumption'] = np.log((fred['PCND']+fred['PCESV'])*10**9/fred['GDPDEF']*100\n",
    "                            /(pop_smooth['CNP16OV']*10**3))\n",
    "\n",
    "dta['Investment'] = np.log((fred['PCDG']+fred['FPI'])*10**9/fred['GDPDEF']*100\n",
    "                           /(pop_smooth['CNP16OV']*10**3))\n",
    "\n",
    "dta['Capital'] = 0*dta['Output']\n",
    "\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",
    "\n",
    "dta['Wages'] = np.log(fred['COMPNFB']/fred['GDPDEF']*100)\n",
    "\n",
    "dta['Interest Rate'] = np.log( (1+fred['TB3MS']/100)**(1/4)/(1+fred['GDPDEF'].pct_change()) )\n",
    "\n",
    "dta['TFP'] = 0*dta['Output']\n",
    "\n",
    "dta['Productivity'] = np.log(fred['PRS85006053']/fred['GDPDEF']*100)-dta['Hours']\n",
    "\n",
    "dta['Price Level'] = np.log(fred['GDPDEF'])\n",
    "\n",
    "dta = dta.dropna()\n",
    "\n",
    "dta.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we will estimate the capital stock based on the Perpetual Inventory Method (PIM) that simply applies the capital accumulation equation:\n",
    "\n",
    "\\begin{align}\n",
    "K_{t+1} = I_{t} + \\left( 1-\\delta \\right) K_{t}\n",
    "\\end{align}\n",
    "\n",
    "We need to \"start\" this procedure from a certain level of initial capital stock, $K_{0}$, which is below calculated basing on the investment trend in the beginning of the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate Capital series using PIM\n",
    "temp = pd.DataFrame()\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",
    "\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.title('Capital stock')\n",
    "plt.show()\n",
    "\n",
    "(temp['Cap']/(fred['GDP']/fred['GDPDEF']*100)).to_period('D').plot(lw=2)\n",
    "plt.title('$K / Y$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We recover the constructed time series for capital and construct the TFP measure from the production function:\n",
    "\n",
    "\\begin{align}\n",
    "Y_{t} &= K_{t}^{\\alpha} \\left( A_{t} L_{t} \\right)^{1-\\alpha} \\\\\n",
    "Y_{t} &= TFP_{t} \\cdot K_{t}^{\\alpha} L_{t}^{1-\\alpha} \\\\\n",
    "TFP_{t} &= Y_{t} / \\left( K_{t}^{\\alpha} L_{t}^{1-\\alpha} \\right) \\\\\n",
    "\\log TFP_{t} &= \\log Y_{t} - \\alpha \\log K_{t} - \\left(1-\\alpha \\right) \\log L_{t}\n",
    "\\end{align}"
   ]
  },
  {
   "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.title('Total Factor Productivity (log)')\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, 70000)\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, 70000]\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), lamb=1600)\n",
    "gdp_cf_cycle, gdp_cf_trend = sm.tsa.filters.cffilter(np.log(RGDP_pc), low=6, high=32)"
   ]
  },
  {
   "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, 70000]\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, 70000]\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 = pd.DataFrame()\n",
    "hp_trend = pd.DataFrame()\n",
    "\n",
    "cf_cycles = pd.DataFrame()\n",
    "cf_trend = pd.DataFrame()\n",
    "\n",
    "for col in dta.columns:\n",
    "    hp_cycles[col], hp_trend[col] = sm.tsa.filters.hpfilter((100*dta[col]).dropna(), lamb=1600)\n",
    "    cf_cycles[col], cf_trend[col] = sm.tsa.filters.cffilter((100*dta[col]).dropna(), low=6, high=32)"
   ]
  },
  {
   "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, axs = plt.subplots(3, 3, figsize=(20, 15), sharex=False, sharey=False)\n",
    "\n",
    "cf_cycles[['Output','Consumption']].to_period('D').plot(ax=axs[0, 0], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Investment / 4']].to_period('D').plot(ax=axs[0, 1], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Capital']].to_period('D').plot(ax=axs[0, 2], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Hours']].to_period('D').plot(ax=axs[1, 0], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Wages']].to_period('D').plot(ax=axs[1, 1], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Interest Rate']].to_period('D').plot(ax=axs[1, 2], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','TFP']].to_period('D').plot(ax=axs[2, 0], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Productivity']].to_period('D').plot(ax=axs[2, 1], style=['k','r'], lw=2)\n",
    "\n",
    "cf_cycles[['Output','Price Level']].to_period('D').plot(ax=axs[2, 2], style=['k','r'], lw=2)\n",
    "\n",
    "for ax_a in axs:\n",
    "    for ax in ax_a:\n",
    "        ax.set_xlabel('')\n",
    "\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, lw=2)\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": [
    "rho_z    = resid_TFP.params[0]\n",
    "stderr_z = resid_TFP.resid.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Dynare import *"
   ]
  },
  {
   "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.7517, 'rho':rho_z}\n",
    "# param_values = {'alpha':0.33, sy.symbols('beta'):0.99, 'delta':0.025, 'phi':1.7517, 'rho':0.9614040418996899}\n",
    "\n",
    "model = ('-1/c + betta/c(+1)*(1+r(+1))',\n",
    "         '-phi/(1-h) + w/c',\n",
    "         '-y + z*k(-1)^alpha*h^(1-alpha)',\n",
    "         '-r + alpha*y/k(-1) - delta',\n",
    "         '-w + (1-alpha)*y/h',\n",
    "         '-k + i + (1-delta)*k(-1)',\n",
    "         '-y + c + i',\n",
    "         '-log(z) + rho*log(z(-1)) + e',\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=stderr_z)\n",
    "# rbc.SimulatedMoments(hp_filter=1600, shocks_stderr=0.0077758029493171866)"
   ]
  },
  {
   "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[:10])\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": "markdown",
   "metadata": {},
   "source": [
    "## Homework: role of Frisch elasticity\n",
    "\n",
    "A key parameter that affects how volatile are hours worked over the business cycle is the so called [Frisch elasticity of labor supply](https://en.wikipedia.org/wiki/Frisch_elasticity_of_labor_supply).\n",
    "\n",
    "It turns out, that our log-log utility function that we have used has a certain Frisch elasticity \"baked in\". For example, if the average hours worked are 1/3 of total time available, then it necessarily implies that the Frisch elasticity is equal to 2.\n",
    "\n",
    "We might want to be able to change the value of this parameter. Therefore, we will use a slightly different utility function. The household will then solve the following problem:\n",
    "\n",
    "\\begin{align*}\n",
    "\\max\\quad & U_{t}=\\mathrm{E}_{t}\\left[\\sum_{i=0}^{\\infty}\\beta^{i}\\left(\\log c_{t+i}-\\phi\\frac{h_{t+i}^{1+\\eta}}{1+\\eta}\\right)\\right]\\\\\n",
    "\\text{subject to}\\quad & c_{t}+a_{t+1}=w_{t}h_{t}+\\left(1+r_{t}\\right)a_{t}+d_{t}\n",
    "\\end{align*}\n",
    "\n",
    "where the parameter $\\eta$ (`eta`) is the inverse of the Frisch elasticity of labor supply.\n",
    "\n",
    "**Task 1**: Derive the first order condition of the households and the optimality condition for the consumption-labor choice. Modify accordingly the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_z    = 0.9614040418996899\n",
    "stderr_z = 0.0077758029493171866\n",
    "\n",
    "var = 'y c i k h w r z yh R'\n",
    "varexo = 'e'\n",
    "\n",
    "param_values_Frisch = {'alpha':0.33, sy.symbols('beta'):0.99, 'delta':0.025, 'phi':4.5511, 'rho':rho_z, 'eta':0.5}\n",
    "\n",
    "model_Frisch = ('-1/c + betta/c(+1)*(1+r(+1))',\n",
    "#          'Here goes the new consumption-labor choice equation',\n",
    "         '-y + z*k(-1)^alpha*h^(1-alpha)',\n",
    "         '-r + alpha*y/k(-1) - delta',\n",
    "         '-w + (1-alpha)*y/h',\n",
    "         '-k + i + (1-delta)*k(-1)',\n",
    "         '-y + c + i',\n",
    "         '-log(z) + rho*log(z(-1)) + e',\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)\n",
    "\n",
    "rbc_Frisch = Dynare(var, varexo, param_values_Frisch, model_Frisch, initval)\n",
    "\n",
    "rbc_Frisch.system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2**: Verify that the above model produces exactly the same steady state and very similar simulated moments as the basic RBC model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your codes go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 3**: Use the alternative parametrization below. Under the new parametrization, compute the simulated moments of the model and compare the relative standard deviation of hours to output with the results we got for the basic RBC model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "param_values_Frisch_alt = {'alpha':0.33, sy.symbols('beta'):0.99, 'delta':0.025, 'phi':2.7759, 'rho':rho_z, 'eta':0.05}\n",
    "\n",
    "# Your codes go here"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}