A Theory of the Term Structure of Interest Rates under Limited Household Risk Sharing

Abstract

We present a theory in which the interaction between limited sharing of idiosyncratic labor income risk and labor adjustment costs (that endogenously arise through search frictions) determines interest rate dynamics. In the general equilibrium, the interaction of these two ingredients relates bond risk premiums, cross-sectional skewness of income growth, and labor market tightness. Our model rationalizes an upward-sloping average yield curve and predicts a negative relation between labor market tightness and bond risk premiums. We provide evidence for our theory’s mechanism and predictions.

Two patterns of interest rates of default-free bonds are well-known to be challenging for general equilibrium models. The first is an upward-sloping average yield curve (see Figure 1). The second is that bond excess returns are countercyclical (see our literature review for existing evidence). A common rationalization of these patterns builds on the Campbell and Cochrane (1999) model, which assumes a representative agent whose utility depends on an exogenously specified process called habit. Wachter (2006) shows that by appropriately choosing the habit process, it is possible to explain the two patterns.

In contrast to this top-down approach, in this paper, we offer a bottom-up approach in which we construct a production-based model to relate bond risk premiums to firms’ labor policies. Our model quantifies the interaction between two key ingredients: limited sharing of idiosyncratic labor income risk and labor adjustment costs that arise endogenously from search frictions. Our model rationalizes the two interest rate patterns mentioned above. In generating bond excess return predictability, our search-based model makes a new prediction: it predicts labor market tightness and the job finding rate, key variables in models of labor search (see, e.g., Shimer 2010), to be negatively related to bond risk premiums. We provide evidence for our model’s mechanism and predictions.

The two patterns of interest rates impose restrictions for the covariance of the stochastic discount factor (SDF) over the holding period of the bond and its remaining tenor. An upward-sloping average yield curve implies that this covariance is negative on average, that is, the SDF is mean reverting. Countercyclical bond excess returns implies that the covariance is more negative under bad aggregate conditions relative to good conditions. The canonical consumption-based model with perfect risk-sharing faces a challenge in satisfying these restrictions because they translate into restrictions for the growth rate of aggregate consumption which are falsified in the data. These challenges were pointed out by Backus, Gregory, and Zin (1989) in an endowment economy model; they also have been documented in workhorse macroeconomic models featuring production (see our literature review for examples). Our approach addresses these challenges.

Our production-based model features a representative firm which hires labor supplied by a cross-section of individuals. Our model’s first ingredient, limited risk sharing, implies an equilibrium SDF that is a consumption-weighted average of individual marginal utilities. The dynamics of this SDF is largely driven by the marginal utility process of a small fraction of individuals who experience large increases in marginal utility from income loss. The severity of income loss borne by this small fraction of individuals is measured by the cross-sectional skewness of labor income growth, henceforth “income skewness.” Our model’s second ingredient, labor adjustment costs, relates the speed of mean reversion in income skewness to the firm’s hiring policy. The latter depends on aggregate labor market conditions. In the general equilibrium, the interaction of the two ingredients links three quantities: the SDF (and hence bond risk premiums), income skewness, and aggregate labor market conditions. This link allows us to test our model’s mechanism and generates our model’s predictions.

Our model rationalizes an upward-sloping average yield curve. Figure 1 shows the slope of the nominal and real yield curves measured as the difference between the 10- and 1-year yields of U.S. Treasury bonds. From the figure, we see that both the nominal and real slopes are positive on average (see also Pflueger and Viceira (2016) for evidence of a positive average real slope).¹ In our model, an upward-sloping average yield curve results from mean reversion in income skewness. The latter depends on labor adjustment costs which we endogenize through labor search frictions. The labor search framework allows us to map a difficult-to-observe adjustment cost process to the observed processes for aggregate labor market variables, such as labor market tightness (ie, the ratio of job vacancies to unemployment) and the job finding rate. This allows us to discipline the speed of mean reversion of income skewness and hence the SDF. Similar to the habit-based model, the mean reversion in our SDF does not come through mean reversion in aggregate consumption; to leading order, mean reversion in our model’s SDF is determined by mean reversion in income skewness. Our model, therefore, avoids the challenges highlighted by Backus, Gregory, and Zin (1989).

We directly test for mean reversion in income skewness in the data since it is a critical component of our model’s mechanism for generating an upward-sloping average yield curve. We use the income skewness series from Guvenen, Ozkan, and Song (2014) who measure income skewness using administrative data from the U.S. Social Security Administration. We verify that income skewness does indeed mean revert.

Our model rationalizes countercyclical bond risk premiums. Specifically, our search-based model predicts bond risk premiums to be negatively related to aggregate labor market conditions as measured by labor market tightness and the job finding rate. The intuition is as follows. Firms reduce hiring during downturns which leads to an increase in the income risk of a larger than average fraction of individuals. This increase makes the SDF more volatile thereby increasing risk premiums of long-term bonds. Since such periods are also associated with a decline in labor market tightness and a decline in the job finding rate, we obtain a negative relation between labor market conditions and bond risk premiums. We test this prediction for both U.S. nominal bonds and Treasury Inflation Protected Securities and find support. For example, we find that a one-standard-deviation decrease in labor market tightness predicts a 1.05% increase in the excess return of an equal weighted portfolio of 2- through 5-year nominal bonds. While our search-based predictor variables are new, the existing literature includes evidence of bond excess return predictability by principal components of macroeconomic variables which load on labor market variables (see, e.g., Ludvigson and Ng 2009; Joslin, Priebsch, and Singleton 2014; Huang and Shi 2016, 2023; Bianchi, Büchner, and Tamoni 2020). Our theory therefore provides an explanation for these previous findings.²

We also show that yield-based variables predict bond excess returns in our model. Specifically, we run the Fama and Bliss (1987) predictability regressions in our model and show that, as in the data, the forward-spot spread predicts bond excess returns. Here we find that it is crucial to account for the slow recovery of search-based variables, such as tightness and the unemployment rate, in the economic recovery following a recession. While a recent literature (see, e.g., Hall and Kudlyak 2022) has documented this robust empirical pattern and analyzed its macroeconomic implications, our model shows its implications for asset prices.

In quantitatively evaluating our channel we minimize using asset pricing moments as calibration targets and instead use moments of real variables whenever possible.³ This strategy allows us to use labor market data to discipline the contribution of our channel in explaining term structure dynamics. For example, after matching the dynamics of cross-sectional income growth in the data, we find the model-implied loadings of bond excess returns on labor market tightness to be more than 70% of the data counterparts.

We show that the interaction of limited risk sharing of idiosyncratic labor income risk and labor adjustment costs can jointly rationalize key patterns of interest rates that are challenging for equilibrium models to explain. Our paper contributes to three strands of the literature and relates them.

First, we contribute to the literature that explains the dynamics of default-free interest rates. Evidence on the predictability of nominal bond excess returns has been documented by Fama and Bliss (1987), Campbell and Shiller (1991), and Cochrane and Piazzesi (2005), while evidence on the predictability of real bond excess returns has been documented by Pflueger and Viceira (2011, 2016). Backus, Gregory, and Zin (1989) highlight the difficulty for a consumption-based model with perfect risk sharing and power utility preferences to capture an upward-sloping average yield curve and predictability of bond excess returns.⁴ These challenges also have been documented in workhorse macroeconomic models featuring production (e.g., Donaldson, Johnsen, and Mehra 1990; den Haan 1995; Rudebusch and Swanson 2008; van Binsbergen et al. 2012). To the best of our knowledge, all consumption- and production-based models that successfully overcome these challenges assume perfect risk sharing among investors (with one exception that we reference below), but use richer specifications for investors’ preferences.⁵ Examples include habit formation (e.g., Wachter 2006; Chen 2017; Hsu, Li, and Palomino 2020), recursive preferences (e.g., Gallmeyer et al. 2007; Piazzesi and Schneider 2007; Le and Singleton 2010; Rudebusch and Swanson 2012; Bansal and Shaliastovich 2013; Kung 2015), and heterogeneity in investors’ preferences (Schneider 2022). In contrast to these models with perfect risk sharing, we rationalize an upward-sloping average yield curve and countercyclical bond risk premiums in a model in which nondiversifiable labor income risk plays a key role. A related paper is Kogan, Papanikolaou, and Stoffman (2020), who investigate the implications of a different form of imperfect risk sharing (the inability of investors to share displacement risks associated with future technological innovations) for stock returns. Their model also generates an upward-sloping average yield curve. However, it is unclear whether their model can explain the predictability of bond excess returns over the business cycle (they do not report results regarding the cyclical properties of the yield curve).

Second, we contribute to the literature that analyzes the asset-pricing implications of nondiversifiable idiosyncratic labor income risk. The idea that nondiversifiable labor income risk can have a first-order effect on equities goes back to at least Mankiw (1986) and Constantinides and Duffie (1996).⁶ More recently, Constantinides and Ghosh (2017) and Schmidt (2022) highlight the importance of disasters at the individual level for stock returns. While labor income risk is exogenously specified in the prior literature, to the best of our knowledge, we are the first to consider the asset pricing implications of nondiversifiable labor income risk where this risk is derived from firms’ labor market policies. Therefore, our theory provides an explanation for the findings of Ludvigson and Ng (2009), Joslin, Priebsch, and Singleton (2014), Huang and Shi (2016, 2023), and Bianchi, Büchner, and Tamoni (2020) that bond excess returns are predicted by principal components of macroeconomic variables that load on labor market variables (including employment, unemployment, and vacancies, amongst others). We focus on the implications of limited risk sharing for default-free bonds because, in contrast to other asset classes (e.g., equities), they do not have cash flow risk and therefore offer a cleaner test of our economic mechanism.

