Lab_4_NK.ipynb

{
 "cells": [
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Taylor Rule\n",
    "The Taylor rule is an interest rate forecasting model to determine what interest rates should be in order to shift the economy toward stable prices and full employment proposed by John Taylor in 1993.\n",
    "The original Taylor Rule is in the following form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "i_{t}=\\pi_{t}+0.5x_{t}+0.5*(\\pi_{t}-2)+2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $i_{t}$ is the federal funds rate, $\\pi_{t}$ is inflation, and $x_{t}$ is the output gap, actual output less potential output, the amount the economy can sustainably produce when capital and labor are fully employed.\n",
    "\n",
    "The (original) Taylor rule predicts that the FOMC will raise the federal funds rate (tighten monetary policy) by one-half percentage point:\n",
    "\n",
    "(1) for each percentage point that inflation rises relative to the Fed’s target, assumed to be 2 percent; or\n",
    "\n",
    "(2) for each percentage point that that output rises relative to its potential.\n",
    "\n",
    "The Taylor rule also predicts that when inflation is at target and output is at potential (the output gap is zero), the FOMC will set the real federal funds rate at 2 percent—about its historical average."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's see how original Taylor fits to the actual federal funds rate. We will download following variables from Federal Reserve Bank of St. Louis Economic Database (FRED) https://fred.stlouisfed.org/.\n",
    "\n",
    "- nominal Gross Domestic Product in billions of current dollars - `GDP`\n",
    "- nominal Potential Gross Domestic Product in billions of current dollars - `NGDPPOT`\n",
    "- Implicit Price Deflator for Gross Domestic Product - `GDPDEF`\n",
    "- Personal Consumption Expenditures: Chain-type Price Index - `PCEPI`\n",
    "- Core PCE, Personal Consumption Expenditures Excluding Food and Energy - `PCEPILFE`\n",
    "- Federal Funds Effective Rate - `FEDFUNDS`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fred = DataReader(['GDP', 'NGDPPOT', 'GDPDEF','PCEPI','PCEPILFE','FEDFUNDS'],\n",
    "                  'fred', start='1993-01', end='2023-12')\n",
    "fred = fred.resample('QS').mean()\n",
    "fred = fred.dropna()\n",
    "fred.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we are constructing the series that we need to test original Taylor rule:\n",
    "- `Output` is real GDP \n",
    "- `Potential Output` is the real potential GDP\n",
    "- `Output Gap` is the log real GDP less log real potential GDP\n",
    "- `Inflation` is annualized quarterly percentage change in GDP implicit price deflator\n",
    "- `PCE Inflation` is annualized quarterly percentage change in PCE\n",
    "- `Core PCE Inflation` is annualized quarterly percentage change in Core PCE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dta = pd.DataFrame()\n",
    "\n",
    "dta['Output'] = np.log(fred['GDP']*10**9/fred['GDPDEF']*100)\n",
    "dta['Potential Output'] = np.log(fred['NGDPPOT']*10**9/fred['GDPDEF']*100)\n",
    "dta['Output_Gap'] = 400*(dta['Output']-dta['Potential Output'])\n",
    "dta['Inflation'] = 400*fred['GDPDEF'].pct_change()\n",
    "dta['PCE_Inflation'] = 400*fred['PCEPI'].pct_change()\n",
    "dta['Core_PCE_Inflation'] = 400*fred['PCEPILFE'].pct_change()\n",
    "dta['FFR'] = fred['FEDFUNDS']\n",
    "\n",
    "dta.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the original Taylor formula to predict the interest rate. We will use three different inflation index to calculate the predicted FFR:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta['Taylor FFR PCE']     = dta['PCE_Inflation'] + 0.5*dta['Output_Gap'] + 0.5*(dta['PCE_Inflation']-2) + 2\n",
    "dta['Taylor FFR COREPCE'] = dta['Core_PCE_Inflation'] + 0.5*dta['Output_Gap'] + 0.5*(dta['Core_PCE_Inflation']-2) + 2\n",
    "\n",
    "fig, axs = plt.subplots(2, 1, figsize=(20, 15), sharex=False, sharey=False)\n",
    "dta[['FFR','Taylor FFR PCE']].plot(ax=axs[0], lw=2, style=['k-','r--'])\n",
    "dta[['FFR','Taylor FFR COREPCE']].plot(ax=axs[1], lw=2, style=['k-','r--'])\n",
    "\n",
    "axs[0].legend([\"FED Funds Rate\", \"Taylor Rule using PCE\"]);\n",
    "axs[1].legend([\"FED Funds Rate\", \"Taylor Rule using Core PCE\"]);\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that original Taylor rule does not predict well the behaviour of the interest rates. Now, let's estimate the following Taylor rule:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "    i_{t} = \\rho_{i}i_{t-1} + (1-\\rho_{i})*(\\phi_{\\pi}(\\pi_{t}-\\pi^{*})+\\phi_{x}x_{t})\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Taylor_PCE = smf.ols(formula='FFR ~ 0 + FFR.shift(1) + PCE_Inflation + Output_Gap', data=dta).fit()\n",
    "Taylor_COREPCE = smf.ols(formula='FFR ~ 0 + FFR.shift(1) + Core_PCE_Inflation + Output_Gap', data=dta).fit()\n",
    "\n",
    "Taylor_COREPCE.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's recover the parameter values from Taylor_COREPCE model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_pi=Taylor_COREPCE.params['Core_PCE_Inflation']/(1-Taylor_COREPCE.params['FFR.shift(1)'])\n",
    "phi_x=Taylor_COREPCE.params['Output_Gap']/(1-Taylor_COREPCE.params['FFR.shift(1)'])\n",
    "\n",
    "print('')\n",
    "print('phi_pi    =', phi_pi)\n",
    "print('phi_x    =', phi_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta['Taylor_PCE_fitted'] = Taylor_PCE.fittedvalues\n",
    "dta['Taylor_COREPCE_fitted'] = Taylor_COREPCE.fittedvalues\n",
    "\n",
    "fig, axs = plt.subplots(2, 1, figsize=(20, 15), sharex=False, sharey=False)\n",
    "dta[['FFR','Taylor_PCE_fitted']].plot(ax=axs[0], lw=2, style=['k-','r--'])\n",
    "dta[['FFR','Taylor_COREPCE_fitted']].plot(ax=axs[1], lw=2, style=['k-','r--'])\n",
    "\n",
    "axs[0].legend([\"FED Funds Rate\", \"Taylor Rule using PCE\"]);\n",
    "axs[1].legend([\"FED Funds Rate\", \"Taylor Rule using Core PCE\"]);\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3-equation New Keynesian Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Economy\n",
    "The representative household maximizes the following expected\n",
    "discounted sum of utilities over possible paths of consumption and\n",
    "labor:\n",
    "\n",
    "\\begin{align}\n",
    " & U_{t}=\\mathrm{E}_{t}\\sum_{k=0}^{\\infty}\\beta^{k}\\left[\\frac{C_{t+k}^{1-\\sigma}}{1-\\sigma}-\\phi\\frac{N_{t+k}^{1+\\eta}}{1+\\eta}\\right]\n",
    "\\end{align}\n",
    "\n",
    "The nominal period budget constraint is:\n",
    "\n",
    "\\begin{equation}\n",
    "P_{t}C_{t}+B_{t}=W_{t}N_{t}+\\left(1+i_{t-1}\\right)B_{t-1}+P_{t}D_{t}\n",
    "\\end{equation}\n",
    "\n",
    "The final output good is a CES aggregate of a continuum of intermediates:\n",
    "\n",
    "\\begin{equation}\n",
    "Y_{t}=\\left[\\intop_{0}^{1}Y_{t}\\left(i\\right)^{\\tfrac{\\varepsilon-1}{\\varepsilon}}di\\right]^{\\tfrac{\\varepsilon}{\\varepsilon-1}}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Each monopolistic firm $i$ produces a differentiated variety and\n",
    "hires a homogenous type of labor according to the linear production\n",
    "function:\n",
    "\n",
    "\\begin{align}\n",
    "Y_{t}\\left(i\\right)=A_{t}N_{t}\\left(i\\right)^{1-\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "Following Calvo, we assume now that firms adjust their price infrequently and that\n",
    "the opportunity to adjust follows an exogenous Poisson process. Each\n",
    "period there is a constant probability ($1-\\theta$) that the firm\n",
    "will be able to adjust its price, independently of past history. The\n",
    "expected time between price adjustments is therefore $\\dfrac{1}{1-\\theta}$.\n",
    "If the law of large numbers holds this implies that the fraction of\n",
    "retailers not setting prices at $t$ is $\\theta$. The draw is independent\n",
    "of history, so that we do not need to keep track of firms changing\n",
    "prices over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solutions to the above household and firm problems give the following three key equations of the New Keynesian model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New Keynesian Phillips Curve\n",
    "\\begin{align}\n",
    "\\pi_{t} & =\\beta\\mathrm{E}_{t}\\pi_{t+1}+\\kappa x_{t}\n",
    "\\end{align}\n",
    "\n",
    "In the above equation:\n",
    "\\begin{align}\n",
    "\\text{Output gap:} \\quad & x_{t}=y_{t}-y_{t}^{n}\\\\\n",
    "\\text{Output gap coefficient:} \\quad & \\kappa=\\frac{\\left(1-\\theta\\right)\\left(1-\\beta\\theta\\right)}{\\theta}\\left(\\sigma+\\eta\\right)\n",
    "\\end{align}\n",
    "\n",
    "The actual output and natural level of output are equal to:\n",
    "\\begin{align}\n",
    "\\text{Actual output:} \\quad & y_{t} & =a_{t}+(1-\\alpha)n_{t}\\\\\n",
    "\\text{Natural level of output} \\quad & y_{t}^{n} & =\\dfrac{1+\\eta}{\\sigma(1-\\alpha)+\\eta+\\alpha}a_{t}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New Keynesian IS\n",
    "\\begin{align}\n",
    "x_{t} & =\\mathrm{E}_{t}x_{t+1}-\\frac{1}{\\sigma}\\left(i_{t}-\\mathrm{E}_{t}\\pi_{t+1}-r_{t}^{n}\\right)\n",
    "\\end{align}\n",
    "\n",
    "In the above equation:\n",
    "\\begin{align}\n",
    "\\text{Natural real interest rate:} \\quad & r_{t}^{n} & =\\rho+\\sigma\\left(\\dfrac{1+\\eta}{\\sigma(1-\\alpha)+\\eta+\\alpha}\\mathrm{E}_{t}\\triangle a_{t+1}\\right)\n",
    "\\end{align}\n",
    "\n",
    "The technology evolves according to:\n",
    "\\begin{align}\n",
    "\\text{TFP:} \\quad & a_{t}=\\rho_{a}a_{t-1}+\\varepsilon_{a,t}\n",
    "\\end{align}\n",
    "\n",
    "The devation of natural real interest rate from its steady-state is given by:\n",
    "\\begin{align}\n",
    "\\hat{r}_{t}^{n} & =-\\sigma\\left(\\dfrac{1+\\eta}{\\sigma(1-\\alpha)+\\eta+\\alpha}\\right)(1-\\rho_{a})a_{t}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Monetary Policy Rule\n",
    "\\begin{align}\n",
    "i_{t}=\\phi_{\\pi}\\pi_{t}+\\phi_{x}x_{t}+\\nu_{t}\n",
    "\\end{align}\n",
    "\n",
    "The monetary policy shock evolves according to:\n",
    "\\begin{align}\n",
    "\\text{Monetary Shock:} \\quad & \\nu_{t}=\\rho_{\\nu}\\nu_{t-1}+\\varepsilon_{R,t}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above equations, the parameters correspond to:\n",
    "\n",
    "|Parameter | Justification\t\n",
    "|:---|:---|\n",
    "$\\alpha$ |\tShare of capital in production (output)\t\n",
    "$\\beta$\t| Household discount factor\t\t\n",
    "$\\sigma$ |\tRisk aversion (Inverse of intertemporal elasticity of substitution)\t\n",
    "$\\eta$\t| Inverse Frisch-elasticity\t\t\n",
    "$\\theta$\t| quarterly probability of price rigidity\t\t\n",
    "$\\kappa$\t| Coefficient of output gap in NKPC\t\t\n",
    "$\\phi_{\\pi}$\t| Response of monetary rule to inflation\t\t\n",
    "$\\phi_{x}$\t| Response of monetary rule to output gap\t\n",
    "$\\rho_{a}$\t| Coefficient in TFP AR(1) regression\t\t\n",
    "$\\rho_{v}$\t| Coefficient of monetary policy shock AR(1) regression\t\t\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Dynare import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = 'Pi x i y r_r n p nu w y_n r_n a'\n",
    "varexo = 'eps_R'\n",
    "\n",
    "param_values = {'alppha':0.25, sy.symbols('beta'):0.99, 'siggma':2, 'eta':3.77, 'rho_nu':0.5,\n",
    "                'theta':0.75, 'kappa':0.1717, 'phi_pi':1.25, 'phi_x':0.125, 'rho_a':0.9}\n",
    "\n",
    "model = ('-Pi + betta*Pi(+1)+kappa*x',\n",
    "         '-x + x(+1) -1/siggma*(i-Pi(+1)-r_n)',\n",
    "         '-i + phi_pi*Pi+phi_x*x+nu',\n",
    "         '-r_n + -siggma*(1+eta)/(siggma*(1-alppha)+eta+alppha)*(1-rho_a)*a',\n",
    "         '-r_r + i-Pi(+1)',\n",
    "         '-y_n + (1+eta)/(siggma*(1-alppha)+eta+alppha)*a',\n",
    "         '-x + y-y_n',\n",
    "         '-nu + rho_nu*nu(-1)+eps_R',\n",
    "         '-a + rho_a*a(-1)',\n",
    "         '-y + a+(1-alppha)*n',\n",
    "         '-Pi + p-p(-1)',\n",
    "         '-w+p+siggma*y+eta*n')\n",
    "\n",
    "initval = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbc = Dynare(var, varexo, param_values, model, initval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbc.steady()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbc.stoch_simul(irf=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework\n",
    "The other important shocks in the New Keynesian Model are demand and cost push-up shocks. We can incorporate them into the 3-eq New Keynesian model in following way.\n",
    "\n",
    "The New Keynesian IS (NKIS) with demand shock where demand shock follows AR(1) process is given by\n",
    "\n",
    "\\begin{align}\n",
    "x_{t} & =\\mathrm{E}_{t}x_{t+1}-\\frac{1}{\\sigma}\\left(i_{t}-\\mathrm{E}_{t}\\pi_{t+1}-r_{t}^{n}\\right)+u_{t}\\\\\n",
    "u_{t} & = \\rho_{u}u_{t-1}+\\epsilon_{u,t}\n",
    "\\end{align}\n",
    "\n",
    "The New Keynesian Phillips Curve (NKPC) with cost push-up shock where cost push-up shock follows AR(1) process is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\pi_{t} & =\\beta\\mathrm{E}_{t}\\pi_{t+1}+\\kappa x_{t}+e{t}\\\\\n",
    "e_{t} & = \\rho_{e}e_{t-1}+\\epsilon_{e,t}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1**: Insert the **demand shock** given in the above **New Keynesian IS equation** into the New Keynesian Model. Compute impulse response functions (IRFS) of the variables following a **demand shock**. Briefly comment on the responses of inflation, output gap, and nominal interest rate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = 'Pi x i y r_r n p u w y_n r_n a'\n",
    "varexo = 'eps_u'\n",
    "\n",
    "param_values = {'alppha':0.25, sy.symbols('beta'):0.99, 'siggma':2, 'eta':3.77, 'rho_u':0.75,\n",
    "                'theta':0.75, 'kappa':0.1717, 'phi_pi':1.25, 'phi_x':0.125, 'rho_a':0.9}\n",
    "\n",
    "model = ('-Pi + betta*Pi(+1)+kappa*x',\n",
    "#         'Here goes the New Keynesian IS curve with demand shock',\n",
    "         '-i + phi_pi*Pi+phi_x*x',\n",
    "         '-r_n + -siggma*(1+eta)/(siggma*(1-alppha)+eta+alppha)*(1-rho_a)*a',\n",
    "         '-r_r + i-Pi(+1)',\n",
    "         '-y_n + (1+eta)/(siggma*(1-alppha)+eta+alppha)*a',\n",
    "         '-x + y-y_n',\n",
    "#         'Here goes the AR(1) demand shock',\n",
    "         '-a + rho_a*a(-1)',\n",
    "         '-y + a+(1-alppha)*n',\n",
    "         '-Pi + p-p(-1)',\n",
    "         '-w+p+siggma*y+eta*n')\n",
    "\n",
    "initval = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2**: Insert the **cost push-up shock** given in the above **NKPC equation** into the New Keynesian Model. Compute impulse response functions (IRFS) of the variables following a **cost push-up shock**. Briefly comment on the responses of inflation, output gap, and nominal interest rate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = 'Pi x i y r_r n p e w y_n r_n a'\n",
    "varexo = 'eps_e'\n",
    "\n",
    "param_values = {'alppha':0.25, sy.symbols('beta'):0.99, 'siggma':2, 'eta':3.77, 'rho_e':0.75,\n",
    "                'theta':0.75, 'kappa':0.1717, 'phi_pi':1.25, 'phi_x':0.125, 'rho_a':0.9}\n",
    "\n",
    "model = (#         'Here goes the NKPC with cost push-up shock',\n",
    "         '-x +x(+1) -1/siggma*(i-Pi(+1)-r_n)',\n",
    "         '-i + phi_pi*Pi+phi_x*x',\n",
    "         '-r_n + -siggma*(1+eta)/(siggma*(1-alppha)+eta+alppha)*(1-rho_a)*a',\n",
    "         '-r_r + i-Pi(+1)',\n",
    "         '-y_n + (1+eta)/(siggma*(1-alppha)+eta+alppha)*a',\n",
    "         '-x + y-y_n',\n",
    "#         'Here goes the AR(1) cost push-up shock',\n",
    "         '-a + rho_a*a(-1)',\n",
    "         '-y + a+(1-alppha)*n',\n",
    "         '-Pi + p-p(-1)',\n",
    "         '-w+p+siggma*y+eta*n')\n",
    "\n",
    "initval = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optional: Unemployment and Labor Markets over the Business Cycle\n",
    "\n",
    "We will take a look at the data on unemployment and functioning of labor markets in the United States"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we are downloading two datasets. \n",
    "\n",
    "The first dataset, `fred`, contains monthly data on: \n",
    "\n",
    "- whether in a given month the US economy was in a recession state (`1`) or not (`0`) - `USREC`\n",
    "- number of people in labor force in thousands - `CLF16OV`\n",
    "- number of unemployed people in thousands - `UNEMPLOY`\n",
    "- unemployment rate in percent - `UNRATE`\n",
    "- number of job openings (vacancies) in thousands - `JTSJOL`\n",
    "- job vacancy rate - `JTSJOR`\n",
    "\n",
    "The second dataset, `hours`, contains quarterly data on:\n",
    "\n",
    "- real GDP in billions of 2009 dollars - `GDPC1`\n",
    "- total hours worked in the nonfarm business sector (index) - `HOANBS`\n",
    "- average hours worked per employee in the nonfarm business sector (index) - `PRS85006023`\n",
    "- number of employees in the nonfarm business sector (index) - `PRS85006013`"
   ]
  },
  {
   "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": [
    "# Get FRED data\n",
    "fred = DataReader(['USREC', 'CLF16OV', 'UNEMPLOY', 'UNRATE', 'JTSJOL', 'JTSJOR'], \n",
    "                   'fred', start=start, end=end)\n",
    "\n",
    "hours = DataReader(['GDPC1', 'HOANBS', 'PRS85006023', 'PRS85006013'], \n",
    "                    'fred', start=start, end=end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Separate trend and cyclical components of GDP, hours and employment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_cycle = pd.DataFrame()\n",
    "hp_trend = pd.DataFrame()\n",
    "\n",
    "cf_cycle = pd.DataFrame()\n",
    "cf_trend = pd.DataFrame()\n",
    "\n",
    "for col in hours.columns:\n",
    "    hp_cycle[col], hp_trend[col] = sm.tsa.filters.hpfilter(100*np.log(hours[col]).dropna(), lamb=1600)\n",
    "    cf_cycle[col], cf_trend[col] = sm.tsa.filters.cffilter(100*np.log(hours[col]).dropna(), low=6, high=32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare cyclical components of total hours worked vs its components: hours per employee and number of employees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_cycle.columns = ['Output','Total Hours','Hours per Employee','Employment']\n",
    "cf_cycle.columns = ['Output','Total Hours','Hours per Employee','Employment']\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))\n",
    "\n",
    "cf_cycle[['Total Hours','Employment']].to_period('D').plot(ax=ax1, style=['k-','r-'], lw=2)\n",
    "\n",
    "ylim = ax1.get_ylim()\n",
    "\n",
    "ax1.hlines(0, cf_cycle.index[0], cf_cycle.index[-1], linewidth=0.5)\n",
    "ax1.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "l = ax1.legend(loc='upper right')\n",
    "l.get_frame().set_linewidth(0)\n",
    "l.get_frame().set_alpha(1)\n",
    "\n",
    "cf_cycle[['Total Hours','Hours per Employee']].to_period('D').plot(ax=ax2, style=['k-','r-'], lw=2)\n",
    "\n",
    "ylim = ax2.get_ylim()\n",
    "\n",
    "ax2.hlines(0, cf_cycle.index[0], cf_cycle.index[-1], linewidth=0.5)\n",
    "ax2.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "l = ax2.legend(loc='upper right')\n",
    "l.get_frame().set_linewidth(0)\n",
    "l.get_frame().set_alpha(1)\n",
    "\n",
    "# plt.savefig('Hours_CF.pdf', bbox_inches='tight', pad_inches=0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the variance-covariance matrix of total hours worked and its components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_cycle[['Total Hours','Employment','Hours per Employee']].cov()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructing the vacancy rate time series\n",
    "\n",
    "The statistics on job openings (vacancies) from the `JOLTS` program are available only starting from December 2000. However, there are data on `Help Wanted Index`, which were gathered by private companies. Thanks to the work of Regis Barnichon, we can use them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta = fred.copy()\n",
    "dta.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dta[['JTSJOR','UNRATE']]['2000-12':].plot(lw=2)\n",
    "plt.legend(frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in Regis Barnichon's Composite Help Wanted Index and join the two datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hwi = pd.read_csv('data/HWI_index_old.txt', delimiter='\\t', skiprows=5)\n",
    "\n",
    "# Manage dates\n",
    "dates = []\n",
    "for date in hwi['Date ']:\n",
    "    dates.append(date[-2:]+'-01-'+date[0:4])\n",
    "\n",
    "hwi.index = pd.to_datetime(dates)\n",
    "hwi.index.rename('DATE', inplace=True)\n",
    "\n",
    "# Cleanup\n",
    "hwi = hwi.drop('Date ', 1)\n",
    "hwi.columns = ['HWI']\n",
    "hwi.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Join datasets\n",
    "df = dta.join(hwi)\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adjust the index to observed levels and splice the data from two sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Vacancies'] = df['JTSJOL']['2014-01-01'] * df['HWI'] / df['HWI']['2014-01-01']\n",
    "df['Vacancies']['2005-01-01':] = df['JTSJOL']['2005-01-01':]\n",
    "\n",
    "df[['Vacancies','JTSJOL']].plot(lw=2)\n",
    "plt.legend(frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct time series for unemployment and vacancy rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Unemployment rate'] = 100 * df['UNEMPLOY'] / df['CLF16OV']\n",
    "df['Vacancy rate'] = 100 * df['Vacancies'] / df['CLF16OV']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "df['Vacancy rate'].to_period('D').plot(ax=ax, style='k', lw=2)\n",
    "df['Unemployment rate'].to_period('D').plot(ax=ax, style='r', lw=2)\n",
    "\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "l = ax.legend(loc='upper left')\n",
    "l.get_frame().set_linewidth(0)\n",
    "l.get_frame().set_alpha(1)\n",
    "\n",
    "plt.title('US vacancy and unemployment rates (%)')\n",
    "# plt.savefig('VU.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Behavior of unemployment and vacancy rates in the United States\n",
    "\n",
    "Below I plot the scatterplot of unemployment and vacancy rates, with colors reflecting different decades. \n",
    "\n",
    "The resulting negative relationship is known as the Beveridge curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfq = df.resample('QS').mean()\n",
    "\n",
    "plt.plot(dfq['Unemployment rate']['1950-01-01':'1959-12-01'], \n",
    "         dfq['Vacancy rate']['1950-01-01':'1959-12-01'], 'o-', label='1950s')\n",
    "plt.plot(dfq['Unemployment rate']['1960-01-01':'1969-12-01'], \n",
    "         dfq['Vacancy rate']['1960-01-01':'1969-12-01'], 'o-', label='1960s')\n",
    "plt.plot(dfq['Unemployment rate']['1970-01-01':'1979-12-01'], \n",
    "         dfq['Vacancy rate']['1970-01-01':'1979-12-01'], 'o-', label='1970s')\n",
    "plt.plot(dfq['Unemployment rate']['1980-01-01':'1989-12-01'], \n",
    "         dfq['Vacancy rate']['1980-01-01':'1989-12-01'], 'o-', label='1980s')\n",
    "plt.plot(dfq['Unemployment rate']['1990-01-01':'1999-12-01'], \n",
    "         dfq['Vacancy rate']['1990-01-01':'1999-12-01'], 'o-', label='1990s')\n",
    "plt.plot(dfq['Unemployment rate']['2000-01-01':'2009-12-01'], \n",
    "         dfq['Vacancy rate']['2000-01-01':'2009-12-01'], 'o-', label='2000s')\n",
    "plt.plot(dfq['Unemployment rate']['2010-01-01':'2019-12-01'], \n",
    "         dfq['Vacancy rate']['2010-01-01':'2019-12-01'], 'ko-', label='2010s')\n",
    "\n",
    "plt.legend(frameon=False)\n",
    "plt.xlim(2, 12)\n",
    "plt.ylim(1, 5)\n",
    "plt.yticks(np.arange(1, 6))\n",
    "\n",
    "plt.xlabel('Unemployment rate (%)')\n",
    "plt.ylabel('Vacancy rate (%)')\n",
    "\n",
    "plt.title('Shifts in the US Beveridge curve')\n",
    "# plt.savefig('BC.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Separate trend from cycle to eliminate structural shifts to the Beveridge curve, note the adjustment of filtering options to monthly frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_cycle_uv = pd.DataFrame()\n",
    "hp_trend_uv = pd.DataFrame()\n",
    "\n",
    "cf_cycle_uv = pd.DataFrame()\n",
    "cf_trend_uv = pd.DataFrame()\n",
    "\n",
    "for col in ['Vacancy rate','Unemployment rate']:\n",
    "    hp_cycle_uv[col], hp_trend_uv[col] = sm.tsa.filters.hpfilter(100*np.log(df[col]).dropna(), lamb=1600*3**4)\n",
    "    cf_cycle_uv[col], cf_trend_uv[col] = sm.tsa.filters.cffilter(100*np.log(df[col]).dropna(), low=1.5*12, high=8*12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot cyclical components of unemployment and vacancy rates vs cyclical component of output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "cf_cycle_uv.resample('QS').mean().to_period('D').plot(ax=ax, style=['k','r'], lw=2)\n",
    "cf_cycle['Output'].plot(ax=ax, style=['b'], lw=2)\n",
    "\n",
    "ax.set_ylim(-60, 60)\n",
    "ylim = ax.get_ylim()\n",
    "\n",
    "ax.hlines(0, hours.index[0], hours.index[-1], linewidth=0.5)\n",
    "\n",
    "ax.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')\n",
    "\n",
    "ax.set_xlim('1950-01', hours.index[-1])\n",
    "\n",
    "l = ax.legend(loc='upper right')\n",
    "l.get_frame().set_linewidth(0)\n",
    "l.get_frame().set_alpha(1)\n",
    "\n",
    "plt.title('Deviations from Christiano-Fitzgerald trend (%)')\n",
    "# plt.savefig('VU_CF.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run a (very simplified) linear regression on cyclical components of unemployment and vacancy rates, the slope is very close to -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = hp_cycle_uv['Unemployment rate']\n",
    "y = hp_cycle_uv['Vacancy rate']\n",
    "\n",
    "slope, intercept = np.polyfit(x, y, 1)\n",
    "\n",
    "print(slope)\n",
    "print(intercept)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "plt.scatter(x, y, alpha=0.25) #facecolor='none', edgecolor='C0'\n",
    "plt.plot(x, intercept + slope*x, 'r-', lw=2)\n",
    "\n",
    "plt.xlim(-45, 45)\n",
    "plt.ylim(-45, 45)\n",
    "\n",
    "plt.hlines(0, -45, 45, linewidth=0.5)\n",
    "plt.vlines(0, -45, 45, linewidth=0.5)\n",
    "\n",
    "plt.title('Deviations from Hodrick-Prescott trend (%)')\n",
    "plt.xlabel('Unemployment rate')\n",
    "plt.ylabel('Vacancy rate')\n",
    "\n",
    "# plt.savefig('BC_HP.pdf')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the 'estimated' Beveridge curve without structural shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.mean(dfq['Unemployment rate'])\n",
    "v = np.mean(dfq['Vacancy rate'])\n",
    "\n",
    "print(u, v)\n",
    "\n",
    "scale = np.linspace(-40, 60, 101)\n",
    "\n",
    "plt.plot(u*np.exp(scale/100), v*np.exp(slope*scale/100), 'r', lw=2)\n",
    "\n",
    "plt.plot(u, v, 'ko')\n",
    "\n",
    "plt.xlim(2, 12)\n",
    "plt.ylim(1, 5)\n",
    "plt.yticks(np.arange(1, 6))\n",
    "\n",
    "# plt.hlines(v, 2, u, linestyle='--', lw=1)\n",
    "# plt.vlines(u, 1, v, linestyle='--', lw=1)\n",
    "\n",
    "plt.xlabel('Unemployment rate (%)')\n",
    "plt.ylabel('Vacancy rate (%)')\n",
    "plt.title('US Beveridge curve without structural shifts')\n",
    "\n",
    "# plt.savefig('BC_est.pdf')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.mean(dfq['Unemployment rate'])\n",
    "v = np.mean(dfq['Vacancy rate'])\n",
    "\n",
    "print(u, v)\n",
    "\n",
    "scale = np.linspace(-40, 60, 101)\n",
    "\n",
    "plt.plot(u*np.exp(scale/100), v*np.exp(slope*scale/100), 'r', lw=2)\n",
    "plt.scatter(u*np.exp(x/100), v*np.exp(y/100), alpha=0.25) #, marker='.'\n",
    "\n",
    "plt.plot(u, v, 'ko')\n",
    "\n",
    "plt.xlim(2, 12)\n",
    "plt.ylim(1, 5)\n",
    "plt.yticks(np.arange(1, 6))\n",
    "\n",
    "plt.hlines(v, 2, 12, linewidth=0.5)\n",
    "plt.vlines(u, 1, 5, linewidth=0.5)\n",
    "\n",
    "plt.xlabel('Unemployment rate (%)')\n",
    "plt.ylabel('Vacancy rate (%)')\n",
    "plt.title('US Beveridge curve without structural shifts')\n",
    "\n",
    "# plt.savefig('BC_est_2.pdf')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 2
}