# 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