Third, we contribute to the literature that analyzes the importance of labor adjustment costs for asset prices. This importance has been highlighted by Belo, Lin, and Bazdresch (2014) for the cross-section of stock returns. While the existing literature mostly uses a reduced form specification for adjustment costs, these costs arise endogenously in our model from search frictions. The advantage of the labor search framework is that it maps the observed processes for aggregate labor market variables into an otherwise difficult-to-observe adjustment cost process. This disciplines the adjustment cost process which, in turn, places restrictions on the SDF and hence interest rates. The importance of search frictions also has been explored by Petrosky-Nadeau, Zhang, and Kuehn (2018) for the aggregate stock market, and by Kuehn, Simutin, and Wang (2017) for the cross-section of stock returns. More generally, labor market frictions has been shown to be important in jointly accounting for asset pricing and macroeconomic facts. For example, the importance of accounting for wage rigidity has been highlighted by Uhlig (2007) and Favilukis and Lin (2016) for stock returns, and by Favilukis, Lin, and Zhao (2020) for defaultable corporate bonds.

1 Model

This section presents our general equilibrium model of interest rates. Section 1.1 describes the economy, Section 1.2 derives the equilibrium, and Section 1.3 solves for the term structure of interest rates.

1.1 The economy

The economy is set in discrete time, with the horizon being infinite. There is a single aggregate productivity shock whose value $zt$ evolves according to a first-order Markov chain.

1.1.1 The household

The “head of household” is a modeling device that captures the essence of imperfect risk sharing in a heterogenous agent production-based setting, while still preserving the tractability of the representative agent framework. We further discuss this preference assumption in Section 1.2.7.

1.1.2 Limited risk sharing

To make the model tractable, we assume that, after conditioning on $Nt$ and $Ut$⁠, the idiosyncratic employment shocks $eit$ are independent both across individuals and over time. As a result, the employment and unemployment probabilities, $Nt$ and $Ut$⁠, also correspond to the equilibrium aggregate employment and unemployment rates, respectively, whose dynamics are described in Section 1.1.3.

We assume that the head of household can trade a complete menu of state-contingent payoffs for aggregate risks. From the Fundamental Theorem of Asset Pricing (see, e.g., Dybvig and Ross 2003), the absence of arbitrage implies the existence of a SDF, $Mt,t+n$⁠, which prices returns between t and $t+n$⁠, $Rt,t+n$⁠, according to the asset pricing relationship $1=Et[Mt,t+nRt,t+n].$ We characterize the equilibrium SDF in Section 1.2.5.

1.1.3 Labor market search frictions

Employed individuals are paid wages $wt$ determined using a generalized Nash-bargaining protocol, described in more detail in Section 1.2.3, in which employees obtain a fraction $η$ of the surplus. Unemployed individuals are paid the amount b in each period of unemployment; these unemployment benefits are funded by lump-sum taxes.

1.1.4 Timing of events

Figure 2 illustrates the timing of events within each period, which is as follows:

• At the start of period t, there is a mass $Nt∈[0,1]$ of employed individuals, and $Ut=1−Nt$ unemployed individuals. Nature draws aggregate productivity $zt$ according to its law of motion. Nature also draws the idiosyncratic employment shocks $eit$⁠.

• The head of household chooses its policies. This includes (1) the consumption of the employed, $Ce,t$⁠, and the unemployed, $Cu,t$⁠, and (2) portfolio choices.

• The representative firm posts vacancies $Vt$⁠, and labor market matching takes place. Matched individuals become employed at the start of the next period.

• Wages are set via a generalized Nash bargaining rule with employed individuals capturing a fraction $η$ of the surplus.

• Production takes place and output $Yt$ is realized. Wages are then paid, unemployment benefits are collected, and consumption takes place.

• Existing matches (excluding newly formed ones) exogenously separate with probability s.

1.2 Equilibrium

1.2.1 Firm’s problem

The dividend in period t is $Dt=Yt−wtNt−κVt$ where output $Yt$ is given by (6), $wt$ is the wage paid to each employed individual, and $Vt≥0$ is the number of vacancies posted by the firm. Future dividends are discounted using the SDF $Mt,t+k$⁠.

The upper bound on matches $m¯t$ given by Equation (8) limits the degree to which vacancies translate into new matches. Specifically, resources dedicated to vacancies $Vt$ beyond the threshold $V¯t$ do not create additional matches and are wasted, where the threshold $V¯t$ solves $g(V¯t/Ut)V¯t=m¯t$⁠. We assume that the firm recognizes the threshold $V¯t$ induced by the upper bound $m¯t$ and limits vacancy posting to $Vt≤V¯t$⁠.¹⁰

1.2.2 Household’s problem

Finally, lump-sum taxes are levied to exactly finance unemployment benefits. That is, $Tt=bUt$ for all t.

1.2.3 Wages

The derivation of this result is in Section B.1 of our Internet Appendix.

1.2.4 Equilibrium

In equilibrium, all policies and value functions are a function of state variables $zt$⁠, $Φt$⁠, and $Nt$⁠. For example, $Θt=Θ(zt,Φt,Nt)$⁠, $wt=w(zt,Φt,Nt)$⁠, and so on.

1.2.5 Equilibrium SDF

The first component of the SDF (28), $βn(Ct+n/Ct)−γ$⁠, depends only on aggregate consumption growth and is the only term that would appear under perfect risk sharing (see, e.g., Breeden 1979). This property is illustrated by the dashed line in Figure 3: when $Φt=1$ for all t, $ζt=1$ and hence $ζt,t+n=1$ for all t so that idiosyncratic labor income risk does not affect the SDF.

The second component of the SDF (28), $ζt,t+n$⁠, arises as a result of limited risk sharing (ie, $Φt<1$ in some states). In our model, the dispersion in nondiversifiable employment outcomes across individuals becomes larger at lower levels of employment $Nt$ and risk sharing capacity $Φt$⁠, which raises the marginal utility of the head of household through the $ζt$ term (27). These properties are illustrated by Figure 3 which shows that $ζt$ is decreasing in employment $Nt$ (for $Φt<1$⁠) and risk sharing capacity $Φt$ (for $Nt<1$⁠). As a result, $log ζt,t+n=log ζt+n−log ζt$ depends on changes in $Nt$ and $Φt$ between t and $t+n$⁠.

1.2.6 Labor adjustment costs

We use the labor search framework to use search-based variables to discipline labor adjustment costs. The latter is an important determinant of bond risk premiums in our setting—bond risk premiums depend on the $ζt,t+n$ component of the SDF (we provide the details in Section 1.3) whose behavior, in turn, depends on labor adjustment costs.

Equation (29) illustrates how the search framework allows us discipline an otherwise difficult-to-observe adjustment cost process. Specifically, the observed processes for unemployment rates, labor market tightness, and job flow rates allow us to discipline the curvature of the matching function $ι$⁠, vacancy posting costs $κ$⁠, and the equilibrium benefit of hiring an additional worker $Et[Mt,t+1∂Ft+1/∂Nt+1]$⁠. Similarly, when Equation (15c) holds with equality so that the economy is in a slow labor market recovery, the observed rate of labor market recovery allows us to discipline the maximal fractional unemployment decline $1−a$ following recessions.

1.2.7 Discussion of assumptions

We assume (A1) a head of the household with preferences (1) who allocates consumption subject to the constraint (5) and (A2) iid employment shocks (3) after conditioning on aggregate employment. Assumptions (A1) and (A2) make our model as tractable as a representative agent model while capturing the effect of limited sharing of idiosyncratic labor income risk on bond prices.

Taken together, assumptions (A1) and (A2) allow us to avoid having to keep track of the cross-sectional distribution of wealth which determines consumption and portfolio choices in the cross-section. This is especially convenient in our production based model because, unlike an endowment economy where one is able to specify income processes that make the model tractable (e.g., an income process that ensures no-trade among investors as in Constantinides and Duffie 1996), the income process of individual agents in our model is endogenously determined and varies with macroeconomic conditions.

Our assumption that the head’s preference is an aggregate of individual utilities (1) implies that nondiversifiable labor income risk affects the equilibrium SDF. For example, the increase in this SDF in regimes with lower $Φt$⁠, that is, with poorer risk sharing, is largely driven by the increase in marginal utility of a small fraction of individuals who experience income loss, since the head allocates lower consumption to these individuals according to Equation (5). This drives up the consumption-weighted average marginal utility, and hence $ζt$ through Equation (26).

1.3 Term structure of interest rates

The main result of this section is a decomposition of bond risk premiums into a standard component (ie, bond risk premiums under perfect risk sharing) and components that arise only due to limited risk sharing. This allows us to gauge the contribution of limited risk sharing to interest rates.

1.3.1 Real term structure

The time t price of a default-free zero coupon bond which matures in n periods time and pays off a unit in real terms at maturity is given by $Pt(n)=Et[Mt,t+n]$⁠. The corresponding yield to maturity is $yt(n)≡−1nlog Pt(n)$⁠.

1.3.2 Decomposition of real bond risk premiums

The first approximation (37a) shows that $hpxrtH,T,C$ depends only on the dynamics of aggregate consumption growth. Therefore, we refer to this term as the “Aggregate consumption” term. Under perfect risk sharing (ie, $Φt=1$ for all t), only the “Aggregate consumption” term, $hpxrtH,T,C$⁠, shows up in Equation (33). This leads to the bond risk premium puzzle highlighted by Backus, Gregory, and Zin (1989). They show that matching the observed (positive) bond risk premium would then require aggregate consumption growth to have a counterfactually large and negative autocorrelation.

