Unemployment and Labor Markets over the Business Cycle

We will take a look at the data on unemployment and functioning of labor markets in the United States

In [1]:
%matplotlib inline

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

from pandas_datareader.data import DataReader
C:\Users\Marcin\Anaconda3\lib\site-packages\statsmodels\compat\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools
In [2]:
from matplotlib import rcParams

# Restore old behavior of rounding default axis ranges
rcParams['axes.autolimit_mode'] = 'round_numbers'
rcParams['axes.xmargin'] = 0
rcParams['axes.ymargin'] = 0

Below we are downloading two datasets.

The first dataset, fred, contains monthly data on:

  • whether in a given month the US economy was in a recession state (1) or not (0) - USREC
  • number of people in labor force in thousands - CLF16OV
  • number of unemployed people in thousands - UNEMPLOY
  • unemployment rate in percent - UNRATE
  • number of job openings (vacancies) in thousands - JTSJOL
  • job vacancy rate - JTSJOR

The second dataset, hours, contains quarterly data on:

  • real GDP in billions of 2009 dollars - GDPC1
  • total hours worked in the nonfarm business sector (index) - HOANBS
  • average hours worked per employee in the nonfarm business sector (index) - PRS85006023
  • number of employees in the nonfarm business sector (index) - PRS85006013
In [3]:
# Get FRED data
fred = DataReader(['USREC','CLF16OV','UNEMPLOY','UNRATE','JTSJOL','JTSJOR'], 
                   'fred', start='1947-01', end='2020-01')

hours = DataReader(['GDPC1','HOANBS','PRS85006023','PRS85006013'], 
                    'fred', start='1947-01', end='2020-01')

Separate trend and cyclical components of GDP, hours and employment

In [4]:
hp_cycle, hp_trend = sm.tsa.filters.hpfilter(100*np.log(hours[['GDPC1','HOANBS','PRS85006023','PRS85006013']]).dropna())
cf_cycle, cf_trend = sm.tsa.filters.cffilter(100*np.log(hours[['GDPC1','HOANBS','PRS85006023','PRS85006013']]).dropna())

Compare cyclical components of total hours worked vs its components: hours per employee and number of employees

In [5]:
hp_cycle.columns = ['Output','Total Hours','Hours per Employee','Employment']
cf_cycle.columns = ['Output','Total Hours','Hours per Employee','Employment']

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))

cf_cycle[['Total Hours','Employment']].to_period('D').plot(ax=ax1, style=['k-','r-'], lw=2)

ylim = ax1.get_ylim()

ax1.hlines(0, cf_cycle.index[0], cf_cycle.index[-1], linewidth=0.5)
ax1.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

l = ax1.legend(loc='upper right')
l.get_frame().set_linewidth(0)
l.get_frame().set_alpha(1)

cf_cycle[['Total Hours','Hours per Employee']].to_period('D').plot(ax=ax2, style=['k-','r-'], lw=2)

ylim = ax2.get_ylim()

ax2.hlines(0, cf_cycle.index[0], cf_cycle.index[-1], linewidth=0.5)
ax2.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

l = ax2.legend(loc='upper right')
l.get_frame().set_linewidth(0)
l.get_frame().set_alpha(1)

# plt.savefig('Hours_CF.pdf', bbox_inches='tight', pad_inches=0.05)
plt.show()

Calculate the variance-covariance matrix of total hours worked and its components

In [6]:
hp_cycle[['Total Hours','Employment','Hours per Employee']].cov()
Out[6]:
Total Hours Employment Hours per Employee
Total Hours 3.520613 2.870760 0.649864
Employment 2.870760 2.465843 0.404910
Hours per Employee 0.649864 0.404910 0.244972

Constructing the vacancy rate time series

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.

In [7]:
dta = fred.copy()
dta.tail()
Out[7]:
USREC CLF16OV UNEMPLOY UNRATE JTSJOL JTSJOR
DATE
2017-12-01 0 160597.0 6576.0 4.1 5667.0 3.7
2018-01-01 0 161115.0 6684.0 4.1 6228.0 4.0
2018-02-01 0 161921.0 6706.0 4.1 6078.0 3.9
2018-03-01 0 161763.0 6585.0 4.1 6550.0 4.2
2018-04-01 0 161527.0 6346.0 3.9 NaN NaN
In [8]:
dta[['JTSJOR','UNRATE']]['2000-12':].plot()
plt.legend(frameon=False)
plt.show()

