Lab_3_RBC.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",
    "*Recession is a \"business cycle contraction when there is a general decline in economic activity. Recessions generally occur when there is a widespread drop in spending (an adverse demand shock). This may be triggered by various events, such as a financial crisis, an external trade shock, an adverse supply shock, the bursting of an economic bubble, or a large-scale anthropogenic or natural disaster (e.g. a pandemic).\"*\n",
    "\n",
    "*In the United States, it is defined as \"a significant decline in economic activity spread across the market, lasting more than a few months, normally visible in real GDP, real income, employment, industrial production, and wholesale-retail sales\".*\n",
    "\n",
    "*In the United Kingdom, it is defined as a negative economic growth for two consecutive quarters.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
    "*Civilian noninstitutional population is defined as persons 16 years of age and older residing in the 50 states and \n",
    "the District of Columbia, who are not inmates of institutions (e.g., penal and mental facilities, homes for the aged), \n",
    "and who are not on active duty in the Armed Forces.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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   = '2021-06'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = DataReader('USREC', 'fred', start=start, end=end)\n",
    "rec = rec.dropna()\n",
    "rec.tail()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop = DataReader(['B230RC0Q173SBEA', 'CNP16OV'], 'fred', start=start, end=end)\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.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. \n",
    "The time series are composed of three components: seasonal, cyclical, and trend. Assume that we have seasonally adjusted time series and we want to decompose time series into trend and cyclical components.<br> \n",
    "The Hodrick–Prescott filter (also known as Hodrick–Prescott decomposition) is a tool to remove the cyclical component of a time series from raw data.<br> \n",
    "Let $y_{t}$, for $t=1,2,...,T$ denote the logarithms of a time series variable. The series $y_{t}$, is made up of a trend component $\\tau_{t}$, and a cyclical component $c_{t}$ such that $y_{t}=\\tau_{t}+c_{t}$. Given an adequately chosen, positive value of $\\lambda$ , the filter will decompose $y_{t}$ into $\\tau_{t}$ and $c_{t}$.<br> \n",
    "$\\lambda$ is the Hodrick-Prescott smoothing parameter. A value of $1600$ is suggested for quarterly data. Ravn and Uhlig suggest using a value of $6.25$ for annual data and $129600$ for monthly data."
   ]
  },
  {
   "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=1600)\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']['2012-01':'2012-10']))\n",
    "                      /(fred['PRS85006053']*(pop_smooth['CNP16OV']/np.mean(pop_smooth['CNP16OV']['2012-01':'2012-10']))))\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",
    "print(trend.summary())\n",
    "\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 (%) =', round(avg,2))"
   ]
  },
  {
   "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, linear 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, linear 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']].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": "markdown",
   "metadata": {},
   "source": [
    "#  Real Business Cycles (RBC) models \n",
    "\n",
    "Dynamic Stochastic General Equilibrium (DSGE) models aimto replicate business cycle behavior of real-world economies\n",
    "\n",
    "– Dynamic: forward-looking behavior of agents\n",
    "\n",
    "– Stochastic: the economy is subject to shocks\n",
    "\n",
    "– GE: what happens in one market influences other markets\\smallskip\n",
    "\n",
    "We can generate quantitative predictions on short-term movements of macro variables and compare them with the data\\medskip\n",
    "\n",
    "We use those models to\n",
    "\n",
    "– Simulate counterfactual scenarios\n",
    "\n",
    "– Explain past developments (historical decomposition)\n",
    "\n",
    "– Construct forecasts (conditional and uncoditional)\n",
    "\n",
    "– Perform policy experiments\\smallskip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The RBC Model\n",
    "Ramsey model with endogenous labor supplyand stochastic “technology” shocks<br><br>\n",
    "Closed economy with no government<br><br>\n",
    "Perfect competition<br><br>\n",
    "Single final good with price normalized to 1 – all other prices are real<br><br>\n",
    "Two groups of representative agents\n",
    "\n",
    "– Households\n",
    "\n",
    "– Firms\\smallskip\n",
    "\n",
    "Rational expectations – agents make no systematic forecast errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Household\n",
    "\n",
    "A representative household solves the expected utility maximization problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max_{c_{t},a_{t+1}}\\quad & U_{t}=E_{t}\\left[\\sum_{t=0}^{\\infty}\\beta^{t}\\left(\\log c_{t}-\\phi \\log \\left(1-h_{t}\\right)\\right)\\right]\\\\\n",
    "\\text{subject to}\\quad & c_{t}+a_{t+1}=\\left(1+r_{t}\\right)a_{t}+w_{t}h_{t}+d_{t}\n",
    "\\end{align}\n",
    "\n",
    "where:\n",
    "\n",
    "$\\beta$\t- discount factor<br>\n",
    "$c$\t- per capita consumption<br>\n",
    "$\\phi$\t- relative preference for leisure<br>\n",
    "$h$\t- per capita hours (as fraction of total available time)<br>\n",
    "$a$\t- per capita assets (physical capital)<br>\n",
    "$r$\t- real interest rate<br>\n",
    "$w$\t- real wage per hour<br>\n",
    "$d$\t- per capita dividends\n",
    "\n",
    "Write down the Lagrangian:\n",
    "\\begin{align}\n",
    "\\mathcal{L}\t=\\sum_{t=0}^{\\infty}\\beta^{t}\\mathrm{E}_{t}\\left[\\ln c_{t}+\\phi\\ln\\left(1-h_{t}\\right)\\right]+\\sum_{t=0}^{\\infty}\\beta^{t}\\mathrm{E}_{t}\\left[\\lambda_{t}\\left[\\left(1+r_{t}\\right)a_{t}+w_{t}h_{t}+d_{t}-a_{t+1}-c_{t}\\right]\\right]\n",
    "\\end{align}\n",
    "\n",
    "FOCs:\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial c_{t}} &\t=  \\mathrm{E}_{t}\\left[\\frac{1}{c_{t}}\\right]-\\mathrm{E}_{t}\\left[\\lambda_{t}\\right]=0\\to\\quad & \\lambda_{t} & =\\frac{1}{c_{t}}\\\\\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial h_{t}} &\t=  \\mathrm{E}_{t}\\left[-\\frac{\\phi}{1-h_{t}}\\right]+\\mathrm{E}_{t}\\left[\\lambda_{t}w_{t}\\right]=0\\to\\quad & \\lambda_{t} & =\\frac{\\phi}{w_{t}\\left(1-h_{t}\\right)}\\\\\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial a_{t+1}} &\t=  - \\mathrm{E}_{t}\\left[\\lambda_{t}\\right]+\\beta\\mathrm{E}_{t}\\left[\\lambda_{t+1}\\left(1+r_{t+1}\\right)\\right]=0\n",
    "\t\\to\\quad & \\lambda_{t} & =\\beta\\mathrm{E}_{t}\\left[\\lambda_{t+1}\\left(1+r_{t+1}\\right)\\right]\n",
    "\\end{align}\n",
    "\n",
    "Optimality conditions:\n",
    "\\begin{align}\n",
    "\\text{Intertemporal condition }\\left(c+a\\right)\t& :\t\\frac{1}{c_{t}}=\\beta\\mathrm{E}_{t}\\left[\\frac{1}{c_{t+1}}\\left(1+r_{t+1}\\right)\\right]\\\\\n",
    "\\text{Intratemporal condition }\\left(c+h\\right)\t& :\t\\frac{1}{c_{t}} =\\frac{\\phi}{w_{t}\\left(1-h_{t}\\right)}\\to\\quad h_{t}=1-\\phi\\frac{c_{t}}{w_{t}}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Firms\n",
    "A representative firm solves profit (dividend) maximization problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max\\quad&d_{t}=y_{t}-r_{t}^{k}k_{t}-w_{t}h_{t}\\\\\\text{subject to}\\quad&y_{t}=z_{t}k_{t}^{\\alpha}h_{t}^{1-\\alpha}\\\\&r_{t}^{k}=r_{t}+\\delta\n",
    "\\end{align}\n",
    "\n",
    "where:\n",
    "\n",
    "$d$\t- per capita dividends<br>\n",
    "$y$\t- per capita output<br>\n",
    "$r^{k}$\t- capital rental rate<br>\n",
    "$k$\t- per capita physical capital stock<br>\n",
    "$w$\t- real wage per hour<br>\n",
    "$h$\t- per capita hours (as fraction of total available time)<br>\n",
    "$z$\t- stochastic total factor productivity (TFP) level<br>\n",
    "$\\alpha$\t- physical capital share in output<br>\n",
    "$r$\t- real interest rate<br>\n",
    "$\\delta$\t- physical capital depreciation rate<br>\n",
    "\n",
    "Rewritten problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max\\quad d_{t}=z_{t}k_{t}^{\\alpha}h_{t}^{1-\\alpha}-\\left(r_{t}+\\delta\\right)k_{t}-w_{t}h_{t}\n",
    "\\end{align}\n",
    "\n",
    "First Order Conditions:\n",
    "\\begin{align}\n",
    "\\frac{\\partial d_{t}}{\\partial k_{t}}\t=\\alpha z_{t}k_{t}^{\\alpha-1}h_{t}^{1-\\alpha}-\\left(r_{t}+\\delta\\right)=0\\quad\\,\\to\\quad r_{t}=\\alpha z_{t}k_{t}^{\\alpha-1}h_{t}^{1-\\alpha}-\\delta\\\\\n",
    "\\frac{\\partial d_{t}}{\\partial h_{t}}\t=\\left(1-\\alpha\\right)z_{t}k_{t}^{\\alpha}h_{t}^{-\\alpha}-w_{t}=0\\quad\\quad\\to\\quad w_{t}=\\left(1-\\alpha\\right)z_{t}k_{t}^{\\alpha}h_{t}^{-\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "Alternative expressions for factor prices:\n",
    "\\begin{align}\n",
    "r_{t}\t=\\alpha\\frac{y_{t}}{k_{t}}-\\delta\\\\\n",
    "w_{t}\t=\\left(1-\\alpha\\right)\\frac{y_{t}}{h_{t}}\n",
    "\\end{align}\n",
    "\n",
    "Due to perfect competition and CRS economic profits equal zero:\n",
    "\\begin{align}\n",
    "d_{t}=y_{t}-r_{t}^{k}k_{t}-w_{t}h_{t}=y_{t}-\\alpha\\frac{y_{t}}{k_{t}}\\cdot k_{t}-\\left(1-\\alpha\\right)\\frac{y_{t}}{h_{t}}\\cdot h_{t}=0\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General Equilibrium\n",
    "Capital market clears:\n",
    "\\begin{align}\n",
    "a_{t}=k_{t}\n",
    "\\end{align}\n",
    "\n",
    "Households' budget constraint can be written as resource constraint:\n",
    "\\begin{align}\n",
    "a_{t+1}+c_{t}\t& =\t\\left(1+r_{t}\\right)a_{t}+w_{t}h_{t}+d_{t}\\\\\n",
    "k_{t+1}+c_{t}\t& =\t\\left(1+\\alpha\\frac{y_{t}}{k_{t}}-\\delta\\right)k_{t}+\\left(1-\\alpha\\right)\\frac{y_{t}}{h_{t}}\\cdot h_{t}+0\\\\\n",
    "k_{t+1}+c_{t}\t& =\t\\alpha y_{t}+\\left(1-\\delta\\right)k_{t}+\\left(1-\\alpha\\right)y_{t}\\\\\n",
    "k_{t+1}+c_{t}\t& =\ty_{t}+\\left(1-\\delta\\right)k_{t}\n",
    "\\end{align}\n",
    "\n",
    "If we define investment as:\n",
    "\\begin{align}\n",
    "i_{t}=k_{t+1}-\\left(1-\\delta\\right)k_{t}\n",
    "\\end{align}\n",
    "\n",
    "We can rewrite the resource constraint as the GDP accounting equation:\n",
    "\\begin{align}\n",
    "y_{t}=c_{t}+i_{t}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TFP Process\n",
    "TFP evolves according to an AR(1) process:\n",
    "\\begin{align}\n",
    "z_{t}=\\left(1-\\rho_{z}\\right)+\\rho_{z}\\cdot z_{t-1}+\\varepsilon_{t}\n",
    "\\end{align}\n",
    "\n",
    "where $\\rho_{z}<1$ regulates shock persistence and $\\varepsilon$ is zero-mean white noise\n",
    "\n",
    "It is often assumed that \n",
    "\\begin{align}\n",
    "\\varepsilon\\sim\\mathcal{N}\\left(0,\\sigma_{z}^{2}\\right)\n",
    "\\end{align}\n",
    "\n",
    "In the absence of shocks\n",
    "\\begin{align}\n",
    "z\\to1\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full set of equilibrium conditions\n",
    "System of 8 equations and 8 unknowns:\n",
    "$\\left\\{ y,c,i,k,h,w,r,z\\right\\}$\n",
    "\n",
    "\\begin{align}\n",
    "\\text{Euler equation}\t& :\t\\frac{1}{c_{t}}=\\beta\\mathrm{E}_{t}\\left[\\frac{1}{c_{t+1}}\\left(1+r_{t+1}\\right)\\right]\\\\\n",
    "\\text{Consumption-hours choice}\t& :\th_{t}=1-\\phi\\frac{c_{t}}{w_{t}}\\\\\n",
    "\\text{Production function}\t& :\ty_{t}=z_{t}k_{t}^{\\alpha}h_{t}^{1-\\alpha}\\\\\n",
    "\\text{Real interest rate}\t& :\tr_{t}=\\alpha\\frac{y_{t}}{k_{t}}-\\delta\\\\\n",
    "\\text{Real hourly wage}\t& :\tw_{t}=\\left(1-\\alpha\\right)\\frac{y_{t}}{h_{t}}\\\\\n",
    "\\text{Investment}\t& :\ti_{t}=k_{t+1}-\\left(1-\\delta\\right)k_{t}\\\\\n",
    "\\text{Output accounting}\t& :\ty_{t}=c_{t}+i_{t}\\\\\n",
    "\\text{TFP AR(1) process}\t& :\tz_{t}=\\left(1-\\rho_{z}\\right)+\\rho_{z}\\cdot z_{t-1}+\\varepsilon_{t}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution Method\n",
    "\n",
    "• Our model is a system of non-linear difference equations\n",
    "\n",
    "• There exist no closed form solutions for the transitional dynamics except for few unrealistic cases\n",
    "\n",
    "• We can solve easily an approximated version of the system\n",
    "\n",
    "– (log-)linearize by hand around steady-state\n",
    "\n",
    "– let Dynare compute n-th order Taylor expansion\n",
    "\n",
    "• Solving the DSGE model involves transforming the forward looking system into a VAR (backward looking) system\n",
    "\n",
    "– Many good methods: Blanchard-Kahn, Klein, Sims, etc.\n",
    "\n",
    "• Computer software exists that does it for you, e.g. DynareDynare does that for you: stoch_simul;\n",
    "\n",
    "• This is possible thanks to the Rational Expectations assumption"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Steady-State\n",
    "\n",
    "Start with the Euler equation:<br>\n",
    "$\\frac{1}{c}=\\beta\\frac{1}{c}\\left(1+r\\right)\\quad\\to\\quad r=\\frac{1}{\\beta}-1$\n",
    "\n",
    "From the interest rate equation obtain the k/h ratio:<br> \n",
    "$r=\\alpha k^{\\alpha-1}h^{1-\\alpha}-\\delta\\quad\\to\\quad\\frac{k}{h}=\\left(\\frac{\\alpha}{r+\\delta}\\right)^{\\frac{1}{1-\\alpha}}$\n",
    "\n",
    "From the production function obtain the y/h ratio and then wage:<br>\n",
    "$y=k^{\\alpha}h^{1-\\alpha}\\quad\\to\\quad\\frac{y}{h}=\\left(\\frac{k}{h}\\right)^{\\alpha}\\quad\\text{and}\\quad w=\\left(1-\\alpha\\right)\\frac{y}{h}$\n",
    "\n",
    "From investment and output accounting eqns. obtain the c/h ratio:<br>\n",
    "$i=\\delta k\\quad\\to\\quad y=c+\\delta k\\quad\\to\\quad\\frac{c}{h}=\\frac{y}{h}-\\delta\\frac{k}{h}$\n",
    "\n",
    "Get h from the consumption-hours choice. The rest follows from h:<br>\n",
    "$h=1-\\phi\\frac{c}{w}\\quad\\to\\quad1=\\frac{1}{h}-\\phi\\frac{c}{h}\\frac{1}{w}\\quad\\to\\quad h=1/\\left[1+\\phi\\frac{c}{h}\\frac{1}{w}\\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log-linearized Equilibrium\n",
    "\n",
    "Log-linearization means taking the log-deviation around a steady state value. Define the log-deviation of variable $x_{t}$ from its steady state $x^{*}$ as\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{x}_{t}=\\ln x_{t}-\\ln x^{*}\\approx\\dfrac{x_{t}-x^{*}}{x^{*}}\n",
    "\\end{align}\n",
    "\n",
    "The log-linearized system is given by \n",
    "\n",
    "\\begin{align}\n",
    "\\text{Euler equation}\t& :\t\\hat{c}_{t}=\\mathrm{E}_{t}\\hat{c}_{t+1}-\\dfrac{1}{1+r^{*}}\\mathrm{E}_{t}\\left[\\hat{r}_{t+1}\\right]\\\\\n",
    "\\text{Consumption-hours choice} &\t:\t\\hat{h}_{t}=\\dfrac{\\phi c^{*}}{w^{*}h^{*}}\\left(\\hat{w}_{t}-\\hat{c}_{t}\\right)\\\\\n",
    "\\text{Production function} &\t:\t\\hat{y}_{t}=\\hat{z}_{t}+\\alpha\\hat{k}_{t}+\\left(1-\\alpha\\right)\\hat{h}_{t}\\\\\n",
    "\\text{Real interest rate}\t& :\t\\hat{r}_{t}=\\alpha\\dfrac{y^{*}}{k^{*}}\\left(\\hat{y}_{t}-\\hat{k}_{t}\\right)\\\\\n",
    "\\text{Real hourly wage}\t& :\t\\hat{w}_{t}=\\hat{y}_{t}-\\hat{h}_{t}\\\\\n",
    "\\text{Investment}\t& :\t\\hat{k}_{t+1}=\\dfrac{i^{*}}{k^{*}}\\hat{i}_{t}+\\left(1-\\delta\\right)\\hat{k}_{t}\\\\\n",
    "\\text{Output accounting}\t& :\t\\hat{y}_{t}=\\dfrac{c^{*}}{y^{*}}\\hat{c}_{t}+\\dfrac{i^{*}}{y^{*}}\\hat{i}_{t}\\\\\n",
    "\\text{TFP AR(1) process}\t& :\t\\hat{z}_{t}=\\rho_{z}\\hat{z}_{t-1}+\\varepsilon_{t}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameter values\n",
    "• We need to specify parameter values tell Dynare the parameter values, e.g. $\\alpha = 0.33$;\n",
    "\n",
    "• There is a variety of approaches on how to obtain those values\n",
    "\n",
    "• Two most widely used are:\n",
    "\n",
    "– Calibration – picking parameter values to fit certain long-run (average) features of data. For example, we might want to pick the parameters so that the model's investment share in GDP matches the average share in the data\n",
    "\n",
    "– Estimation – Dynare allows us to easily run a Bayesian estimation procedure on real data. It still needs as an input prior estimates of parameter values and their confidence intervals, which makes the calibration exercise very useful\n",
    "\n",
    "• Most models in recent papers are estimated\n",
    "\n",
    "• Today's toy model is calibrated\n",
    "\n",
    "|Parameter | Value justification\t|Mean\n",
    "|:---|:---|:---|\n",
    "$\\alpha$ |\tCapital income share of GDP\t| 0.33\t\n",
    "$\\beta$\t| From average real interest rate\t| 0.99\t\n",
    "$\\delta$ |\tFrom investment share of GDP\t| 0.025\t\n",
    "$\\phi$\t| Work for 1/3 of time endowment\t| 1.75\t\n",
    "$\\rho_{z}$\t| Coefficient in TFP AR(1) regression\t| 0.97\t\n",
    "$\\sigma_{z}$\t| Error term in TFP AR(1) regression\t| 0.007\t\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full set of equilibrium conditions\n",
    "System of 8 equations and 8 unknowns:\n",
    "$\\left\\{ y,c,i,k,h,w,r,z\\right\\}$\n",
    "\n",
    "\\begin{align}\n",
    "\\text{Euler equation}\t& :\t\\frac{1}{c_{t}}=\\beta\\mathrm{E}_{t}\\left[\\frac{1}{c_{t+1}}\\left(1+r_{t+1}\\right)\\right]\\\\\n",
    "\\text{Consumption-hours choice}\t& :\th_{t}=1-\\phi\\frac{c_{t}}{w_{t}}\\\\\n",
    "\\text{Production function}\t& :\ty_{t}=z_{t}k_{t}^{\\alpha}h_{t}^{1-\\alpha}\\\\\n",
    "\\text{Real interest rate}\t& :\tr_{t}=\\alpha\\frac{y_{t}}{k_{t}}-\\delta\\\\\n",
    "\\text{Real hourly wage}\t& :\tw_{t}=\\left(1-\\alpha\\right)\\frac{y_{t}}{h_{t}}\\\\\n",
    "\\text{Investment}\t& :\ti_{t}=k_{t+1}-\\left(1-\\delta\\right)k_{t}\\\\\n",
    "\\text{Output accounting}\t& :\ty_{t}=c_{t}+i_{t}\\\\\n",
    "\\text{TFP AR(1) process}\t& :\tz_{t}=\\left(1-\\rho_{z}\\right)+\\rho_{z}\\cdot z_{t-1}+\\varepsilon_{t}\n",
    "\\end{align}"
   ]
  },
  {
   "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.9598499596182404}\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": {},
   "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.008016122502811389)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Standard Deviations')\n",
    "print(round(hp_cycles.std(),4))\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', round(hp_cycles[dta.columns.values[i]].autocorr(),4))\n",
    "\n",
    "print('')\n",
    "print('Correlations')\n",
    "print(round(hp_cycles.corr(method='pearson'),4))"
   ]
  },
  {
   "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.9598499596182404\n",
    "stderr_z = 0.008016122502811389\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 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}