Limited risk sharing (ie, $Φt<1$ for some t) leads to two new terms, $hpxrtH,T,ζ$ and $hpxrtH,T,cross$⁠, both of which depends on the time-series properties of $ζt,t+n$⁠. The second approximation (37b) shows that $hpxrtH,T,ζ$ depends on the cross-sectional dispersion in consumption growth. This is because $log ζt,t+n=log ζt+n−log ζt$ is the growth in $ζt$⁠, and $ζt$ depends on the cross-sectional dispersion in consumption levels through Equation (26). The premise of our paper is that the source of this dispersion in consumption growth is the inability of investors to diversify their idiosyncratic income risk. Therefore, we call this term “Income dispersion.” The third approximation (37c) shows that $hpxrtH,T,cross$ depends on the cross-covariance between aggregate consumption growth and $log ζt,t+n$⁠, and accordingly, we label it as the “Cross-covariance” term.

1.3.3 Nominal term structure

We introduce inflation to obtain nominal bond prices. We do so because using nominal bond prices to test our model’s predictions offers three practical advantages. First, U.S. nominal bond price data has a substantially longer history compared to U.S. TIPS data for real bonds. Second, to extract the real yield curve from TIPS data, it is necessary to account for a liquidity premium (see, e.g., Pflueger and Viceira 2016). Third, using nominal bond prices allows us to avoid the so-called “indexation lag” problem which makes the prices of short maturity TIPS erratic (Gurkaynak, Sack, and Wright 2010, p. 76).¹² Inflation otherwise plays a secondary role in our model.

We provide details for the computation of the nominal bond price (38) in  Appendix B.3.

2 Quantitative Analysis

In this section, we study our model’s quantitative implications. Section 2.1 describes the calibration and Section 2.2 reports the term structure of interest rates. The next section, Section 3, investigates the link between labor market conditions and bond risk premiums.

2.1 Calibration

In calibrating our model, we minimize using asset pricing moments as targets and instead use moments of labor market variables whenever possible. Although this calibration strategy ties our hands in better matching asset pricing moments, it allows us to better discipline the contribution of our channel in explaining asset prices, especially properties of the term structure. The only asset pricing moments that we use are the first two moments of the nominal 1-year yield to determine two preference parameters—the time preference parameter $β$ and the elasticity of substitution $χ$⁠.

We solve our model numerically using global methods. Section B.2 of the Internet Appendix provides details for the numerical algorithm and illustrates firm policies and aggregate variables as a function of the state variables. We simulate our model at monthly frequency using the parameters shown in Table 1. We report our model-implied moments over a quarterly frequency in Table 2, together with their data counterparts. We first discuss the calibration of the risk sharing process in Sections 2.1.1 and 2.1.2. We then discuss the calibration of the remaining parameters in Section 2.1.3. Details for the data used in our calibration procedure are available in  Appendix A.

2.1.1 Calibration of the $ζt$ process

Our model’s SDF depends on cross-sectional differences in consumption risk through the $ζt$ term (27) whose dynamics depend on employment $Nt$ and the risk sharing process $Φt$⁠. The behavior of $Nt=1−Ut$ depends on labor search frictions which we pin down via properties of the unemployment rate $Ut$ in the data.

We calibrate the $Φt$ process so that our model generates realistic cross-sectional differences in consumption risk. Specifically, we choose the $Φt$ process so that our model-implied third central moment of individuals’ consumption growth, $μ3,t,t+n$⁠, is in line with that of the data. In measuring $μ3,t,t+n$ in the data, we use measures of cross-sectional differences in income risk from Guvenen, Ozkan, and Song (2014) and prior estimates for the pass-through from income shocks to consumption. We explain the details for this calibration strategy for $Φt$ in the remainder of this section.

Measuring cross-sectional consumption risk

Here, $Δct+1=ct+1−ct$ is the change in aggregate log consumption $ct≡∫01citdi$⁠, and the consumption shock $ϵi,t+1=Δci,t+1−Δct+1$ determines individual i’s consumption growth relative to that of the aggregate. Equation (41) is an accounting identity¹³ and the consumption shocks $ϵi,t+1$ embed both unemployment risk and the earnings risk of employed workers.

Equations (44a) and (44b) are our measurement equations for $σt,t+T$ and $μ3,t,t+T$⁠, respectively. In our baseline calibration, we choose $α=0.8$ which is the recent estimate for income pass-through from Agarwal and Qian (2014). This choice for the value of $α$ is also in the middle of the range of estimates obtained by Blundell, Pistaferri, and Preston (2008). We use data from Guvenen, Ozkan, and Song (2014), henceforth GOS, to construct $σt,t+Tinc$ and $μ3,t,t+Tinc$⁠. The GOS series for $σt,t+Tinc$ and $μ3,t,t+Tinc$ are available for horizons of $T=12$ months and $T=60$ months and are computed using data from a large cross-section (a 10% sample of all U.S. working-age males from 1978 to 2011 taken from the U.S. Social Security Administration). In Section 2.1.2, we discuss the sensitivity our results to alternative estimates for income pass-through and cross-sectional income risk.

We measure an approximately constant $σt,t+T$ and a procyclical $μ3,t,t+T$ (⁠$μ3,t,t+T$ is negative on average and becomes more negative during recessions); this follows from GOS’s findings for an acyclical $σt,t+Tinc$ and a procyclical $μ3,t,t+Tinc$⁠. These findings imply that $log ζt,t+T$ is predominantly driven by changes in $μ3,t,t+T$ rather than $σt,t+T$ (see equation (42)). For this reason, we calibrate the $Φt$ process based on $μ3,t,t+T$⁠.

Specification of $Φt$

We choose $x¯$ based on consumption decline following unemployment estimated in Ganong and Noel (2019) (henceforth “GN”), who find “… a drop of 6 percent at unemployment onset, a drop of less than 1 percent per month during UI receipt, and a 12 percent drop at UI benefit exhaustion” (p. 2384). GN’s estimates, which are based on the expansionary period 2014—2016, together with an unemployment duration of 3 months,¹⁵ implies a target of $E[Φt|zt=z2]=1−(0.06+3×0.01)=0.91$ in the high productivity regime. We set $x¯=2.557$ based on this target.

We choose $σx=0.15$ to match $σ(μ3,t,t+12)$ in the data. To obtain the data value, we first note that Equation (44b) implies $σ(μ3,t,t+12)=α3σ(μ3,t,t+12inc)$⁠. Next, we use the GOS series to construct the time-series $μ3,t,t+12inc$ from which we estimate $σ(μ3,t,t+12inc)=0.0304$⁠. Using our baseline value of $α=0.8$⁠, we find $σ(μ3,t,t+12)=α3σ(μ3,t,t+12inc)=0.0156$⁠; the corresponding model-implied value is also 0.0156, where we use Equation (45) to compute $μ3,t,t+12$ and its volatility $σ(μ3,t,t+12)$ in the model.

We choose $ρx=0.9112=0.9913$ to approximately match the mean reversion of $μ3,t,t+T$ in the data. To do so, first note that Equation (44b) implies that the autocorrelation coefficients $Corr(μ3,t,t+12,μ3,t+12,t+T)=Corr(μ3,t,t+12inc,μ3,t+12,t+Tinc)$ for all T. We use the GOS series to estimate $Corr(μ3,t,t+12inc,μ3,t+12,t+Tinc)$ for various T; the values are reported in row (2) of Table 7. The mean autocorrelation averaged over $T∈{24,36,48,60,72}$ months is $−0.15$ in the data; the corresponding model-implied value is $−0.22$⁠, where we use Equation (45) to compute $μ3,t,t+12$ and $Corr(μ3,t,t+12,μ3,t+12,t+T)$ in the model.

2.1.2 Discussion: Measuring cross-sectional consumption risk

Use of income data

In using the measurement Equation (44), we do not directly use cross-sectional consumption data to estimate $σt,t+T$ and $μ3,t,t+T$⁠. This is because consumption data does not allow us to directly estimate individual consumption risk for horizons longer than a year,¹⁶ which is needed for estimating the contribution of limited risk sharing to the risk premiums of long-term bonds (see Section 4).

The use of GOS income data in implementing Equation (44) offers some additional advantages. First, the administrative nature of the GOS data set avoids issues with measurement errors associated with survey-based data of individual consumption (see the discussions in Koijen, Nieuwerburgh, and Vestman [2014, p. 309] and Guvenen, Ozkan, and Song [2014, p. 622]). Second, the large sample size in GOS helps avoid the robustness concerns associated with estimating higher moments using small samples (see, e.g., Kim and White 2004).

A potential concern with using GOS data, however, is that it may overstate the amount of income risk faced by households. This is because GOS data measures income risk faced by working-age males which ignores within-household smoothing of individual income shocks. In Section B.3.1 of the Internet Appendix, we show that our model’s implications for bond risk premiums are robust to adjusting for within-household smoothing after recalibrating preference parameters according to the procedure in Section 2.1.3.

Finally, we note that it is reassuring that direct estimates of $μ3,t,t+T$ using consumption data results in similar properties for $μ3,t,t+T$ compared to our income-based estimate. For example, Constantinides and Ghosh (2017) use the Consumption Expenditure Survey to estimate $μ3,t,t+T$ for a horizon of $T=3$ months. They find $μ3,t,t+T$ to be procyclical and negative on average.