Read in Regis Barnichon's Composite Help Wanted Index and join the two datasets

In [9]:
hwi = pd.read_csv('data/HWI_index_old.txt', delimiter='\t', skiprows=5)

# Manage dates
dates = []
for date in hwi['Date ']:
    dates.append(date[-2:]+'-01-'+date[0:4])

hwi.index = pd.to_datetime(dates)
hwi.index.rename('DATE', inplace=True)

# Cleanup
hwi = hwi.drop('Date ', 1)
hwi.columns = ['HWI']
hwi.tail()
Out[9]:
HWI
DATE
2014-08-01 92.02
2014-09-01 87.76
2014-10-01 90.36
2014-11-01 91.35
2014-12-01 89.19
In [10]:
# Join datasets
df = dta.join(hwi)
df.tail()
Out[10]:
USREC CLF16OV UNEMPLOY UNRATE JTSJOL JTSJOR HWI
DATE
2017-12-01 0 160597.0 6576.0 4.1 5667.0 3.7 NaN
2018-01-01 0 161115.0 6684.0 4.1 6228.0 4.0 NaN
2018-02-01 0 161921.0 6706.0 4.1 6078.0 3.9 NaN
2018-03-01 0 161763.0 6585.0 4.1 6550.0 4.2 NaN
2018-04-01 0 161527.0 6346.0 3.9 NaN NaN NaN

Adjust the index to observed levels and splice the data from two sources

In [11]:
df['Vacancies'] = df['JTSJOL']['2014-01-01'] * df['HWI'] / df['HWI']['2014-01-01']
df['Vacancies']['2005-01-01':] = df['JTSJOL']['2005-01-01':]

df[['Vacancies','JTSJOL']].plot()
plt.legend(frameon=False)
plt.show()
C:\Users\Marcin\Anaconda3\lib\site-packages\ipykernel_launcher.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  

Construct time series for unemployment and vacancy rates

In [12]:
df['Unemployment rate'] = 100 * df['UNEMPLOY'] / df['CLF16OV']
df['Vacancy rate'] = 100 * df['Vacancies'] / df['CLF16OV']
In [13]:
fig, ax = plt.subplots()

df['Vacancy rate'].to_period('D').plot(ax=ax, style='k', lw=2)
df['Unemployment rate'].to_period('D').plot(ax=ax, style='r', lw=2)

ylim = ax.get_ylim()

ax.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

l = ax.legend(loc='upper left')
l.get_frame().set_linewidth(0)
l.get_frame().set_alpha(1)

plt.title('US vacancy and unemployment rates (%)')
# plt.savefig('VU.pdf')
plt.show()

Behavior of unemployment and vacancy rates in the United States

Below I plot the scatterplot of unemployment and vacancy rates, with colors reflecting different decades.

The resulting negative relationship is known as the Beveridge curve

In [14]:
dfq = df.resample('QS').mean()

plt.plot(dfq['Unemployment rate']['1950-01-01':'1959-12-01'], 
         dfq['Vacancy rate']['1950-01-01':'1959-12-01'], 'o-', label='1950s')
plt.plot(dfq['Unemployment rate']['1960-01-01':'1969-12-01'], 
         dfq['Vacancy rate']['1960-01-01':'1969-12-01'], 'o-', label='1960s')
plt.plot(dfq['Unemployment rate']['1970-01-01':'1979-12-01'], 
         dfq['Vacancy rate']['1970-01-01':'1979-12-01'], 'o-', label='1970s')
plt.plot(dfq['Unemployment rate']['1980-01-01':'1989-12-01'], 
         dfq['Vacancy rate']['1980-01-01':'1989-12-01'], 'o-', label='1980s')
plt.plot(dfq['Unemployment rate']['1990-01-01':'1999-12-01'], 
         dfq['Vacancy rate']['1990-01-01':'1999-12-01'], 'o-', label='1990s')
plt.plot(dfq['Unemployment rate']['2000-01-01':'2009-12-01'], 
         dfq['Vacancy rate']['2000-01-01':'2009-12-01'], 'o-', label='2000s')
plt.plot(dfq['Unemployment rate']['2010-01-01':'2019-12-01'], 
         dfq['Vacancy rate']['2010-01-01':'2019-12-01'], 'ko-', label='2010s')

