{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\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": [
    "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": [
    "# Get FRED data\n",
    "fred = DataReader(['USREC','CLF16OV','UNEMPLOY','UNRATE','JTSJOL','JTSJOR'], \n",
    "                   'fred', start='1947-01', end='2020-01')\n",
    "\n",
    "hours = DataReader(['GDPC1','HOANBS','PRS85006023','PRS85006013'], \n",
    "                    'fred', start='1947-01', end='2020-01')"
   ]
  },
  {
   "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, hp_trend = sm.tsa.filters.hpfilter(100*np.log(hours[['GDPC1','HOANBS','PRS85006023','PRS85006013']]).dropna())\n",
    "cf_cycle, cf_trend = sm.tsa.filters.cffilter(100*np.log(hours[['GDPC1','HOANBS','PRS85006023','PRS85006013']]).dropna())"
   ]
  },
  {
   "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()\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()\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, hp_trend_uv = sm.tsa.filters.hpfilter(100*np.log(df[['Vacancy rate','Unemployment rate']]).dropna(), \n",
    "                                             lamb=1600*3**4)\n",
    "\n",
    "cf_cycle_uv, cf_trend_uv = sm.tsa.filters.cffilter(100*np.log(df[['Vacancy rate','Unemployment rate']]).dropna(), \n",
    "                                            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(-45, 45)\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",
    "plt.scatter(x, y, 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(intercept+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.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(intercept+slope*scale/100), 'r', lw=2)\n",
    "plt.scatter(u*np.exp(x/100), v*np.exp(y/100), marker='.')\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",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