Alternative values of income pass-through

The literature that estimates the pass-through of income shocks to consumption $α$ arrives at a range of estimates for this parameter. In addition to our baseline choice of the value of $α$⁠, we examine the robustness of our main results as we vary $α$ over the range $α=0.35$ and $α=0.94$ that have been estimated in the literature. At the low end, $α=0.35$ corresponds to the midpoint of estimates from Arellano, Blundell, and Bonhomme (2017, p. 717), who find “[o]n average, the estimated [passthrough] parameter lies between 0.3 and 0.4.” At the high end, $α=0.94$ corresponds to the largest estimate for $α$ in Blundell, Pistaferri, and Preston (2008, table 6, row for “Partial insurance perm. shock”).

We perform two exercises to investigate the effect of $α$ on bond risk premiums. In the first exercise, we vary $α$ while calibrating all other model parameters according to our calibration strategy. In particular, for each choice of $α$⁠, we also simultaneously calibrate the elasticity of substitution across individuals’ consumption $χ$ to match the volatility of the 1-year yield, the latter being our calibration strategy to choose $χ$ (see Section 2.1.3). We find our results to be robust and report them in Section B.3.1 of the Internet Appendix. While it might seem surprising that the term structure of yields and yield volatilities do not change much as $α$ is varied over the wide range of $α=0.35$ to $α=0.94$⁠, the results are robust because lower values of $α$ are associated with lower values of the fitted preference parameter $χ$ (from targeting the volatility of the 1-year nominal rate). The latter prevents bond risk premiums and the slope of the term structure from decreasing.¹⁷

The second exercise is a sensitivity analysis: we vary $α$ over the same range 0.35 to 0.94, but do not additionally recalibrate $χ$ to target the volatility of the 1-year nominal yield. We find that the model implied volatility of yields (both nominal and real), the slope of the term structure of yields (both nominal and real), the magnitude of the slope coefficients of the bond predictability regressions using labor market variables, and the Fama-Bliss regression slope coefficients all increase when $α$ increases. This effect works through imperfect risk sharing—the increase is purely through the larger income dispersion (35) and cross-covariance (36) terms for larger values of $α$⁠. We report these results in Section B.3.2 of the Internet Appendix.

2.1.3 Calibration of remaining parameters

We model the exogenous productivity process $zt$ as a two-regime Markov chain for our baseline analysis. Having a simple two-regime specification for productivity makes it easier to explain the intuition of our model, and all results presented in our main text are for this two-regime specification. We show in  Appendix C that our results are robust to increasing the number of states for $zt$⁠. Specifically, we present results for the case in which aggregate productivity follows an AR(1) process; the main findings remain unchanged.

We follow Barton, David, and Fix (1962) and, without loss of generality, parameterize the transition probabilities for the two regimes to be $pij=(1−λ)pj+λ$ if $i=j$⁠, and $pij=(1−λ)pj$ if $i≠j$⁠, where $pj$ is the long-run probability of regime j, and the parameter $λ$ captures the persistence of aggregate regimes (⁠$λ=0$ corresponds to i.i.d regimes). We choose $p1=1/6$⁠, $p2=5/6$⁠, and $λ=0.9$⁠. These choices imply an average duration of 1 and 5 years for the low and high productivity regimes, respectively, which are in line with the durations of NBER recessions and expansions over the period 1964Q1-2016Q4. Our sample stops at the end of 2016 because the updated Barnichon (2010) Help-Wanted-Index, which we use to construct labor market tightness, stops then. We normalize log productivity in the expansionary regime to $log z2=0$⁠. We then set the difference in log-productivity between the two regimes to be $log z2−log z1=0.0355$ to match the unconditional volatility of U.S. gross domestic product (GDP) per capita growth over the period 1964Q1-2016Q4 which is $0.81%$⁠.

There are six labor market parameters. We estimate the monthly job separation probability s and the curvature of the matching function $ι$ directly from the data over the period 1964Q1—2016Q4. We set $s=3.4%$ to match the average monthly job separation rate (constructed following Elsby, Michaels, and Solon 2009). We obtain a curvature of $ι=1.24$ for the matching function (7) by minimizing the sum square error of the difference between the empirically observed job finding rate (constructed following Elsby, Michaels, and Solon 2009) and its model counterpart. The latter is obtained by feeding the observed times-series for labor market tightness into the model’s specification for the job finding probability (10). We choose the vacancy posting cost $κ=0.1067$ and the unemployment benefit parameter $b=0.9362$ to match the mean unemployment rate of 6.09% in the United States over the period 1964Q1—2016Q4 and the volatility of $0.78%$ of the (HP filtered) series over this period. Our choice of the value for b is within the range of estimates of b used in the literature which ranges between 0.4 used by Shimer (2005) to 0.955 used by Hagedorn and Manovskii (2008). We choose the bargaining power of employed individuals to be $η=0.312$ to match the ratio of the volatility of the total wage bill to the volatility of output of 0.87 both in the data (see, e.g., Favilukis and Lin 2016, Table 3) and in our model. Finally, we set $a=exp (−0.1/12)=0.9917$ so that the annualized rate of decline of the unemployment rate following recessions is 10% of the prior level which is in line with the finding in Hall and Kudlyak (2022).

Although not target moments, our model implied moments for labor market-tightness and its correlation with the unemployment rate are close to their data counterparts. Mean labor market tightness is 0.58 in the data and 0.90 in our model. The volatility of the HP filtered series for labor market tightness is 0.14 in the data and 0.19 in our model. The correlation of the unemployment rate and vacancies is -0.90 in the data and -0.91 in our model.

As discussed in Section 1.3.3, our yield-curve-based calibration targets and predictions focus on moments of the nominal yield curve instead of the real yield curve. Our model’s predictions for the nominal term structure are based on the monthly ARMA(1,1) inflation process (39). We obtain parameter estimates for this process using maximum likelihood estimation over the monthly sample 1964m1—2016m12. As inputs into our estimation procedure, we measure inflation using the monthly Consumption Price Index, and obtain the data counterpart to the consumption growth shocks (40) using residuals from an AR(1) model fitted to log aggregate consumption growth data. The resultant parameter estimates are: an average monthly inflation of $μπ=0.00325$⁠, an autocorrelation coefficient of $ρπ=0.81$⁠, a loading on the innovation to aggregate consumption of $ξπ=−0.035$⁠, a moving-average coefficient of $νπ=−0.338$⁠, and a volatility of $σπ=0.00245$ for the residuals $ϵπ,t$⁠.

We choose the head of household’s coefficient of relative risk aversion $γ=2$ to approximately match the volatility and autocorrelation of aggregate consumption growth.¹⁸ The model-implied autocorrelation and volatility of quarterly consumption growth are 0.32 and $0.79%$⁠, respectively; these values are close to their data counterparts which are 0.32 and $0.67%$⁠, respectively.¹⁹

The remaining two preference parameters which include the head of household’s time preference parameter $β$ and the elasticity of substitution across individuals’ consumption $χ$ are difficult to estimate using real variables. We therefore follow the standard approach in equilibrium asset pricing²⁰ and choose these two preference parameters by targeting asset pricing moments. Specifically, we choose the head of household’s time preference parameter $β=0.9982$ and the elasticity of substitution across individuals’ consumption $χ=0.2611$⁠, to match the mean and volatility of the 1-year nominal yield of 5.28% and 3.32%, respectively.

The equity premium is not a targeted moment of our term structure model. Table 2 shows that our model generates a reasonable annual equity premium of 4.82% compared to $6.33%$ in the data. In computing this premium, we define equity as a levered claim on the unlevered value of the firm with a leverage-factor of three (Abel 1999).²¹ Our model-implied volatility of excess returns of the stock market is $34.14%$⁠, which is higher than the volatility of $17.62%$ in the data. This indicates that our general equilibrium model, in which the market price of risk is determined by real variables, generates smaller Sharpe ratios than in the data. Our model’s prediction for the (annualized) dividend-to-price (DP) ratio is in line with the data. The mean of the DP ratio is 0.029 and 0.022 in the data and model, respectively. The volatility of the DP ratio is 0.019 in the data and 0.017 in the model.

2.2 Term structure of interest rates

In this section, we explain how our model with imperfect risk sharing generates a yield curve that is upward sloping on average (both nominal and real) without counterfactual implications for the dynamics of aggregate consumption growth. We show that imperfect risk sharing is predominantly responsible for this positive slope. We conclude this section by showing that our model predicts countercyclical bond risk premiums.

2.2.1 Moments of yields

We report average yields for maturities of 1–5 years, as observed in the data and implied by our model, in panel A of Table 3. The data numbers for nominal yields in row 1 are computed from the monthly zero coupon bond yields from the Fama-Bliss data set over the period 1964m1—2016m12. We see in row 2 that our model captures the observed upward slope of the nominal term structure; the slope, as measured by the difference in yields between the 5- and 1-year bond, is 0.69% in the data and 0.32% in our model. Our model therefore explains about half of the observed slope of the nominal yield curve.

Row 3 of panel A reports average real yields as implied by our model. We do not report empirical estimates of average real yields in Table 3. There are several empirical estimates for moments of the real yield curve in the literature; we report them in Section A.1 of our Internet Appendix and summarize the key findings below. The values in row 3 of panel A show that our model implies an upward-sloping average real yield curve: the real slope, as measured by the 5- minus 1-year real yield spread, is 0.32% on average in the model. Our model’s implication of an upward-sloping average real yield curve is in line with empirical estimates from the existing literature. For example, the slope is 0.40% on average based on estimates of the real yield curve in Chernov and Mueller (2012) for the period 1971m1—2002m4. Similarly, the slope is 0.48% based estimates from Haubrich, Pennacchi, and Ritchken (2012) for the period 1982m1—2009m5.²² Our model therefore explains between two-thirds and four-fifth of the average slope of the real yield curve (depending on the estimate).