plt.legend(frameon=False)
plt.xlim(2, 12)
plt.ylim(1, 5)
plt.yticks(np.arange(1, 6))

plt.xlabel('Unemployment rate (%)')
plt.ylabel('Vacancy rate (%)')

plt.title('Shifts in the US Beveridge curve')
# plt.savefig('BC.pdf')
plt.show()

Separate trend from cycle to eliminate structural shifts to the Beveridge curve, note the adjustment of filtering options to monthly frequency

In [15]:
hp_cycle_uv, hp_trend_uv = sm.tsa.filters.hpfilter(100*np.log(df[['Vacancy rate','Unemployment rate']]).dropna(), 
                                             lamb=1600*3**4)

cf_cycle_uv, cf_trend_uv = sm.tsa.filters.cffilter(100*np.log(df[['Vacancy rate','Unemployment rate']]).dropna(), 
                                            low=1.5*12, high=8*12)

Plot cyclical components of unemployment and vacancy rates vs cyclical component of output

In [16]:
fig, ax = plt.subplots()

cf_cycle_uv.resample('QS').mean().to_period('D').plot(ax=ax, style=['k','r'], lw=2)
cf_cycle['Output'].plot(ax=ax, style=['b'], lw=2)

ax.set_ylim(-45, 45)
ylim = ax.get_ylim()

ax.hlines(0, hours.index[0], hours.index[-1], linewidth=0.5)

ax.fill_between(fred.index, ylim[0], ylim[1], fred['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

ax.set_xlim('1950-01', hours.index[-1])

l = ax.legend(loc='upper right')
l.get_frame().set_linewidth(0)
l.get_frame().set_alpha(1)

plt.title('Deviations from Christiano-Fitzgerald trend (%)')
# plt.savefig('VU_CF.pdf')
plt.show()

Run a (very simplified) linear regression on cyclical components of unemployment and vacancy rates, the slope is very close to -1

In [17]:
x = hp_cycle_uv['Unemployment rate']
y = hp_cycle_uv['Vacancy rate']

slope, intercept = np.polyfit(x, y, 1)

print(slope)
print(intercept)

plt.scatter(x, y, facecolor='none', edgecolor='C0')
plt.plot(x, intercept + slope*x, 'r-', lw=2)

plt.xlim(-45, 45)
plt.ylim(-45, 45)

plt.hlines(0, -45, 45, linewidth=0.5)
plt.vlines(0, -45, 45, linewidth=0.5)

plt.title('Deviations from Hodrick-Prescott trend (%)')
plt.xlabel('Unemployment rate')
plt.ylabel('Vacancy rate')

# plt.savefig('BC_HP.pdf')

plt.show()
-0.9751202090453988
-3.785359803459378e-11

Generate the 'estimated' Beveridge curve without structural shifts

In [18]:
u = np.mean(dfq['Unemployment rate'])
v = np.mean(dfq['Vacancy rate'])

print(u, v)

scale = np.linspace(-40, 60, 101)

plt.plot(u*np.exp(scale/100), v*np.exp(intercept+slope*scale/100), 'r', lw=2)

plt.plot(u, v, 'ko')

plt.xlim(2, 12)
plt.ylim(1, 5)
plt.yticks(np.arange(1, 6))

plt.xlabel('Unemployment rate (%)')
plt.ylabel('Vacancy rate (%)')
plt.title('US Beveridge curve without structural shifts')

# plt.savefig('BC_est.pdf')

plt.show()
5.777576017463108 3.049301080368294
In [19]:
u = np.mean(dfq['Unemployment rate'])
v = np.mean(dfq['Vacancy rate'])

print(u, v)

scale = np.linspace(-40, 60, 101)

plt.plot(u*np.exp(scale/100), v*np.exp(intercept+slope*scale/100), 'r', lw=2)
plt.scatter(u*np.exp(x/100), v*np.exp(y/100), marker='.')

plt.xlim(2, 12)
plt.ylim(1, 5)
plt.yticks(np.arange(1, 6))

plt.hlines(v, 2, 12, linewidth=0.5)
plt.vlines(u, 1, 5, linewidth=0.5)

plt.xlabel('Unemployment rate (%)')
plt.ylabel('Vacancy rate (%)')
plt.title('US Beveridge curve without structural shifts')

# plt.savefig('BC_est_2.pdf')

plt.show()
5.777576017463108 3.049301080368294