In panel B of Table 3, we report volatilities for yields for maturities between 1 and 5 years, as observed in the data and as implied by our model. Comparing rows 1 and 2, we see that while our model matches the volatility of the 1-year nominal yield (this was a calibration target), the model-implied volatility of nominal yields declines more rapidly with maturity than in the data. The average volatility of the nominal term structure is 3.17% in the data and 2% in our model.

Our model implied volatility of real yields is in line with the data. For instance, at the short end, the model-implied volatility of the 1-year real yield is 3.10% (see row 3 of panel B in Table 3). Estimates of the corresponding moment in the data range from 1.76% to 3.08% (see Section A.1 of the Internet Appendix). Our model also captures the decline in volatility of real yields with maturity. For the 5-year yield, our model-implied volatility is 1.26%, while in the data this value ranges from 1.07% to 1.61% (see Section A.1 of our Internet Appendix).

Comparing rows 5 and 6 of panel C, we see that our model-implied nominal bond risk premiums are predominantly driven by real bond risk premiums. Our model’s predicted bond risk premiums are in line with estimates from the literature. For example, our model-implied real risk premium for the 5-year bond is 52 basis points compared to the estimate of 54 basis points from Haubrich, Pennacchi, and Ritchken (2012, p. 1618). Our model captures the fact that the mean and volatility of excess returns are increasing in maturity. However, we see that excess returns are more volatile at shorter maturities in the model relative to the data, with the fit improving at higher maturities. The average of the volatilities of excess returns between maturities 2 through 5 years is 3.7% in the data and 4.7% in the model.

2.2.2 Contribution of limited risk sharing

We use the decomposition (33) to study the contribution of limited risk sharing to the positive slope implied by our model. The 5- minus 1-year unconditional real slope is related to the unconditional risk premium of an equal-weighted portfolio of 2- through 5-year real bonds through the identity (47) (with $H=1$ year and $n=5$⁠).

Column 1 of Table 4 reports the decomposition of the unconditional risk premium for the equal-weighted portfolio (over a holding period of a year). The mean bond risk premium is 0.416%, with 88% of this premium due to imperfect risk sharing (ie, the sum of the “Income dispersion” (35) and the “Cross-covariance” (36) terms). Aggregate consumption risk (ie, the “Aggregate consumption” term (34)) contributes 12% of the total risk premium. Similarly, columns 2 and 3 show that the majority of the conditional bond risk premium in each of the two productivity regimes is due to labor income dispersion. Because a significant portion of the risk premium arises due to nondiversifiable idiosyncratic risk, our model provides a potential resolution of the puzzle highlighted by Backus, Gregory, and Zin (1989). That is, our model predicts long-term bonds to have a positive and empirically realistic bond risk premium without predicting a counterfactually large negative autocorrelation of aggregate consumption growth.

The quantitative effect of imperfect risk sharing on bond risk premiums depends on the values of the preference parameter $χ$ and the pass-through parameter $α$⁠. This can be seen, for instance, from the approximation of the SDF in Equation (42), where $χ$ appears explicitly and $α$ implicitly through higher moments of cross-sectional consumption growth $σt,t+T$ and $μ3,t,t+T$ (see Equation (44)). Equation (42) shows that, all else equal, imperfect risk sharing becomes more important for higher values of $α$ or higher values of $χ−1$⁠. While we followed the standard approach in equilibrium asset pricing and calibrated the value of the preference parameter $χ$ to target the volatility of the 1-year nominal yield, it will be interesting for future research to use joint data on individuals’ portfolio holdings and income (or, more generally, characteristics that predict individuals’ income risk) to provide more direct micro-level estimates of parameters, such as $χ$⁠, to determine the strength with which imperfect sharing of income risk affects term structure dynamics, or asset prices in general.

In Section B.4 of the Internet Appendix, we report results of a sensitivity analysis that illustrates the role of the parameters for the risk sharing process (ie, $x¯$⁠, $ρx$ and $σx$⁠) and the preference parameter $χ$ for bond risk premiums. We find that lowering $x¯$ or $χ$ increases bond risk premiums by reducing risk sharing $Φt$ and by making $ζt,t+T$ more volatile (see, e.g., the approximation (42)), respectively. In turn, this leads to a more volatile SDF. Increasing $ρx$ or $σx$ increases bond risk premiums by increasing the unconditional volatility of $Φt$ which makes the SDF more volatile. We additionally find that changes in $x¯$⁠, $ρx$⁠, $σx$⁠, and $χ$ have similar effects for the average equity risk premium. Section B.4 of the Internet Appendix also present additional comparative static results that investigate how limited risk sharing and labor adjustment costs influence the term structure of interest rates.

2.2.3 Countercyclical bond risk premiums

Columns 2 and 3 of Table 4 show that bond risk premiums are countercyclical in our model. For example, the bond risk premium for the equal-weighted portfolio is 1.12% on average when productivity is in the low regime versus 0.28% when productivity is in the high regime.

Risk premiums increase in a downturn for two reasons. First, a decline in the firm’s hiring (during a downturn) leads to an increase in the income risk of a larger than average fraction of individuals. Second, risk sharing $Φt$ worsens in a recession (⁠$σx>0$ in Equation (46)). As a result, the marginal utility of the head of the household increases through the $ζt$ term in the SDF (28).

In addition, labor market variables are procyclical. For example, the mean of tightness is $E[Θt|zt=z1]=0.61$ and $E[Θt|zt=z2]=0.95$ in recessions and expansions, respectively. Similarly, the job finding rate averages $E[ft|zt=z1]=0.42$ and $E[ft|zt=z2]=0.55$ in recessions and expansions, respectively. This generates a natural prediction: labor market variables, such as tightness and the job finding rate, can be expected to negatively predict bond excess returns. We provide evidence for this prediction in Section 3.1.

3 Labor Market Conditions and Bond Risk Premiums

In this section we focus on time variation in bond risk premiums. In Section 2.2.3, we showed that our search-based model predicts bond risk premiums to be negatively related to labor market tightness and the job finding rate, both of which are key variables in search models (see, e.g., Shimer 2010). We provide evidence for this prediction in Section 3.1. To the best of our knowledge, this evidence is new. In Section 3.2, we show that our model is able to rationalize the Fama-Bliss predictability regressions. There, we show the importance of accounting for a robust stylized fact of U.S. labor markets (see, e.g., Hall and Kudlyak 2022)—slow recoveries following the end of recessions—in generating the Fama-Bliss predictability results in our model.

3.1 Tightness predicts bond excess returns

The left-hand-side variable is the log excess return of a $n∈{24,36,48,60}$ month nominal bond, over a holding period of 12 months. The predictive variable is either labor market tightness (⁠$Xt=Θt$⁠) or the job finding rate (⁠$Xt=ft$⁠). Our sample consists of monthly observations of annual returns for the period 1964m1—2016m12 (see  Appendix A for details regarding the data). The summary statistics for the right-hand side variables appearing in regression (48) are in Table A1; the summary statistics for bond excess returns are in panels C and D of Table 3. Afterwards, we compare our model-implied values for $β(n)$ to their data counterparts in order to assess the quantitative performance of our model.

Row 1 of panel A of Table 5 shows that, in line with our model’s prediction, labor market tightness is negatively correlated with future bond excess returns in the data. We report Newey-West standard errors with 12 lags to account for the overlapping windows in the regressions. The estimated coefficients imply that a one-standard-deviation decrease in labor market tightness of 0.26,²⁴ see Table A1, is associated with an increase of 0.47% for the 2-year bond to 1.65% for the 5-year bond; the associated $R2$ ranges between 7.3% for the two year bond and 9.3% for the five year bond. The loadings are statistically significant (see row 2 of same panel). As a summary measure of the average response of bond excess returns across maturities ranging from 2–5 years, consider an equal-weighted portfolio of 2- through 5-year bonds with a 1-year holding period. The time $t+12$ excess return of this portfolio is $14(rxt+12(24)+rxt+12(36)+rxt+12(48)+rxt+12(60))$⁠. The loading of this portfolio on tightness is the average of the four loadings of the 2- through 5-year bonds and equals -4.05. This translates into a 1.05% increase in future bond excess returns for a one-standard-deviation decrease in labor market tightness.

As in the data, our model-implied coefficients for $βΘ(n)$ is negative. The model-implied value of $βΘ(n)$ for the equal-weighted portfolio is –2.93 which is 72% of the data counterpart. In addition, our model captures the fact that the magnitude of $βΘ(n)$ is increasing in maturity, although the profile of $βΘ(n)$ as a function of maturity is flatter in the model compared to that of the data.

Row 1 of panel B of Table 5 reports bond excess return predictability results for the job finding rate $ft$⁠. We construct $ft$ following Elsby, Michaels, and Solon (2009). We see that the job finding rate negatively forecasts bond excess returns, with the regression coefficients increasing in magnitude as a function of the maturity of the bond. Both of these patterns are in line with our model’s predictions although the model-implied slope coefficient is higher than its data counterpart—the slope coefficient for the equal-weighted portfolio is –9.44 in our model versus -6.31 in the data. Section A.2.2 of our Internet Appendix shows that the job separation rate does not forecast bond returns. Our modeling assumption of a constant job separation rate in our baseline model reflects this finding. We report results of an extension of our model which relaxes this assumption in  Appendix C.

More broadly, Ludvigson and Ng (2009), Joslin, Priebsch, and Singleton (2014), Huang and Shi (2016, 2023), and Bianchi, Büchner, and Tamoni (2020) show a robust link between bond excess returns and principal components of a wide set of real macroeconomic variables, which include labor market variables, such as employment, unemployment, and vacancies. In particular, Ludvigson and Ng (2009) find bond excess returns to be predicted by a principal component of real macroeconomic variables that has a large loading on employment and hours. Our theory provides an explanation for the ability of labor market variables to predict bond excess returns (as we remark in  appendix 2, further investigating the “spanning hypothesis” is beyond the scope of our paper).

We conclude our empirical analysis with two robustness checks for the predictive regressions (48). First, we rerun our analysis but additionally include the federal funds rate, the supply of outstanding government bonds (measured by the maturity-weighted debt to GDP series from Greenwood and Vayanos 2014), and the level of inflation (based on the CPI index) as controls. Second, we rerun the univariate regression (48) using the excess returns on Treasury Inflation-Protected Securities (TIPS) as the left-hand-side variable. The results for our first robustness check are shown in Table A2 in  Appendix A. The results for our second robustness check are shown in Section A.2.1 of our Internet Appendix. Our findings remain unchanged: labor market tightness and the job finding rate continue to negatively forecast bond excess returns.

3.2 Fama-Bliss regressions

We conclude this section by showing that yield-based variables predict bond excess returns in our model. Specifically, we run the Fama and Bliss (1987) predictive regression in our model which is Equation (48) with the forward-spot spread $Xt=ft(n)−yt(12)$ as the predictive variable. The forward rate is the interest rate between times $t+n−12$ and $t+n$ and is given by $ft(n)=log Pt(n−12)−log Pt(n)$⁠.

Panel A of Table 6 shows results for our baseline model. Row 4 shows that the forward-spot spread predicts bond excess returns in our model. The coefficients are positive and statistically significant as in the data (see row 1). While the magnitude of our model-implied coefficients are smaller than their data counterparts, they are comparable to those in the literature on production-based asset pricing models of the term structure (see, e.g., Jermann 2013; Kung 2015).

We now show that in order to generate the Fama-Bliss results, it is important to account for the slow recovery of labor market variables (e.g., the unemployment rate) recently documented in Hall and Kudlyak (2022). For this purpose, we calibrate a model that is otherwise identical to our baseline model but with $m¯t=∞$ so that the matching function (7) becomes $m(Ut,Vt)=UtVt/(Utι+Vtι)1ι$⁠. This modification of the matching function from our baseline model turns off slow recovery of labor market variables following a switch from the low to the high regime; henceforth we refer to this model as the “model without slow recovery.” We provide details for this model in  Appendix C.

We see from rows 4 and 5 in panel B of Table 6 that the model without slow recovery is unable to generate the Fama-Bliss results: both the loadings $β(n)$ and the $R2$ are essentially zero across maturities 12–60 months. To better understand the role of slow recovery in generating the Fama-Bliss results, it is useful to decompose the forward-spot spread as follows: $ft(n)−yt(12)=Lt(Mt,t+n−12)+Lt(Mt,t+12)−Lt(Mt,t+n)−Et[log Mt+n−12,t+n−log Mt,t+12]$⁠,²⁵ where we have ignored the role of inflation since it plays a secondary role in our setting. Just like bond risk premiums (see Equation (31)), the $Lt(Mt,t+n−12)+Lt(Mt,t+12)−Lt(Mt,t+n)$ term is countercyclical in the baseline model as well as the model without slow recovery and contributes to a positive $β(n)$ for the Fama-Bliss regression.

The difference in the loadings $β(n)$ across the two models arises from the $Et[log Mt+n−12,t+n−log Mt,t+12]$ term. This term measures expected changes in the log SDF over time. In the baseline model with slow recovery, employment grows slowly when the economy transitions out of a recession (see the solid line in panel A of Figure 4) so that the log SDF changes slowly along the transition path (see the solid line in panel B). As a result, the $Et[log Mt+n−12,t+n−log Mt,t+12]$ term is small in magnitude. In contrast, in the model without slow recoveries, employment recovers quickly following a recession (see the dashed line in panel A) so that the log SDF changes rapidly along the transition path (see the dashed line in panel B). As a result, the $Et[log Mt+n−12,t+n−log Mt,t+12]$ term becomes large in magnitude. This, in turn, attenuates the loadings $β(n)$ in the model without slow recoveries.

While it is important to account for slow recovery of labor markets to generate the Fama-Bliss predictability results, the central results of our paper are robust to whether or not we incorporate slow recovery. That is, limited sharing of income risk implies (1) an upward-sloping yield curve and (2) a negative relation between labor market tightness and bond risk premiums in both economies with and without slow recoveries; we include these results in  Appendix C. In Section B.5 of the Internet Appendix, we additionally show that (1) our model captures the ability for the term spread to predict bond excess returns, and (2) similar to the models in Wachter (2006) and Chen (2017), real yields are countercyclical in our setting.

4 Evidence for the Mechanism

In this section, we use estimates of the third central moment of income growth $μ3,t,t+Tinc$ from Guvenen, Ozkan, and Song (2014), henceforth GOS, to provide evidence that its dynamics implies an upward-sloping average yield curve and countercyclical bond risk premiums in the United States.

Mean reversion in $log ζt,t+T$

An upward-sloping average yield curve requires bond risk premiums to be positive on average (see the identity (47)). A positive bond risk premium requires mean reversion of the (log) SDF; this can be seen from the log-normal approximation for the risk premium of a long-term bond (30), $hpxrtH,T≈−Covt(log Mt,t+H, log Mt+H,t+T)$ (we derive this approximation in  Appendix B.2). This requirement of mean reversion of the (log) SDF is true in models with and without perfect risk sharing. However, as discussed in Section 1.3.2, with perfect risk sharing, bond risk premiums depends only on the “Aggregate consumption” term (37a) which leads to the Backus, Gregory, and Zin (1989) puzzle. Instead, under limited risk sharing, the presence of the additional “Income dispersion” (37b) and “Cross-covariance” (37c) terms make it possible to match bond risk premiums without counterfactual implications for aggregate consumption growth if $log ζt,t+T$ mean reverts.

To empirically examine mean reversion in $log ζt,t+T$⁠, we use (i) Equation (42) along with its subsequent discussion which shows that $log ζt,t+T$ is approximately affine in the third central moment of cross-sectional consumption growth $μ3,t,t+T$⁠, and (ii) the assumption of a constant pass-through of income shocks to consumption which implies Equation (44b). Together, (i) and (ii) imply that the requirement of mean reversion in $log ζt,t+T$ translates to the requirement of mean reversion in $μ3,t,t+Tinc$⁠.

The GOS series for $μ3,t,t+Tinc$ are available for 1- and 5-year horizons (ie, $T=12$ and $T=60$ months). This allows us to estimate $Corr(μ3,t,t+12inc,μ3,t+12,t+12ninc)$ for $n=2$ and $n=6$⁠, which are relevant for the risk premiums of 2- and 6-year bonds, respectively. Row 1 of Table 7 reports the results. We see that the correlations are 0.21 and -0.49 for $n=2$ and $n=6$ years, respectively. That is, the third central moment of income growth mean reverts at the 6-year horizon. This implies the “Income dispersion” term positively contributes to the unconditional risk premium of 6-year bonds.

We use two proxies for $μ3,t,t+12ninc$ to estimate $Corr(μ3,t,t+12inc,μ3,t+12,t+12ninc)$ for intermediate maturities $n=3,4,5$⁠. The first proxy $μ^3,t,t+12ninc$ is constructed from the 1-year GOS series $μ3,t,t+12inc$ according to $μ^3,t,t+12ninc≡μ3,t,t+12inc+μ3,t+12,t+24inc+⋯+μ3,t+12(n−1),t+12ninc$⁠. This proxy assumes that idiosyncratic labor income shocks are conditionally independent after conditioning on aggregate variables. The GOS data allows us to check the validity of this assumption for $n=5$ years—we find a correlation of 0.96 between the actual series $μ3,t,t+60inc$ and the proxy $μ^3,t,t+60inc$⁠. Therefore, we deem $μ^3,t,t+12ninc$ to be a reasonable proxy for the unavailable $μ3,t,t+12ninc$ series and use $Corr(μ^3,t,t+12inc,μ^3,t+12,t+12ninc)$ to estimate $Corr(μ3,t,t+12inc,μ3,t+12,t+12ninc)$ for $n=3$⁠, 4, and 5 years. The results are shown in row 2 of Table 7. From this table we see that the third central moment of income growth mean reverts for n between 3 and 5 years, becoming more negative for longer maturities. This implies that the “Income dispersion” term positively contributes to the unconditional risk premiums of 3- to 5-year bonds. The mean autocorrelation averaged over $n∈{2,3,4,5,6}$ years is $−0.15$⁠.

Our second proxy for $μ3,t,t+Tinc$ uses unemployment rate changes between t and $t+T$⁠, $Ut+T−Ut$⁠. Figure 5 verifies that $Ut+T−Ut$ is highly correlated with $μ3,t,t+Tinc$ for the two horizons over which the GOS series is available. The correlation between income skewness and unemployment rate changes over the sample period 1979-2011 is -0.77 and -0.72 over horizons of 1 and 5 years, respectively. This high correlation is not surprising since unemployment risk is a part of total income risk in the GOS series.

Row 3 of Table 7 reports values for $Corr(Ut+12−Ut,Ut+12n−Ut+12)$ over horizons ranging between $n=2$ and $n=6$ years, and we see that the correlation is negative. Once again, this implies that the “Income dispersion” term positively contributes to unconditional bond risk premiums.

Mean reversion in $log ζt,t+T$ over the business cycle

Next, we provide evidence that $log ζt,t+T$ mean reverts at a faster rate following recessions. This implies that the “Income dispersion” term (37b) contributes to a countercyclical bond risk premium.

We use the unemployment rate change proxy for changes in $μ3,t,t+T$⁠, since it is available at a monthly frequency. The low frequency of the annual GOS series makes it and its proxy $μ^3,t,t+T$ unsuitable for this investigation, and therefore we do not use them for this exercise. Table 8 reports the covariance $Covt(Ut+12−Ut,Ut+12n−Ut+12)$ conditional on business cycle conditions at time t. The estimates are for horizons ranging between $n=2$ and $n=6$ years. We see that the covariance is more negative during NBER-recessions. This result holds robustly across all horizons. For instance, the covariance relevant for a 6 year bond is -0.49 when estimated over NBER-recessions and -0.08 over expansions. These results imply a countercyclical contribution to bond risk premiums from the “Income dispersion” term, and this pattern holds across the term structure.

We have thus far provided evidence for the mechanism in the context of the SDF from our model, which is based on the ratio of cross-sectional averages of individuals’ marginal utilities. The same evidence also applies to the SDF in the model of Constantinides and Duffie (1996) which is equal to the cross-sectional average of individuals’ marginal rates of substitution (ie, $Mt,t+n=βn∫01(Ci,t+n/Ci,t)−γdi$⁠). This is because the Constantinides and Duffie (1996) SDF also can be written in the form (28) with a $log ζt,t+n$ term having the exact same form of dependance on $μ3,t,t+n$ as that in Equation (42) (we demonstrate this connection in  Appendix D). Our evidence of mean reversion in $μ3,t,t+T$⁠, and its implications for interest rate dynamics, therefore, also applies to models whose SDF has the same form as the SDF in Constantinides and Duffie (1996).

5 Conclusion

We present a theory in which limited risk sharing of idiosyncratic labor income risk and labor market adjustment costs (endogenously derived from search frictions) play a key role in determining the dynamics of interest rates. In the general equilibrium, the interaction of these two ingredients relates three quantities: bond risk premiums, cross-sectional skewness of labor income growth, and aggregate labor market conditions as measured by labor market tightness and the job finding rate.

Our model rationalizes two patterns of interest rates that are challenging for equilibrium models: an upward-sloping average yield curve and countercyclical bond risk premiums. It predicts labor market tightness and the job finding rate, key variables in labor search models, to negatively forecast bond excess returns. We find supporting evidence for this prediction.

Yield based variables also predict bond excess returns in our model. We show that it is important to account for the recent evidence for a slow labor market recovery following recessions in order to rationalize the Fama and Bliss (1987) predictability regressions in our setting.

Our mechanism relies on mean reversion in cross-sectional skewness of income growth. We find supporting evidence for such mean reversion using data from Guvenen, Ozkan, and Song (2014). As panel data with joint data on individuals’ portfolio holdings and income (or, more generally, characteristics that predict individuals’ income risk) become available, it will be interesting to use that data to provide more direct micro-level estimates of the strength with which imperfect sharing of income risk affects term structure dynamics, or asset prices in general.

Acknowledgement

We thank Ralph Koijen (the editor), two anonymous referees, Jack Bao, George Constantinides, Andres Donangelo, Jack Favilukis, Mike Gallmeyer (discussant), Shiyang Huang, Mete Kilic (discussant), Leonid Kogan, Jack Liebersohn, Yang Liu, Yukun Liu (discussant), Sydney Ludvigson, Andrey Malenko, Thomas Maurer, Jianjun Miao, Francisco Palomino, Juan Passadore, Carolin Pflueger(discussant), Juan Rubio-Ramirez, Andres Schneider (discussant), Andrea Tamoni, Haoxiang Zhu, and seminar participants at the AFA, AFR Summer Institute, Atlanta Fed, Bloomberg, City U of Hong Kong, Cornell, Delaware, EEA-ESEM, European Finance Association, Federal Reserve Board, Florida State, Georgia, HKU, Houston, Labor and Finance Group, Maryland, MFA, Michigan, NASMES, New York Fed, NFA, Office of Financial Research, Penn State, Rochester, Temple, Texas A&M, and the WFA for helpful comments and discussions. The views expressed here are ours and not necessarily those of the Federal Reserve Bank of Atlanta or the Federal Reserve System. Supplementary data can be found on The Review of Financial Studies web site.

Footnotes

Author notes

Authors have furnished an Internet Appendix, which is available on the Oxford University Press Web site next to the link to the final published paper online.

References

Appendix A Empirical Appendix

Data

We obtain nominal bond price data from the Fama-Bliss Discount Bonds series (available from the Center for Research in Security Prices) and from Gurkaynak, Sack, and Wright (2007). We use the former data set to construct yields and holding period log excess returns for nominal bonds with maturities of 5 years and less. We use the latter data set to construct the 10- minus 1-year nominal slope displayed in panel A of Figure 1. We use 10-year TIPS from Gurkaynak, Sack, and Wright (2010) to construct the 10 minus 1-year TIPS slope shown in panel B of Figure 1; this data set is available from 1999 onward. Since this data set does not report TIPS yields for maturities less than 2 years, we proxy for the 1-year real yield using the 1-year nominal yield minus inflation expectations. The latter is measured using the University of Michigan: Inflation Expectation series (MICH), which we downloaded from Federal Reserve Economic Data (FRED). The liquidity adjusted TIPS slope in panel B of Figure 1 additionally corrects for the liquidity of 10-year TIPS following the procedure in Pflueger and Viceira (2016). This procedure uses three controls for the liquidity differential between 10-year TIPS and its nominal counterpart: the relative transaction volume between 10-year TIPS and 10-year nominal bonds, the spread between 10-year off-the-run and on-the-run treasuries, and the 10-year synthetic minus cash spread (ie, difference between inflation swap rates and break-even inflation). We obtain Primary Dealers’ transaction volume data from the Federal Reserve Bank of New York, and download data for on-the-run treasuries and inflation swap rates from Bloomberg.

Our estimates for the monthly inflation process (39) is based on CPI inflation constructed using the Consumer Price Index (CPIAUCSL) series, which we downloaded from FRED. We deflate the Personal Consumption Expenditures (PCE) series using the Personal Consumption Expenditures: Chain-type Price Index (PCEPI). The resultant monthly series is then used to construct the consumption growth innovations (40).

We use the quarterly real personal consumption expenditures per capita (A794RX0Q048SBEA) series, retrieved from FRED, to compute the data moments for aggregate consumption growth reported in Table 2. We measure unemployment rates using the seasonally adjusted civilian unemployment rate series, downloaded from FRED. We measure vacancies using the composite Help Wanted Index from Barnichon (2010), which we downloaded from Regis Barnichon’s website. We construct job finding and separation rates following the procedure in Elsby, Michaels, and Solon (2009). This procedure makes use of the following series from FRED: civilian labor force size (CLF16OV), unemployment level (UNEMPLOY), and number of civilians unemployed for less than 5 weeks (UEMPLT5).

The estimates for moments of equity excess returns in Table 2 are based on data for the annual Fama-French factors which we downloaded from Kenneth French’s website. We calculate the D/P series using CRSP data for the monthly ex and cum-dividend returns of the aggregate stock market; we annualize the monthly D/P ratio by multiplying by 12.

Additional return predictability results

Table A2 reports results for the return predictability regression (48) when the federal funds rate, inflation and maturity weighted debt to GDP are included as additional controls.

Appendix B Model Appendix

B.1 Derivation of the SDF

The SDF in Equation (24) can be obtained by combining the first-order conditions characterizing optimal portfolio choice, $ΛtPtS=Et[∂Jt+1∂φtS]$ and $ΛtPt(1)=Et[∂Jt+1∂φtB]$⁠, with the envelope conditions $∂Jt+1∂φtS=Λt+1(Dt+1+Pt+1S)$ and $∂Jt+1∂φtB=Λt+1$⁠.

The complementary slackness condition $Ψt(ΦtCe,t−Cu,t)=0$ implies that the sum of $Ce,t$ times (B2) and $Cu,t$ times (B3) is equal to $Λt(NtCe,t+UtCu,t)=(NtCe,t1−χ−1+UtCu,t1−χ−1)1−γ1−χ−1.$ Substituting in Equation (21) for aggregate consumption, we obtain $Λt=Ct−γ(Nt(Ce,t/Ct)1−χ−1+Ut(Cu,t/Ct)1−χ−1)1−γ1−χ−1,$ which is Equation (25).

B.2 Approximate Bond Risk Premiums

B.3 Computing Nominal Bond Prices

The inflation process (39) is of the form

The second term on the right-hand side of Equation (B5) is given by $P($,B,n)(Δπ(1),επ)=exp (−a($,B,n)−b($,B,n)Δπ(1)−c($,B,n)επ),$ where the coefficients can be computed recursively: $a($,B,n)=a($,B,n−1)+μπ(1−ρπ)(1+b($,B,n−1))−σπ22(1+b($,B,n−1)+c($,B,n−1))2$⁠, $b($,B,n)=ρπ(1+b($,B,n−1))$⁠, $c($,B,n)=νπ(1+b($,B,n−1))$⁠, starting from the initial condition $a($,B,0)=b($,B,0)=c($,B,0)=0.$

B.4 Derivation of Model-Implied $μ3,t,t+Hmodel$

To derive Equation (45), begin by writing cross-sectional consumption (4) in logs as $cit=ct+ϵit$ where $ct$ is the log of aggregate consumption (21) and $exp(ϵit)=(Φt(1−Nt)+Nt)−1(Φt+(1−Φt)eit)$ with $eit$ being the individual employment shock (3). The third central moment of cross-sectional consumption growth, conditional on the path of aggregate state variables, is equal to the third derivative $μ3,t,t+Hmodel=mΔϵ,t→t+H″′(0)$⁠, where $mΔϵ,t→t+H(α)≡ log Et+H[exp (α[ϵi,t+H−ϵit])]$ is the cumulant-generating function for relative consumption growth $(ci,t+H−cit)−(ct+H−ct)=ϵi,t+H−ϵit$⁠. Using the fact that $eit$ is independent conditional on aggregate variables, we obtain $mΔϵ,t→t+H(α)=mϵ,t+H(α)+mϵ,t(−α)$⁠, where $mϵ,t(α)≡ log Et[exp(αϵit)]=log (Nt+(1−Nt)Φtα)−α log (Nt+(1−Nt)Φt)$⁠. Equation (45) follows from $mΔϵ,t→t+H″′(0)=mϵ,t+H″′(0)−mϵ,t″′(0)$ and $mϵ,t″′(0)=Nt(1−Nt)(2Nt−1)(log Φt)3$⁠.

Appendix C Additional Model Results

In this section, we demonstrate that our baseline model’s results for bond risk premiums continue to hold under alternative model specifications. We also demonstrate the importance of slow recoveries for Fama-Bliss regressions.

C.1 Alternative Model Specifications

No slow recoveries

In this specification, we switch off slow recoveries. This is implemented by ignoring the upper bound on the number of matches that could be formed. That is, $m¯t=∞$ so that the matching function (7) is $m(Ut,Vt)=UtVt/(Utι+Vtι)1ι$ instead.

Additive preference aggregator

AR(1) productivity process.

In contrast to our baseline model, in which productivity takes two values, the AR(1) specification (C1) allows for many more states.²⁶

Time-varying job separation rate

In this specification, we allow the job-separation rate to vary over the business cycle. We do so by allowing $st=s(zt)$ to depend on aggregate productivity. All laws of motion involving job separation rates are modified accordingly (e.g., the upper bound on the number of matches (8) becomes $m¯t=s(zt)+(1−s(zt)−a)Ut$ while the law of motion for employment (12) becomes $Nt+1=(1−s(zt))Nt+g(Θt)Vt$⁠).

C.2 Calibration of Alternative Models

We group parameters into two categories: (1) parameters whose values are identical across model specifications, and (2) parameters whose values are different across model specifications. We describe our choice of parameters for each of these two categories below.

Parameters that are identical across model specifications

To facilitate comparison across models, we fix many parameters at their baseline values in Table 1. These fixed parameters are the curvature of the matching function $ι$⁠, the long-run mean and persistence of the risk sharing process, $x¯$ and $ρx$⁠, respectively, as well as the parameters governing the inflation process (ie, $μπ$⁠, $ρπ$⁠, $ξπ$⁠, $νπ$⁠, and $σπ$⁠).

Parameters that are different across model specifications

Table C1 summarizes the difference in parameter values across models. To facilitate comparison, column 1 reports the parameter values from our baseline model in the main text.

For the “AR(1) productivity process” model, we set $ρz=0.9$ based on the autocorrelation of the labor productivity series in the data,²⁷ and set $σz=0.0063$ to match an unconditional volatility of quarterly output growth of 0.81%. For all other models, we use the transition probabilities from Table 1 and choose the value of productivity in the recessionary regime $log z2$ based on the same target for the unconditional volatility of output growth.

The model with time-varying job separation rates sets $s(z1)=0.037$ and $s(z2)=0.033$ in the low- and high-productivity regimes, respectively. These values are the average of monthly job separation rates over NBER recessions and expansions from 1964 to 2016, respectively. Job separation rates in all other models remain unchanged from that of the baseline model.

When we change our baseline model to the alternative specifications, some targeted moments change substantially. Therefore, we recalibrate a subset of parameters following the same calibration procedure for the baseline model. First, we recalibrate the volatility of the risk sharing process $σx$ so that $σ(μ3,t,t+12model)=0.0156$⁠. Second, we recalibrate $β$ and $χ$ to target the mean and volatility of the 1-year nominal rate. The “additive preference aggregator” model restricts $γ=χ−1$ while all remaining models fix $γ=2$ at its baseline value. Third, we recalibrate b, $κ$⁠, and $η$ to jointly target the first two moments of the unemployment rate and the volatility of the total wage bill relative to that of output.

Moments across model specifications

Table C2 compares the moments of the alternative models to that of the baseline model. Comparing across models, we see that most model-implied moments are near identical.²⁸ The main difference across models lies in their implications for the equity premium. In particular, the model-implied equity premium is smaller in the “no slow recovery” and “additive preference aggregator” models. Other than this difference in the models’ implications for the equity premium, we show in the next section that the models’ implications for term structure of interest rates are very similar.

C.3 Term Structure under Alternative Model Specifications

We find that the alternative models described in  Appendix C.1 all imply upward-sloping yield curves and positive bond risk premiums. We provide details below.

Term structure of interest rates

Panel A of Table C3 reports the average of model-implied yields for 1- through 5-year nominal bonds. Like our baseline model, the yield curve is upward sloping across all alternative model specifications.

Panel B of Table C3 reports the corresponding yield volatilities. The term structure of volatilities for all alternative models display properties similar to that of the baseline model. For example, the volatility term structures are downward sloping across all models.

Panels C and D of Table C3 report the means and volatilities for the real term structure, respectively. The implications for the real yield curve are similar across alternative model specifications.

Decomposition of real bond risk premiums

Table C4 reports the decomposition of real bond risk premiums (defined in Equation (33)) across the alternative model specifications that we consider. The decomposition is for an equal-weighted portfolio of 2- through 5-year real bonds, over a holding period of a year. Although the precise values of the decomposition differ across models, we see that the main finding from the baseline model continues to hold: the income dispersion and cross-covariance terms are responsible for the majority of real bond risk premiums, with the aggregate consumption persistence term playing a small role.

Labor markets and bond risk premiums

Next, we turn to the return predictability regression (48) for the 1-year log holding period excess return of nominal bonds. Table C5 reports the regression coefficient $β(n)$ across models, with the regressor being labor market tightness (⁠$Xt=Θt$⁠) and the job finding rate (⁠$Xt=ft$⁠) in panels A and B, respectively. Similar to the baseline model, labor market conditions negatively predict bond excess returns across all alternative model specifications, with the magnitude of the regression coefficients being increasing in maturity.

Yield-based predictability regressions

Table C6 reports results for the Fama-Bliss predictability regressions across the alternative model specifications that we consider. We see that except for the no slow recoveries and the additive preference aggregator models, the forward-spot spread predicts bond excess returns with similar loadings.

Appendix D Relating $ζt,t+n$ to Cross-Sectional Moments

In this section, we relate the SDF to moments of the cross-sectional distribution of consumption growth. We do so for two classes of SDFs:

• The SDF in our model that is based on the ratio of the cross-sectional average of individuals’ marginal utilities whose $ζt,t+n$ is defined according to Equations (26) and (28).

• The SDF in Constantinides and Duffie (1996) that is based on the cross-sectional average of individuals’ marginal rates of substitution, $Mt,t+n=βn∫01(Ci,t+n/Ci,t)−γdi$⁠, which when written in the form of Equation (28), implies a $ζt,t+n$ of form

  $ζt,t+n≡∫01(Ci,t+n/Ct+nCi,t/Ct)−γdi.$

  (D1)

where the coefficients A and B depend on preference parameters. This expression corresponds to Equation (42). We provide details below.

As an intermediate step to deriving Equation (D2), we begin by relating $ζt,t+n$ to the idiosyncratic shocks $ϵi,t,t+n$⁠.

$ζt,t+n$ defined according to Equations (26) and (28)

$ζt,t+n$ defined according to Equation (D1)

Relating the two forms of $ζt,t+n$ to cross-sectional moments

In both cases (D9) and (D8), $log ζt,t+n$ is related to the cumulant-generating function of $ϵi,t,t+n$⁠, with parameter $1−χ−1$ in Equation (D8) and $−γ$ in Equation (D9). We obtain Equation (D2) by expanding $log ζt,t+n$ to third order in powers of its central moments, where we used the fact that the cross-sectional average over $ϵi,t,t+n$ is zero, the second central moment is the variance $σt,t+n2≡Et[ϵi,t,t+n−Et[ϵi,t,t+n]]2$ and the third central moment is $μ3,t,t+n≡Et[ϵi,t,t+n−Et[ϵi,t,t+n]]3$⁠. The constants in Equation (D2) are different for the two classes of SDFs. $A=γ2/2$ and $B=−γ3/6$ for the $ζt,t+n$ in Equation (D9), while $A=(γ−1)(χ−1−1)/2$ and $B=−(γ−1)(χ−1−1)2/6$ for the $ζt,t+n$ in Equation (D8).

Code Availability: The replication code is available in the Harvard Dataverse at https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/NH6EEF