Robust Bond Risk Premia Predictability Test in The Quantiles
Robust Bond Risk Premia Predictability Test in The Quantiles
Robust Bond Risk Premia Predictability Test in The Quantiles
and
Xinjue Li
School of Economics, Yunnan University
and
Qingliang Fan ∗
Department of Economics, The Chinese University of Hong Kong
October 7, 2024
Abstract
Different from existing literature on testing the macro-spanning hypothesis of
bond risk premia, which only considers mean regressions, this paper investigates
whether the yield curve represented by CP factor (Cochrane and Piazzesi, 2005)
contains all available information about future bond returns in a predictive quantile
regression with many other macroeconomic variables. In this study, we introduce the
Trend in Debt Holding (TDH) as a novel predictor, testing it alongside established
macro indicators such as Trend Inflation (TI) (Cieslak and Povala, 2015), and macro
factors from Ludvigson and Ng (2009). A significant challenge in this study is the in-
validity of traditional quantile model inference approaches, given the high persistence
of many macro variables involved. Furthermore, the existing methods addressing this
issue do not perform well in the marginal test with many highly persistent predictors.
Thus, we suggest a robust inference approach, whose size and power performance
are shown to be better than existing tests. Using data from 1980–2022, the macro-
spanning hypothesis is strongly supported at center quantiles by the empirical finding
that the CP factor has predictive power while all other macro variables have negligible
predictive power in this case. On the other hand, the evidence against the macro-
spanning hypothesis is found at tail quantiles, in which TDH has predictive power
at right tail quantiles while TI has predictive power at both tails quantiles. Finally,
we show the performance of in-sample and out-of-sample predictions implemented by
the proposed method are better than existing methods.
∗
Corresponding author: Q. Fan (E-mail: michaelqfan@gmail.com).
1
1 Introduction
Treasury bonds are a crucial component in numerous investment portfolios. Therefore,
understanding the risk and return dynamics of this asset class is of fundamental importance
from an economic perspective. One focal topic in the study of U.S. Treasury bond returns
is the famous macro-spanning hypothesis/puzzle on whether the yield curve contains all
1
available information about future bond risk premia. Numerous studies have explored
the predictive power of macro variables on bond risk premiums while controlling for the
CP factor (Cochrane and Piazzesi, 2005) representing the information of the yield curve.
These studies include Cooper and Priestley (2008), Ludvigson and Ng (2009), Bansal and
Shaliastovich (2013), Joslin et al. (2014), Greenwood and Vayanos (2014), Cieslak and Po-
vala (2015), Ghysels et al. (2018), Bauer and Hamilton (2018), Zhao et al. (2021). The
literature mentioned above employs mean regression to assess the macro-spanning hypoth-
esis, which suggests that there is no predictability of bond returns beyond the information
provided by the yield curve. However, it is uncertain whether these results are consistent
across the entire distribution or mainly pertain to the mean return. This paper studies the
predictability of bond risk premia in a predictive quantile regression framework. First, the
quantile model explores the heterogeneous bond return predictability, which enrichs the
study of macro-spanning hypothesis. In the literature, Andreasen et al. (2021) and Borup
et al. (2024) have examined the heterogeneity of bond risk premia predictability, taking
into account the business cycle and economic uncertainty. However, the heterogeneous
predictability of bond risk premia at various quantile levels remains unexplored. Second,
quantile model helps predict bond risk premia in different scenarios, e.g., investors often
pay more attention to the tail risk of bonds. Third, the estimation and inference in mean
regression are significantly affected by outliers, whereas the quantile regression is robust
to the outliers. There are notable challenges to conduct the test of spanning hypothesis in
a quantile model. First, Bauer and Hamilton (2018) noted that the predictors commonly
employed in the literature exhibit high persistence, leading to size distortion in predictive
quantile regression test statistics and potentially causing spurious findings. Second, the
conventional test statistics suffer size distortion at tail quantiles due to inaccurate den-
1
See the beginning of Section 3.1 for the mathematical representation of the hypothesis.
2
sity estimator of quantile regression errors. Third, when investigating the macro-spanning
hypothesis, it is common for previous empirical studies to utilize four or more predictors.
These predictors typically include the CP factor, two LN factors (Ludvigson and Ng, 2009),
2
and at least one additional macro predictor. In summary, testing the spanning hypothesis
at the quantiles needs a robust marginal inference approach (meaning testing the predic-
tive power of one variable while controlling for others, as opposed to conducting a joint
test) in a predictive quantile regression framework that includes a large number of highly
persistent predictors. Unfortunately, existing inference methods suffer size distortion in
marginal tests with many highly persistent predictors. In the literature, Lee (2016) ex-
tended the IVX method (Phillips and Magdalinos, 2009) in mean regression to the quantile
regression framework (IVX-QR) to test the predictability at various quantile levels, with
a focus on joint tests but suffers severe size distortion at tail quantiles due to inaccuracy
of nonparametric density estimator. To avoid size distortions at tail quantiles, Fan and
Lee (2019) combined the moving block bootstrap technique and IVX-QR approach for the
joint test, which circumvents estimating this density. Liu et al. (2023) developed a uni-
fied predictability test for univariate predictive quantile regression with highly persistent
predictors. All three methods above are unsuitable for the marginal tests in multivariate
models. Cai et al. (2023) (hereafter referred to as CCL2023) constructed the instrumental
variable (IV) estimator based on a double-weighted method for predictive QR (DW-QR) to
conduct the joint and marginal tests in multivariate models. However, the size performance
of CCL2023 is still not very satisfactory at tail quantiles and with many predictors.
Our new inference procedure fits the task of testing the macro-spanning hypothesis for
bond risk premia. First, we construct a consistent IV estimator under both the null and the
alternative hypotheses by a two-step quantile regression. Nonetheless, the test statistics
constructed by the IV estimator above continue to experience size distortions with many
highly persistent predictors, which is induced by two higher-order terms detailed in Section
3.2. Second, we improve the size and power performance of the test. To address the size
2
See the details of the variables description in Section 2.
3
distortion issue, we adopt the sample splitting method proposed by Liao et al. (2024) to
eliminate one higher-order term. We select a conservative tuning parameter in constructing
instrumental variables to reduce the size distortion induced by another higher-order term.
To amend the power loss due to the conservative tuning parameter choice, we enhance the
power of the aforementioned test by adding a modified test statistic from the conventional
√
test, which converges to zero at the rate T under the null hypothesis. This way, we fully
utilize the good power performance of the traditional test without the bundling of its size
distortion effect. Third, when constructing the final test statistic, we introduce a new and
accurate density estimator of the error term, which utilizes a simulation-based estimator
3
instead of the inefficient nonparametric estimator. Next, we apply the above procedure
for the macro-spanning hypothesis by testing the predictive power of the LN factors, TDH,
and TI at various quantiles while controlling for the CP factor. The main empirical findings
4
are summarized as follows. On the one hand, the CP factor has significant predictive
power, while all macro variables, including the LN factors, have negligible predictive power
at center quantiles, which strongly supports the macro-spanning hypothesis. On the other
hand, strong empirical evidence against the macro-spanning hypothesis is found in tail
quantiles. In particular, TDH has predictive power for bonds of all maturities at right tail
quantiles, while it only has predictive power for 2- and 3-year bonds at left tail quantiles.
Moreover, for bonds of all maturities, TI has significant predictive power at all quantiles
except for center quantiles. Potential economic explanations for the findings are discussed
in detail in Section 4.2. Our contributions are threefold.
2. We provide a novel and reliable inference approach whose size and power performance
3
See CCL2023 for more details about the procedure using an irrelevant (but suitably picked) auxiliary
variable, such that the resulting test do not depend on regressor persistence under the null.
4
In our empirical study, the variable has predictive power if its rejection rate is less than 1%.
4
are significantly better than the literature for marginal tests in predictive quantile
regression with many highly persistent predictors. As a result, our approach is more
suitable to conduct macro-spanning hypothesis than those in existing literature.
The rest of this paper is organized as follows. Section 2 describes the data characteristics
of bond risk premia and its predictors. Section 3 introduces our inference approach in
the predictive quantile regression with highly persistent predictors. Section 4 presents
the heterogeneous predictability of bond risk premia at various quantiles. Section 5 shows
the excellent in-sample and out-of-sample performance compared with CCL2023. Section 6
concludes the paper. The online appendix includes an algorithm for the inference procedure,
additional theoretical and numerical results, and an application on the left and the right
d p d
tail risk indicators. Throughout this paper, the standard notations ⇒, →
−, →
− and = are
used to represent weak convergence, convergence in distribution and in probability, and
equivalence in distribution, respectively. All limits are for T → ∞ in all limit theories, and
Op (1) is asymptotically bounded while op (1) is asymptotically negligible.
2 Data Characteristics
In this section, we illustrate the variables that are used to test the macro-spanning
hypothesis. Following Cochrane and Piazzesi (2005), we refer to the difference of the
yield-to-maturity (YTM) of the n-year and 1-year maturity of U.S. discount bond as the
dependent variable, bond risk premia rx(n), where n = 2, 3, 4, 5. The five predictors
are CP factor, LN factors (LN1 and LN2), and TI, all from previous literature, and trend
in debt holding (TDH), which is new in the literature. The n-year U.S discount bond data
rx(n) is available on Center for Research in Security Prices (CRSP) and the macroeco-
nomic factors are from the paper of Ludvigson and Ng (2009). We compute the CP factor
following Cochrane and Piazzesi (2005). The consumer price index (CPI) and the debt
5
5
holding data are from the CEIC database. The sample used in this study consists of
monthly data from January 1980 to December 2022. We refer to quantile levels close to
0.5 as center quantiles and quantile levels much less (greater) than 0.5 as left (right) tail
6
quantiles. Next, we explain the predictors in detail. First, the CP factor (Cochrane and
Piazzesi, 2005) represents the information of the yield curve, which is the linear combi-
nation of the forward rate Fi,t−1 , i = 1, 2, 3, 4, 5. Second, the LN factors, LN1 and LN2
constructed by Ludvigson and Ng (2009) are fitted values of 14 5n=2 rx(n)t regressing on
P
macroeconomic factors which are the first eight principal components from a large dateset
of 132 macroeconomic indicators. Third, the trend inflation (TI) is proposed by Cieslak and
Povala (2015) to test the predictive power of highly persistent expected inflation dynamics
on bond excess returns. Specifically, TI at period (t-1) is equal to (1 − w̌) t−1 i
P
i=1 w̌ π̌t−i ,
where π̌t−i is inflation rate in the core CPI and w̌ is 0.9 in this paper, w̌i is the ith power
7
of w̌. Besides the above popular predictors in existing literature, we introduce the trend
in debt holding (TDH), defined as the trend in U.S. federal debt held by the public and
by various government agencies (inter-governmental holdings). TDH represents the rate
at which the federal government’s borrowed funds increase in order to address the out-
standing balance of expenses over time. It contains the information of the demand shift for
the national bonds such as treasury bonds, treasury inflation-protected securities (TIPS)
and discount bonds. Specifically, we construct TDH at period (t-1) in the same manner
iˇ ˇ
as TI, which is equal to (1 − w̌) t−1
P
i=1 w̌ dt−i and dt−i is the increment of debt holding
at time t − i. Figure 2 shows that the demand shift for U.S. treasury bonds represented
by TDH frequently precedes fluctuations in bond returns, which indicates the potential
predictability. The statistical and economic significance of the TDH predictor is further
explained in later sections. We first demonstrate the correlation coefficients between bond
risk premia rx(n) and one-period lagged predictors in Figure 1. The size of the colored
square (maximum=1, colors indicate the values of positive/negative correlations) means
the absolute value of the correlation coefficient. It shows that the one-period lagged CP
factor, which only contains the information of yield curve, has the strongest correlations
5
See more details in https://fiscaldata.treasury.gov/americas-finance-guide/national-debt/.
6
We give specific definitions of the left (right) tail quantiles in the illustrative Figure 2.
7
We also set w̌ to other values such as 0.88 and 0.92, and obtain similar simulation and empirical results.
To save space, these results are not shown here but are available upon request.
6
CP LN1 LN2 TI TDH
1
rx(2) 0.8
0.6
0.4
rx(3)
0.2
0
rx(4) −0.2
−0.4
−0.6
rx(5)
−0.8
−1
with bond risk premia while the other macroeconomic predictors have much weaker corre-
lations with bond risk premia. The correlation coefficients only reveal linear relationship
at the center of the distribution. Next, we demonstrate the relationship between bond risk
premia and one-period lagged predictors, specifically TDH and TI, at different quantiles in
their time series plot in Figure 2. For convenience of discussion, we first define the left tail
risk period in Figure 2 as the period in which one of bond risk premia rx(n) is less than or
8
equal to its unconditional quantile at 0.05 level, for n=2,3,4,5. E.g., Jul. 1980–Sep. 1981
is the left tail risk period by above definition. During periods of left tail risk, the difference
between YTM of 1-year and that of n-year bonds (for n=2,3,4,5) is significantly smaller
compared to normal conditions. In the extreme instance of the left tail risk period, such
as Aug. 2022–Dec. 2022, the inverted yield curve (yields on short-term bonds are higher
than those on long-term bonds) occurs. In the left tail quantiles of risk premia for bonds
with maturities of two to five years, the relatively high YTM of 1-year bonds suggests that
investors are concerned about short-term risks in the bond market. This perception likely
stems from current economic indicators or market volatility, which typically leads investors
to demand higher yields for assuming additional risk over the near term. Consequently, in
such market conditions, investors tend to favor bonds with longer maturities (two to five
years) over 1-year bonds. They perceive these longer-term bonds as offering a more at-
tractive balance of risk and return, particularly when short-term forecasts appear unstable.
This preference is especially pronounced at the left tail quantiles, where risk sensitivity is
8
This criterion of 0.05 is not essential for our test in the following sections.
7
heightened. Likewise, we define the right tail risk period as the period in which one of bond
risk premia rx(n) is greater than or equal to its unconditional quantile at 0.95 level, for
n=2,3,4,5. E.g., Jun. 1982–Dec. 1982, is the right tail risk period. Correspondingly, YTM
of n-year bonds (for n=2,3,4,5) is much bigger than that of 1-year bonds in the right tail
risk period. In the right tail quantiles of YTM of n-year bonds (for n=2,3,4,5), there is an
indication that market risks are expected to emerge over the medium term rather than the
immediate future. This higher YTM reflects a risk premium that investors demand due
to anticipated economic uncertainties or rising interest rates affecting these longer matu-
rities. Consequently, in these scenarios, investors tend to prefer 1-year bonds over those
with longer durations (two to five years). This preference is driven by the perceived safety
and greater liquidity of shorter-term bonds, especially when facing potential medium-term
market fluctuations. Some interesting patterns are observable in Figure 2. First, in the left
tail risk periods, e.g., Jul. 1980–Sept. 1981, Dec. 1981–Mar. 1982 and Aug. 2022–Dec.
2022, one-period lagged TI and TDH have an opposite trend of n-year bond risk premia,
for n=2,3,4,5. Second, in the right tail risk periods, e.g., Dec. 2001–Sept. 2002 and Oct.
2008–Dec. 2008, one-period lagged TI and TDH have the same ascending trend of n-year
bond risk premia, for n=2,3,4,5. This suggests that TDH and TI may have predictive
power at tail quantiles. Granted, the correlations of risk premia with predictors shown in
Figures 1 and 2 are insufficient to infer the macro-spanning hypothesis, it motivates us to
test the predictability of bond returns using CP factor, LN1, LN2, TDH and TI, not only
at center quantiles but also at tail quantiles. Table 1 shows that the AR(1) coefficients
of CP factor, LN1, LN2, TDH and TI are close to 1, indicating that these predictors are
highly persistent. To obtain a reliable test about the macro-spanning hypothesis, which
requires marginal test in a multivariate model with many highly persistent predictors, we
present the new procedure in the following section.
8
TDH, TI and 2−year Bond Risk Premia rx(2)
5.0
2.5
0.0
−2.5
1982 1986 1990 1994 1998 2002 2006 2010 2014 2018 2022
−3
1982 1986 1990 1994 1998 2002 2006 2010 2014 2018 2022
−3
1982 1986 1990 1994 1998 2002 2006 2010 2014 2018 2022
Figure 2: TDH, TI and Bond Risk Premia. Gray/Pink area: left/right tail risk period
9
3 Statistical Model
Denote the bond risk premia at period t as yt and its τ th conditional quantile is
Qyt (τ | Ft−1 ) such that Pr [yt ≤ Qyt (τ | Ft−1 ) | Ft−1 ] = τ ∈ (0, 1), where Ft−1 is the in-
formation set available at time t − 1. For simplicity, a linear conditional quantile of yt is
imposed as follows.
where xt−1 = (x1,t−1 , x2,t−1 · · · , xK,t−1 )⊤ is the vector of predictors (K = 5 in this study)
including CP factor, LN factors (LN1 and LN2), TDH and TI, representing Ft−1 . And
βτ = (β1τ , β2τ , · · · , βKτ )⊤ is a K-dimensional coefficient vector given any τ . The macro-
spanning hypothesis essentially means that only the β1τ of the CP factor is nonzero, and
the βiτ ’s of all other macro variables, for i = 2, 3, 4, 5, are zeroes at all quantile levels τ .
We define the quantile measurement error utτ ≡ yt − Qyt (τ |Ft−1 ) and the quantile score
function ψτ (utτ ) ≡ τ − 1(utτ < 0). Since the AR(1) coefficients of predictors in Table 1 are
close to one (highly persistent), it is common in the literature to model these predictors as
follows.
Remark 1. DGP of innovations of bond risk premia yt in (1) includes the conditional
heteroscedasticity case such as GARCH process, which is also described by Liu et al.
(2023). For instance, yt = µτ + x⊤
t−1 βτ + utτ , utτ = ut − σt Qζ (τ ), ut = σt ζt and
9
We also permit xi,t−1 with i = 1, 2, · · · , K to have mixed types of persistence, that is to say, α could
be a function of i. In this setting, the theoretical results for the test statistics are almost the same.
10
Pq Pr
σt2 = µσ + i=1 ai u2t−i + 2
j=1 bj σt−j , where ζt is a sequence of independent and identi-
cally distributed random vectors with means zero and variances one. P [ζt ≤ Qζ (τ )] = τ
implies P (utτ ≤ 0|Ft−1 ) = τ , thus Qutτ (τ |Ft−1 ) = 0 and Qyt (τ |Ft−1 ) = µτ + x⊤
t−1 βτ .
Following Lee (2016) and CCL2023, we impose the following general weakly dependent
structure for innovation {vt } in (2).
P∞
Assumption 1. Assume that vt follows a linear process given by vt = j=0 Fxj εt−j , where
εt is a martingale difference sequence (MDS) with E(εt |Ft−1 ) = 0 and var(εt ε′t |Ft−1 ) = Σε
for Σε > 0 and E∥εt ∥2+ν < ∞ for some ν > 0 and E[ψτ (utτ )εt ] is a constant. Here,
Fx0 = IK , and ∞
P P∞ P∞ j
j=0 j∥Fxj ∥ < ∞ and Fx (1) = j=0 Fxj > 0, where Fx (z) = j=0 Fxj z .
Lee (2016) and CCL2023 state the functional central limit theorem (FCLT) for {ψτ (utτ ), vt }
as follows.
⌊rT ⌋
1 X ψτ (utτ ) Bψτ (r) τ (1 − τ ) Σψτ v
√ ⇒ = BM , (3)
T t=1 vt Bv (r) Σψτ v Ωvv
The conventional test statistics tend to find spurious predictability with SD predictors
and non-zero contemporary correlation between the error term of bond risk premia and
predictors. This occurs because their asymptotic distribution contains nuisance parameters
that cannot be estimated consistently, which is shown in Lee (2016) and CCL2023.
11
3.2 The Inference Procedure and Theoretical Results
We improve IVX-QR (Lee, 2016) to reduce the size distortions arising from the follow-
ing sources: 1, the inconsistency of IV estimator under alternative hypothesis; 2, density
function estimation for the measurement error; 3, the higher-order terms in the test statis-
tic causing the bias. In the following subsection 3.2.1, we develop the consistent estimator.
In subsection 3.2.2, we address the other two issues of size distortions.
To offer an IV estimator that is consistent under both the null and the alternative
hypothesis, we extend the two-step mean regression to quantile regression. We follow Lee
(2016) and define the instrumental variable zt = (z1,t , z2,t · · · , zK,t )⊤ as follows,
where zi,0 = xi,0 − xi,−1 , ∆xi,t = xi,t − xi,t−1 , ρz = 1 + cz /T δ , cz < 0 and 1/2 < δ < 1 (here
we set δ = 0.95). Next, the new two-step regression is conducted as follows.
By this approach, we decompose predictors xt−1 into two orthogonal parts: the fitted
value of (5) x̃t−1 = µ̂x + θ̂zt−1 , and its residual ṽt−1 = xt−1 − x̃t−1 .
where ρτ (u) = u[τ − 1(u < 0)] is referred to as check function in the literature.
The intuition of the above two-step IV estimator is the following. First, by the equations
xt−1 = x̃t−1 + ṽt−1 and (1), it follows that
12
It is clear that both β̂τ and γ̂τ in (6) are consistent estimators of βτ under both the null and
the alternative hypothesis, since they are the estimates of the same true βτ in (7). Second,
the OLS property guarantees that ṽt−1 is orthogonal to x̃t−1 and the intercept term. The
distribution of β̂τ dependents on x̃t−1 rather than ṽt−1 , and thus, it only dependents on zt−1
since x̃t−1 = µ̂x + θ̂zt−1 . Since zt−1 is mildly integrated with SD predictors, the asymptotic
mixture normal distribution of β̂τ is guaranteed. To obtain the asymptotic distribution of
β̂τ , we first establish the Bahadur representation in the following theorem.
Since the first step decomposes predictors xt−1 to orthogonal components, i.e., Tt=1 x̃t−1 ṽt−1
⊤
P
=
PT
0 and t=1 ṽt−1 = 0, the following theorem holds by Proposition 1 and the asymptotic
property of instrumental variable zt−1 (Kostakis et al., 2015).
−1 1
where Ωzz = τ (1 − τ )Ωvv /(−2cz ) and Ωzx = −c−1 ⊤
dJxc (r) Jxc (r) ⊤ for SD
R
z E(vt vt ) − cz 0
13
Equation (8) in Theorem 1 reveals that β̂τ constructed by the new two-step regression
is the instrumental variable estimator in the quantile regression. It is straightforward to
construct the Wald type test statistic Qivx−qr by self-normalization for the null hypothesis
H0 : Rβτ = rτ where R is a J ×K predetermined matrix with rank J, rτ is a predetermined
vector with dimension K.
⊤ h i−1
Qivx−qr = Rβ̂τ − rτ [ β̂τ )R⊤
R Avar( Rβ̂τ − rτ , (9)
⊤ PT PT
where Avar(
[ β̂τ ) = 1
Ω̂−1 Ω̂
f˜uτ (0)2 zx zz
Ω̂−1
zx , Ω̂zx = ⊤
t=1 z̄t−1 xt−1 and Ω̂zz = τ (1−τ ) t=1 z̄t−1
⊤
z̄t−1 . f˜uτ (0) is some consistent estimator of fuτ (0), e.g., the one obtained by the nonpara-
metric method of Lee (2016). Moreover, we define the t-test statistic for the one-sided
marginal test with J = 1 such as H0 : βiτ = 0, for i = 1, 2, ...K.
Rβ̂τ − rτ
Q̌ivx−qr = h i1/2 . (10)
[ β̂τ )R⊤
R Avar(
Proposition 2. Under Assumptions 1 and 2 and the null hypothesis H0 : Rβτ = rτ , one
can show that the limiting distribution of the t-test statistic Q̌ivx−qr with J = 1 and those of
Wald type test statistic Qivx−qr are the standard normal distribution and the χ2 -distribution
with J degrees of freedom, respectively.
Although the asymptotic distributions of Qivx−qr and Q̌ivx−qr are known, they still
suffer from size distortions in finite sample due to the following two reasons. First, both
the size distortions of Qivx−qr and Q̌ivx−qr arise due to two higher-order terms BT and
10
CT similar to those in mean regression models. Second, the nonparametric estimator
f˜uτ (0) is inaccurate at tail quantiles and when the number of macro variables is large. This
inaccuracy causes both Qivx−qr and Q̌ivx−qr to suffer from size distortion in these cases.
We now address the size distortion issues of the preliminary test statistic Qivx−qr and
Q̌ivx−qr in the following aspects. First, we employ a new simulation-based estimator for
the density function of quantile regression error fuτ (0) to reduce the size distortion at tail
10
For more details on the higher-order terms BT and CT , refer to Proposition 3 in the subsequent section
and the mean regression model counterpart discussed in Liao et al. (2024).
14
quantiles and with a large number of predictors. Second, we apply a sample splitting
procedure introduced by Liao et al. (2024) to eliminate one of the higher-order terms
CT . Third, we choose a conservative value for the tuning parameter cz to reduce the size
distortion effect of another higher-order term BT and point out its negative effect on its
power performance. The details are shown below.
The simulation-based approach to estimate fuτ (0) is as follows. We repeat the following
two steps for i = 1, 2, · · · , M1 , where M1 ∈ N+ is a predetermined positive integer.
(i)
Step 1: We randomly generate i.i.d. multivariate time series {ξt }Tt=1 , which follows
N (0M2 , IM2 ) and M2 ∈ N+ .
(i)
Step 2: By (1), it follows that Qyt (τ |Ft−1 ) = µτ + x⊤ ⊤
t−1 βτ + [ξt−1 ] 0M2 . So we are able to
(i)
where ˆlτ is the estimator of 0M2 .
(i)
Since ξt−1 is randomly generated, it is independent of the main sample, then the following
(i)
asymptotic distribution of ξt−1 holds. As the sample size T → ∞,
( T
)−1 T
√ (i) 1 1 X (i) 2 1 X (i)
T ˆlj,τ = [ξ ] √ ξj,t−1 ψτ (utτ ) + op (1) (12)
fuτ (0) T t=1 j,t−1 T t=1
T
1 1 X (i)
= √ ξ ψτ (utτ ) + op (1).
fuτ (0) T t=1 j,t−1
M1 XM2 h
1 X √ (i) i2 P
h i
(1) 1
ˆ → Avar ˆl1,τ =
T lj,τ FT − τ (1 − τ ), (13)
M i=1 j=1 fuτ (0)2
M1 XM2 h
1 X √ (i) i2 P h i
(1) 1
T ˆlj,τ −
→ Avar ˆl1,τ = τ (1 − τ ),
M i=1 j=1 fuτ (0)2
15
as M1 → ∞ and T → ∞ and M2 is constant. Therefore, the simulation-based estimator
of fuτ (0) is defined as follows
( M1 XM2 h
)−1/2
1 X √ (i) i2 P
fˆuτ (0) := T ˆlj,τ [τ (1 − τ )]1/2 −
→ fuτ (0), (14)
M i=1 j=1
of the appendix imply that the accuracy of fˆuτ (0) could be improved by enlarging M. There-
fore, applying fˆuτ (0) in constructing the test statistics could improve the size performance
at tail quantiles and with a large number of predictors.
Higher-order Terms
Another source of size distortion of Qivx−qr and Q̌ivx−qr is the two higher-order terms,
CT and BT , shown in the following Proposition 3. To demonstrate the concept of higher-
11
order terms, we present them using the univariate case of Q̌ivx−qr . A similar conclusion
can be drawn for Qivx−qr since Qivx−qr = Q̌2ivx−qr with J = 1.
16
J = K = 1, we have the following results for SD predictors.
Q̌ivx−qr = ZT + BT + CT + op T (δ−1)/2 ;
T
1 X P
ZT = Ω−1/2
zz zt−1 ψτ (utτ ) −
→ N(0, 1);
T 1/2+δ/2 t=1
T
1 X P
BT = ϖb Ω−1/2
zz zt−1 ψτ (utτ ) −
→ 0; (15)
T 1/2+δ/2 t=1
√
T (1−δ)/2 E (BT ) → −ρvψ / −2cz ; (16)
" T
#−1/2 T T
(1−δ)/2 1 X
⊤ 1 X 1 X
T CT = τ (1 − τ ) zt−1 zt−1 zt−1 √ ψτ (utτ ) (17)
T 1+δ t=1
T 1/2+δ t=1 T t=1
−1/2
p
where ρvψ = Σvv Σψτ v / τ (1 − τ ) is the correlation coefficient between ψτ (utτ ) and vt and
h i
−1/2 1 PT ⊤ −1/2
ϖb = − 21 τ (1 − τ )Ωzz T 1+δ z z Ω
t=1 t−1 t−1 zz − 1 .
1
PT
The slow convergence of T t=1 zt−1 in Q̌ivx−qr induces CT , while the correlation be-
tween the numerator and the denominator of Q̌ivx−qr leads to BT . To eliminate CT , we
follow Liao et al. (2024) to define a new instrumental variable z̃t−1 based on IV zt−1 in (4).
Specifically, define z̃t−1 = (IK −Sa )zt−1 when 1 ≤ t ≤ T0 = λ T 12 and z̃t−1 = (IK −Sb )zt−1
−1
when T0 + 1 ≤ t ≤ T , where Sa = T1 Tt=1 zt−1 T10 Tt=1
P P 0 ⊤ 1 PT0 ⊤ 1 PT0
zt−1 T0 t=1 zt−1 T0 t=1 zt−1
−1
and Sb = T1 Tt=1 zt−1 T −T
P 1
PT ⊤
1 PT ⊤ 1
PT
0 t=T0 +1 zt−1 T −T0 t=T0 +1 zt−1 T −T0 t=T0 +1 zt−1 . Since
the mean of the IV z̃t−1 is exactly zero, the source of CT is eliminated. Hereafter, we
apply the IV z̃t−1 to run the two-step regression introduced in (5) and (6) to obtain the
consistent IV estimator of βτ under both the null and the alternative hypothesis. Specially,
we run regression (µ̂lx , θ̂l ) = arg minµx ,θ Tt=1 (xt−1 − µx − θz̃t−1 )⊤ (xt−1 − µx − θz̃t−1 ) to
P
obtain the fitted value x̃lt−1 = µ̂lx + θ̂l z̃t−1 and the residual ṽt−1
l
= xt−1 − x̃lt−1 in the first
⊤
step, and µ̂lτ , (β̂τl )⊤ , (γ̂τl )⊤ = arg minµτ ,βτ ,γτ Tt=1 ρτ yt − µτ − βτ⊤ x̃lt−1 − γτ⊤ ṽt−1
P l
in the
second step. Notice β̂τl uses the full sample and will not lose finite sample efficiency. The
following theorem gives the asymptotic distribution of β̂τl .
12
Following Liao et al. (2024), we set λ = 0.5 to evenly split the full sample in half.
17
Theorem 2. Under Assumptions 1 and 2, for SD and WD predictors, it follows that
DT β̂τl − β = DT β̂τl0 − β + op (1) (18)
T
!−1 T
1 X X
= DT−1 z̃t−1 x⊤ −1
t−1 DT DT−1 z̃t−1 ψτ (utτ ) + op (1).
fuτ (0) t=1 t=1
h i
d d
→
− Zβ = MN 0, Avar(β̂τl )
⊤
where Avar(β̂τl ) = 1
fuτ (0)2
(Σzx )−1 Σzz [Σ−1 ⊤
zx ] , Σzz = λ(IK −S̃a )Ωzz (IK −S̃a ) +(1−λ)(IK −S̃b )
Rλ c R1 c
Ωzz (IK −S̃b )⊤ , Σzx = −c−1 c ⊤ −1 c ⊤
z (IK −S̃a ) 0 dJx (r) Jx (r) −cz (IK −S̃b ) λ dJx (r) Jx (r) −cz [IK
−1
−λS̃a −(1−λ)S̃b ] E(vt vt⊤ ) for SD predictors and Σzx = IK −λS̃a −(1−λ)S̃b E xt−1 x⊤
t−1 for
−1
WD predictors. Sa ⇒ S̃a and Sb ⇒ S̃b , where S̃a = Jxc (1)Jxc (λ)⊤ Jxc (λ)⊤ Jxc (λ)
and S̃b =
−1
Jxc (1) [Jxc (1) − Jxc (λ)]⊤ [Jxc (1) − Jxc (λ)]⊤ [Jxc (1) − Jxc (λ)]
for SD predictors and S̃a =
−1
Bv (1)Bv (λ)⊤ Bv (λ)⊤ Bv (λ) and S̃b = Bv (1) [Bv (1) − Bv (λ)]⊤ {[Bv (1) − Bv (λ)]⊤ Bv (1)
Then the test statistic Ql−qr is constructed for the null hypothesis H0 : Rβτ = rτ ,
⊤ h i−1
Ql−qr ≡ Rβ̂τl − rτ [ β̂τl )R⊤
R Avar( Rβ̂τl − rτ , (19)
over, we construct the t-test statistic Q̌l−qr when J = 1 for right-sided test H0 : βi = 0 vs
Har : βi > 0 and left-sided test H0 : βi = 0 vs Hal : βi < 0.
Rβ̂τl − rτ
Q̌l−qr ≡ q . (20)
l
R Avar(β̂τ )R
[ ⊤
Next, the local power of the test statistics Q̌l−qr and Ql−qr are shown as follows.
18
Proposition 5. Under Assumptions 1 and 2 and the local alternative hypothesis Ha :
Rβτ −rτ = bτ /DT and bτ is a constant vector with dimension J, for SD and WD predictors,
as M1 → ∞ and T → ∞, it follows that
h i−1/2
d d
− Q̌∗l−qr = N(0, 1) + R Avar(β̂τl )R⊤
Q̌l−qr → bτ , J = 1; (21)
h i−1
d d
Ql−qr →− Q∗l−qr = χ2J + 2ZQ + b⊤
τ R Avar( β̂ l
τ )R ⊤
bτ , J ≥ 1, (22)
h i−1
where ZQ = Zβ⊤ R Avar(β̂τl )R⊤ bτ .
By the definition Ωzz = τ (1 − τ )Ωvv /(−2cz ) and Theorem 2, it follows that Avar(β̂τl ) =
⊤
−cz 2fτ (1−τ )
uτ (0)
∗ −1 ∗
2 (Σzx ) Σzz [(Σ∗zx )−1 ] is proportional to −cz with SD predictors, where Σ∗zz =
Rλ
λ(IK −S̃a )Ωvv (IK −S̃a )⊤ +(1−λ)(IK −S̃b )Ωvv (IK −S̃b )⊤ and Σ∗zx = (IK −S̃a ) 0 dJxc (r) Jxc (r)⊤ +
R1
(IK −S̃b ) λ dJxc (r) Jxc (r)⊤ + [IK −λS̃a − (1 − λ)S̃b ] E(vt vt⊤ ). Thus, the larger the magnitude
of |cz |, the smaller the absolute values of the test statistics Q̌l−qr and Ql−qr under the local
alternative hypothesis Ha , resulting in reduced power for both Q̌l−qr and Ql−qr . On the
other hand, Proposition C.1 in the appendix shows that while the higher-order term CT
vanishes in both Q̌l−qr and Ql−qr , the higher order term BTl , as defined in Proposition C.1 of
the appendix persists for the same reasons as BT . Moreover, the smaller the magnitude of
|cz |, the more pronounced the size distortions induced by BTl becomes. Thus, we mitigate
the size distortion induced by the higher-order term BT by selecting a conservative tuning
parameter cz = −8 − 2K. However, Proposition 5 suggests that this conservative setting
of cz improves the size performance at the expense of reduced power.
To amend the loss of power due to the conservative choice of turning parameter cz , we
construct the following linear combination test statistic Q̌m−qr when J = 1 for one-sided
marginal test.
1/(1−δ)
1 Q̌o−qr
Q̌m−qr = Q̌l−qr + √ sign(Q̌o−qr ), (23)
T q̌0.999
where q̌0.999 is 99.9% quantile of the standard normal distribution and 1/(1 − δ) = 20
here. And Q̌o−qr is the conventional test statistic defined as Q̌o−qr = (Rβ̂τo − rτ ) τ (1 −
−1 ⊤ −1/2 ⊤
fˆuτ (0), where x̄t−1 = xt−1 − T1 Tt=1 xt−1 and µ̂oτ , (β̂τo )⊤ =
PT ⊤
P
τ )R t=1 x̄t−1 x̄t−1 R
19
PT
ρτ yt − µτ − βτ⊤ xt−1 . We put q̌0.999 here to ensure that, under the null
arg minµτ ,βτ t=1
hypothesis, most realizations of |Q̌o−qr | remain below one. Consequently, this facilitates
1/(1−δ)
1 Q̌o−qr
the second term T q̌0.999
√ sign(Q̌o−qr ) in maintaining the good size performance of
Q̌m−qr in finite sample without exacerbating it. Note that
Remark 2. Theorem 3 shows that Q̌m−qr not only inherits the good size performance of
Q̌l−qr but also the good power performance of Q̌o−qr . The two key factors contributing to
the good size and the good power performance of Q̌m−qr are outlined as follows. First, under
the null hypothesis H0 : Rβτ = rτ , Q̌l−qr plays the key part in Q̌m−qr while Q̌o−qr converges
√
to zero with a rate T . As result, the higher-order term of Q̌m−qr is approximately equal
to that of of Q̌l−qr under the null hypothesis H0 : Rβτ = rτ and thus inherits the good size
1/(1−δ)
performance of Q̌l−qr . Second, both Q̌l−qr and Q̌q̌0.999
o−qr
sign(Q̌o−qr ) play a role under
the local alternative hypothesis. Thus Q̌m−qr inherits the good power performance of Q̌o−qr .
In summary, the weighting technique employed in Q̌m−qr enables it to achieve both good
size and power performance simultaneously.
Next, we construct a linear combination test statistic Qm−qr for two-sided tests using
20
the similar procedure to that used for constructing Q̌m−qr .
1/(1−δ)
1 Qo−qr
Qm−qr = Ql−qr + √ , (26)
T q0.999
where q0.999 is 99.9% quantile of the χ2J distribution and Qo−qr is the conventional test statis-
−1 ⊤ −1
(Rβ̂τo − rτ )⊤ fˆuτ (0)−2 .
PT ⊤
tic defined as Qo−qr ≡ (Rβ̂τo − rτ ) τ (1 − τ )R t=1 x̄ t−1 x̄ t−1 R
Note that
d R1 −1 −1 ∗ ⊤ ∗∗ d R1
where Q∗o−qr = βτ∗ τ (1−τ )R 0 J¯xc (r)J¯xc (r)⊤ R⊤ βτ , Qo−qr = bτ τ (1−τ )R 0 J¯xc (r)
−1 −1 d ¯x (r)J¯xc (r)⊤ −1 R⊤ −1 b⊤
R1 c
J¯xc (r)⊤ R⊤ b⊤ ∗∗∗ ∗
τ and Qo−qr = βτ τ (1 − τ )R 0
J τ with SD pre-
d d −1 d
dictors; Q∗o−qr = χ2J , Q∗∗ −1 ⊤
b⊤ ∗∗∗ ∗
o−qr = bτ τ (1 − τ )R Var(xt−1 ) R τ and Qo−qr = βτ τ (1 −
−1
τ )R Var(xt−1 )−1 R⊤ b⊤ τ for WD predictors. We can then establish the following theorem,
as M1 → ∞ and T → ∞.
As indicated in Remark 2, Theorem 4 shows that Qm−qr not only retains the favorable
size performance of Ql−qr but also adopts the good power performance of Qo−qr .
21
for our test results in Section 4.2. We compare our results with those of CCL2023, the only
feasible method in existing literature that is designed to do marginal test in the quantile
predictive model with many highly persistent predictors.
22
CP (Qm−qr) CP (CCL2023)
0
−2
−4
−6
0
−2
−4
−6
TI (Qm−qr) TI (CCL2023)
0
−2
−4
−6
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
τj
23
has significant predictive power at right tail quantiles for bonds of all maturities, a finding
not reported by CCL2023. This capability of TDH is particularly valuable for predicting
tail risks associated with bond risk premia. Second, we find that TI has predictive power at
both tails quantiles for bonds of all maturities, while CCL2023 fails to find the predictive
power of TI at right tail quantiles for 4- and 5-year bonds. Third, we find the significant
predictive power of TDH at left tail quantiles for 2- and 3-year bonds, which CCL2023 does
not find. Considering the simulations presented in the appendix, which demonstrate that
our size and power performance surpass those of CCL2023 across all quantiles, the most
plausible explanation for the discrepancies between our findings and those of CCL2023 is
that CCL2023 likely commits type II errors (failing to detect true predictive power) at
tail quantiles, whereas our analysis does not. The predictability at tail quantiles (high-risk
region) is practically beneficial, e.g., for constructing the tail risk indicator in the appendix
Section D.2.
The macro-spanning hypothesis states that the yield curve contains all available in-
formation about future bond risk premia, which provides economic interpretation about
the empirical findings at the center quantiles. The CP factor contains the information of
the yield curve so it has significant predictive power while LN1 and LN2 have negligible
predictive power at all quantiles. Moreover, the spanning hypothesis also explains the test
results that all other variables, including TDH and TI, have no predictive power at center
quantiles. The most interesting empirical finding of this paper is the significant predictive
power of TDH at right tail quantiles and the significant predictive power of TI at both tails
quantiles for bonds of all maturities. The predictive power of TDH at left tail quantiles for
2- and 3-year bonds are also found. Subsequently, we offer potential economic explanations
for these findings. Since investors prefer 1-year bonds over n-year bonds (n=2,3,4,5) at
15
right tail quantiles, this preference results in a more pronounced increase in demand
for 1-year bonds compared to n-year bonds when the TDH rises. Consequently, as TDH
increases, the price of 1-year bonds experiences a sharper increase relative to the prices
15
See more discussions of this point in Section 2.
24
of n-year bonds (n=2,3,4,5). This sharper price increase leads to a more substantial de-
crease in the YTM of 1-year bonds compared to that of n-year bonds. Given that the risk
premium of n-year bonds (n=2, 3, 4, 5) is defined as the difference between the YTM of
n-year and 1-year bonds, it follows that the risk premia for n-year bonds will increase if
TDH rises at the right tail quantiles. In other words, TDH has positive predictive power for
bonds of all maturities at right tail quantiles. On the other hand, when investors exhibit
a preference for bonds with maturities of two and three years over 1-year bonds at the
left tail quantiles, the demand for these 2- and 3-year bonds significantly increases. This
heightened demand occurs particularly if the TDH rises. Consequently, the prices of 2- and
3-year bonds rise more sharply compared to 1-year bonds under the same conditions. This
increase in price leads to a more pronounced decrease in the YTM of 2- and 3-year bonds
than that observed in 1-year bonds. As a result, the bond risk premia, which represent
the difference between the returns of 2- and 3-year bonds and 1-year bonds, decrease if the
TDH increases at the left tail quantiles. This effect highlights the sensitivity of bond risk
premia to shifts in investor preference and market conditions at the left tail quantiles. In
other words, TDH has negative predictive power for both 2- and 3-year bonds at left tail
quantiles. Besides that, the demand for n-year (n=2,3,4,5) and 1-year bonds may rise com-
parably if TDH increases within the center quantiles. Additionally, at the left tail quantiles
of the risk premia for bonds with longer maturities (n=4,5 years), the risk associated with
1-year bonds is high, while bonds with maturities of four and five years continue to face
elevated risks due to their longer durations. Consequently, if the TDH increases at the left
tail quantiles, the demand for both 1-year bonds and longer-term bonds (n=4,5 years) may
rise comparably. As a result, the prices and YTM for both 1-year and longer-term bonds
(n=4,5 years) may move in the same direction and by approximately the same magnitude,
leading to relatively stable bond risk premia. That is, TDH has no predictive power at
center quantiles for bonds of all maturities and left tail quantiles for 4- and 5-year bonds.
Next, we discuss the predictive power of TI found at all quantiles except center quantiles.
If TI increases, the real return of investors decrease, leading to a reduced demand for all
1- to 5-year bonds. However, the demand for bonds with maturities of two to five years
decreases more sharply than that for 1-year bonds when TI increases since investors prefer
25
16
1-year bonds at right tail quantiles. As a result, the prices of n-year (n=2,3,4,5) bonds
decrease more sharply than that of a 1-year bond when TI increases. Thus, the YTM of
these longer-term bonds rise more significantly than those of 1-year bonds. Given that
the risk premia for bonds with maturities of two to five years are defined as the difference
between their YTM and that of 1-year bonds, the risk premia for these longer-term bonds
(n=2, 3, 4, 5) increase as TI rises. That is, TI has positive predictive power for all bonds
at right tail quantiles. Following the same logic, TI has negative predictive power for all
bonds at left tail quantiles.
16
See more discussions of this point in Section 2.
26
quantile τj , where j = 1, 2, · · · , J, which are defined in the beginning of Section 4. Finally,
we obtain in-sample prediction value ŷtτj = µ̇ ˆτj + β̇ˆτ⊤ ẋjt−1 and the in-sample error estimator
j
ûtτj = yt − ŷtτj . Then the qw-CRPS, which is calculated by applying weights to ŷtτ for
J − 1 different quantiles, is computed as follows
J−1 T
2 X 1X
qw −CRPSt = W (τj ) ρτj (ûtτj ), qwc = qw −CRPSt (29)
J − 1 j=1 T t=1
where W (τj ) is a deterministic weighting function. The considered weighting functions are:
W (τj ) = τj (1 − τj ), which gives more weight to the center of the predictive distribution;
W (τj ) = (1 − τj )2 (and W (τj ) = τj2 ), which assigns greater weight to the left (right) tail to
emphasize downside (upside) risks; and W (τj ) = (2τj − 1)2 which places emphasis on the
prediction performance at both tails. Using the qw-CRPS criteria defined in (29), we obtain
the following in-sample performance of our proposed approach and CCL2023 in Table 2.
17
It shows that the in-sample performance of the proposed test Qm−qr is much better
than CCL2023 at center quantiles, left tail quantiles, right tail quantiles and both tails
quantiles, since all of its qwc values are much less than those of CCL2023. The excellent
in-sample performance of Qm−qr arises from its success of finding the significant predictive
power of TDH and TI at right tail quantiles. This result is also shown in Table 2, that
the difference between qwc of Qm−qr and CCL2023 when the weight W (τj ) emphasizing in
right tail quantiles is much larger than those in left tail quantiles and center quantiles.
Table 2: In-sample Performance qwc (in 10−3 ) in 1980-2022
27
and Tm ≤ T̃ ≤ T , we first use the sample during period 1 − (T̃ − 1) to run the follow-
ˆτj (T̃ ), β̇ˆτj (T̃ )⊤ ⊤ =
ing quantile regression to obtain the real time coefficient estimator µ̇
PT̃ −1
arg minµτj ,βτj t=1 ρτj (yt − µτj − β̇τ⊤j ẋjt−1 ), where ẋjt−1 is the one used in (29). Then the
out-of-sample predicted value of bond risk premia at quantiles τj is ŷτoj (T̃ ) = µ̇ ˆτj (T̃ ) +
β̇ˆτj (T̃ )⊤ ẋjT̃ −1 and the out-of-sample quantile regression error is ûoτj (T̃ ) = yT̃ − ŷτoj (T̃ ). The
out-of-sample performance in the period from T0 to T is evaluated by qw-CRPS. Specifi-
2
PJ−1 o
cally, the out-of-sample evaluation criteria are qw −CRPST̃ = J−1 j=1 W (τj ) ρτj [ûτj (T̃ )],
and qwc (Tm ) = T −T1m +1 TT̃ =Tm qw −CRPST̃ . By this criteria, we obtain the out-of-sample
P
The out-of-sample performance of Qm−qr is better than those of CCL2023, as qwc (T0 )
of Qm−qr is much less than those of CCL2023 in all quantiles except for rx(2) at center
quantiles and left tail quantiles where they are equivalent. The excellent out-of-sample
performance of Qm−qr also arises from its success of finding the significant predictive power
of TDH and TI at right tail quantiles. The difference between qwc (T0 ) of Qm−qr and
CCL2023 is much larger at right tail quantiles than those at left tail quantiles, center
quantiles and both tails quantiles.
6 Conclusion
We examine the predictability of bond risk premia using predictive quantile regression
with many highly persistent predictors, employing a reliable inference procedure tailored for
this analysis. Our results lend support to the macro-spanning hypothesis at center quan-
tiles, while evidence at tail quantiles contradicts this hypothesis. This dichotomy helps
clarify inconsistencies observed in previous studies that utilized mean regressions. Specifi-
cally, all predictors, with the exception of the CP factor, exhibit negligible predictive power
28
at the center quantile. Additionally, TDH demonstrates predictive power at the right tail
quantiles, and TI shows predictive power at both left and right tail quantiles. Furthermore,
we demonstrate that both in-sample and out-of-sample prediction performance using our
proposed inference method significantly outperforms existing methods.
References
Andreasen, M. M., Engsted, T., Møller, S. V., and Sander, M. (2021). The Yield Spread and
Bond Return Predictability in Expansion and Recessions. Review of Economic Studies,
34(6):2773–2812.
Bauer, M. D. and Hamilton, J. D. (2018). Robust Bond Risk Premia. Review of Financial
Studies, 31(2):399–448.
Baumeister, C., Huber, F., and Marcellino, M. (2024). Risky oil: It’s all in the tails.
Working Paper 32524, National Bureau of Economic Research.
Borup, D., Eriksen, J. N., Kjær, M. M., and Thyrsgaard, M. (2024). Predicting Bond
Return Predictability. Management Science, 70(2):931–951.
Cai, Z., Chen, H., and Liao, X. (2023). A new robust inference for predictive quantile
regression. Journal of Econometrics, 234(1):227–250.
Cieslak, A. and Povala, P. (2015). Expected Returns in Treasury Bonds. Review of Finan-
cial Studies, 28(10):2859–2901.
Cochrane, J. H. and Piazzesi, M. (2005). Bond Risk Premia. American Economic Review,
95(1):138–160.
Cooper, I. and Priestley, R. (2008). Time-Varying Risk Premiums and the Output Gap.
Review of Financial Studies, 22(7):2801–2833.
29
Fan, R. and Lee, J. (2019). Predictive Quantile Regressions under Persistence and Condi-
tional Heteroskedasticity. Journal of Econometrics, 213(1):261–280.
Ghysels, E., Horan, C., and Moench, E. (2018). Forecasting through the Rearview Mir-
ror: Data Revisions and Bond Return Predictability. Review of Financial Studies,
31(2):678–714.
Greenwood, R. and Vayanos, D. (2014). Bond Supply and Excess Bond Returns. Review
of Financial Studies, 27(3):663–713.
Joslin, S., Priebsch, M., and Singleton, K. J. (2014). Risk Premiums in Dynamic Term
Structure Models with Unspanned Macro Risks. Journal of Finance, 69(3):1197–1233.
Kostakis, A., Magdalinos, T., and Stamatogiannis, M. (2015). Robust Econometric Infer-
ence for Stock Return Predictability. Review of Financial Studies, 28(5):1506–1553.
Lee, J. (2016). Predictive Quantile Regression with Persistent Covariates: IVX-QR Ap-
proach. Journal of Econometrics, 192(1):105–118.
Liao, X., Li, X., and Fan, Q. (2024). Robust Inference for Multiple Predictive Regressions
with an Application on Bond Risk Premia. arXiv preprint arXiv:2401.01064.
Liu, X., Long, W., Peng, L., and Yang, B. (2023). A Unified Inference for Predictive Quan-
tile Regression. Journal of the American Statistical Association, 119(546):1526–1540.
Ludvigson, S. C. and Ng, S. (2009). Macro factors in bond risk premia. Review of Financial
Studies, 22(12):5027–5067.
30
Xiao, Z. (2009). Quantile Cointegrating Regression. Journal of Econometrics, 150(2):248–
260.
Zhao, F., Zhou, G., and Zhu, X. (2021). Unspanned Global Macro Risks in Bond Returns.
Management Science, 67(12):7825–7843.
31
Supplementary Material to “Robust Bond Risk
Premia Predictability Test in the Quantiles”
In this supplementary material, we provide the following parts: Appendix A gives the
details of computation algorithm to construct our test statistics used in the main text.
Appendix B collects all the results for Monte Carlo simulations. The numerical results
show that it is necessary to use the new method introduced in the main text Section 3 for
the macro-spanning hypothesis test. Appendix C includes additional asymptotic properties
for the higher-order terms that appear in the test statistics construction, whose results we
applied in the main text Section 3. Appendix D provides additional empirical results. A
new tail risk indicator is also introduced there, which further showcases the usefulness of
our test. Appendix E presents selected proofs for main theorems.
Explanation about subscript: For any matrix A, A1/2 is a matrix defined as A1/2 =
1/2 1/2 1/2
QL diag(λ1 , λ2 , · · · , λK )Q⊤
L , in which QL is the matrix whose ith column is the
1
Algorithm 1 The instructions for the construction of the test statistics Q̌m−qr and Qm−qr
• Use the simulation-based estimator fˆuτ (0) for density fuτ (0) in (14), which is more accurate than
that of CCL2023.
• Set the conservative tuning parameter cz = −8 − 2K to reduce the size distortion induced by the
higher-order term BT .
6: Amend the loss of power due to the conservative tuning parameter cz = −8 − 2K while keeping a good
1/(1−δ)
Qo−qr
size performance: Construct Qm−qr = Ql−qr + √1T q0.999 in (23) for two-sided test.
7: if Focus on the one-sided marginal test H0 : βi = 0 vs Har : βi > 0 and H0 : βi = 0 vs Hal : βi < 0
then
1/(1−δ)
Q̌o−qr
8: Construct the t-test statistic Q̌m−qr = Q̌l−qr + √1 sign(Q̌o−qr ).
T q̌0.999
9: end if
2
for two-sided test with multiple predictors. Moreover, we show the power performance of
Qm−qr and Q̌m−qr and the size performance of Q̌m−qr in multivariate models. Additionally,
we show the result of the joint test by Lee (2016) in multivariate models. In all experi-
ments, the DGP of xt−1 is equation (2) while the DGP of yt is the same as equation (1) of
Liu et al. (2023), which is described in Remark 1. We set µτ = 1 and µσ = 1. The DGP
of the error vi,t is as follows: vi,t = γi ζt + v̌i,t , where (ζt , v̌i,t )⊤ ∼ i.i.d. N(0K , IK ). Thus
the contemporaneous correlation coefficient between ζt and vi,t is γi (1 + γi2 )−1/2 , which is
the source of the size distortion shown in Proposition 3. We report the simulation results
for a i.i.d. model with σt = 1 for t = 1, 2, · · · , T and the sample size T = 750 with
18
the nominal size 5%. Simulation is repeated 5000 times for each setting. Hereafter,
Lee (2016), Fan and Lee (2019) and Liu et al. (2023) are referred to as L2016, FL2019
and LLPY2023, respectively. Experiment A1: In experiment A1, we conduct the test
in univariate model with γi = −3 such that corr(ut , vt ) = −0.95. First, the size per-
formance (in %) of two-sided test in univariate model of L2016, FL2019, CCL2023 and
LLPY2023 and the proposed test statistic Qm−qr for two-sided test H0 : βτ = 0 are shown
in Table B.1, in which τ = 0.05, 0.25, 0.5, 0.75, 0.95 and predictors are SD (c = 0), SD
(c = −5), SD (c = −15) and WD (c = −0.05). Second, the power performances of L2016,
FL2019, CCL2023 and LLPY2023 and the proposed test statistic Qm−qr for two-sided test
H0 : βτ = 0 vs Ha : βτ ̸= 0 are shown in Figure B.1, in which τ = 0.5 and predictors
are SD (c = 0), SD (c = −5), SD (c = −15) and WD (c = −0.05). Third, the size
performance (in %) and the power performance of the proposed test statistic Q̌m−qr for
one-sided test H0 : βτ = 0 vs Har : βτ > 0 in univariate model are shown in Table B.2,
in which τ = 0.05, 0.25, 0.5, 0.75, 0.95 and predictors are SD (c = 0), SD (c = −5), SD
(c = −15) and WD (c = −0.05). We have the following findings from Tables B.1 and B.2
18
The absolute value of the correlation coefficient between utτ and vi,t in i.i.d. model is significantly
greater than that in the GARCH model. This difference leads to a larger size distortion in all approaches
when compared to their counterparts in the GARCH model. Consequently, it is sufficient to demonstrate
the advantages of the proposed test statistics using only the i.i.d. model results, rather than including
those from the GARCH model. To maintain brevity and convey the main message, we have omitted the
results for the GARCH and ARCH models with sample sizes T = 250, T = 500, and T = 750, as well as
for the i.i.d. model with sample sizes T = 250 and T = 500. Detailed results for these models are available
upon request.
3
and Figure B.1. First, the size performance of the proposed test statistic Qm−qr surpasses
that of L2016, FL2019 and CCL2023 and is slightly better than LLPY2023. The proposed
test statistic Qm−qr is free of size distortion even at tail quantiles. Meanwhile, LLPY2023
tends to slightly over-reject when τ = 0.05 and under-reject when τ = 0.95. Second, the
power performance of the proposed test statistic Qm−qr is better than L2016, FL2019 and
CCL2023 while is comparable to LLPY2023. In scenarios SD (c = 0) and SD (c = −5),
the power of the proposed test statistic Qm−qr is lower than LLPY2023 when βτ is small.
Conversely, Qm−qr exhibits greater power than LLPY2023 when βτ is large. For case SD
(c = −15), the power of the proposed test statistic Qm−qr is the same as LLPY2023. For
case WD (c = −0.05), the power of the proposed test statistic Qm−qr is greater than that of
LLPY2023. Third, the size and the power performance of the proposed test statistic Q̌m−qr
is very well for one-sided test. In summary, the performance of the proposed test Q̌m−qr
is comparable to LLPY2023 and superior to L2016, FL2019 and CCL2023 for two-sided
test. Meanwhile, the proposed test Q̌m−qr performs well for one-sided test, which is not
considered in the literature.
Note that the black dash line in Figure B.1 represent 5% nominal size.
SD (c=0) SD (c=−5)
1.00 1.00
0.75 0.75
0.50 0.50
0.25 0.25
Reject Rate
0.00
0.00 0.01 0.02 0.03 0.00 0.01 0.02 0.03
SD (c=−15) WD (c=−0.05)
1.00 1.00
0.75 0.75
0.50 0.50
0.25 0.25
0.00 0.00
0.00 0.01 0.02 0.03 0.00 0.01 0.02 0.03
βτ
5
Table B.2: Right-sided Test Results (%) of Q̌m−qr in the Univariate
Model
Experiment A2: First, Table B.3 compares the size performance of the proposed test
statistic Qm−qr and CCL2023 with K=8 for both the joint test H0 : βτ = 0 and the two-
sided marginal tests with H0 : βiτ = 0, for i = 1, 2, · · · , 8. Second, the power performances
of two-sided test results of Qm−qr are shown in Table B.4, in which K=8 and β1τ = β2τ =
· · · = βKτ = β̃τ are set to 0, 0.02, 0.04, 0.06, 0.08 and 0.1, respectively. Third, we
compare the size performances of L2016, CCL2023 and Qm−qr for the joint test H0 : βτ = 0
with K = 2, 3, 4, 5, 6, 7, 8 in Table B.5. Fourth, the size and the power performances of
Q̌m−qr for right-sided test in Table B.6, in which K=8 and β1τ = β2τ = · · · = βKτ = β̃τ
19
are set to 0, 0.02, 0.04, 0.06, 0.08 and 0.1, respectively. In these tables, we set τ =
19
We do not show the power performance Fan and Lee (2019) in experiment 2 since it does not offer the
details of key procedure to construct empirical CDF in multivariate models. Also, We do not show the
power performance Liu et al. (2023) in experiment 2 since the stationarity of predictors, i.e. the value of
α, must be known in its procedure but we skip the test of α here.
6
0.05, 0.25, 0.5, 0.75, 0.95 and (ρ1 , ρ2 , · · · , ρK )⊤ = (ρ)K , where
And (γ1 , γ2 , · · · , γK )⊤ = (Γ)K , Γ = (−3, 2, 1, 3, 1, −0.833, 0.667, 0.5)⊤ and thus the contem-
2 −1/2 ⊤
poraneous correlation coefficient γ1 (1 + γ12 )−1/2 , · · · , γK (1 + γK
) between ut and vt
are (Γ̃)K and Γ̃ = (−0.949, 0.894, 0.707, 0.949, 0.707, 0.64, 0.555, 0.447)⊤ . Several findings
can be observed from Tables B.3, B.4, B.5 and B.6. First, for two-sided tests H0 : βτ = 0
and H0 : βiτ = 0, i = 1, 2, · · · , 8, with K = 8, the proposed test statistic Qm−qr is free
of size distortion, unlike the CCL2023, which exhibits size distortions at quantile levels
τ = 0.05, 0.25, 0.5, 0.75, 0.95. Second, the power of the proposed test statistic Qm−qr is
satisfactory for two-sided tests H0 : βτ = 0 and H0 : βiτ = 0, i = 1, 2, · · · , 8, with K = 8.
Third, at quatile levels τ = 0.25, 0.5, 0.75, the size performance of L2016 and CCL2023 is
acceptable when the number of predictors is less than 4 and 5, respectively, for joint test
H0 : βτ = 0. However, L2016 and CCL2023 suffer size distortion for other cases, such
as the case when the number of predictors is less than 4 and 5 respectively and at tail
quantiles τ = 0.05, 0.95. This issue becomes more pronounced as the number of predictors
increases. Fourth, the proposed test statistic Qm−qr outperforms L2016 and CCL2023 for
joint test H0 : βτ = 0, since it is free of size distortion for all cases. Fifth, the proposed
test statistic Q̌m−qr perform well in terms of size and power for one-sided test, which is
not considered in literature. To sum up, the size and power of the proposed test are com-
parable to those found in the existing literature for univariate models, and it significantly
outperforms previous works in multivariate contexts.
Panel A: CCL2023
τ H0 : βτ = 0 i=1 i=2 i=3 i=4 i=5 i=6 i=7 i=8
0.05 26.2 12.8 11.9 12.2 12.2 11.6 12.1 11.6 11.8
0.25 9.5 7.2 7.4 7.3 7.5 7.5 7.6 7.9 7.0
0.5 8.6 7.6 7.5 8.2 7.4 8.0 7.5 7.4 7.6
0.75 9.5 7.9 8.5 8.4 8.0 8.1 8.7 8.0 7.9
0.95 26.6 13.0 12.0 12.2 13.2 12.2 11.9 12.4 12.3
Panel B: Qm−qr
Continued on next page
7
Table B.3 – continued from previous page
τ H0 : βτ = 0 i=1 i=2 i=3 i=4 i=5 i=6 i=7 i=8
τ H0 : βτ = 0 i=1 i=2 i=3 i=4 i=5 i=6 i=7 i=8
0.05 4.4 4.8 4.8 5.2 5.1 4.7 4.6 4.8 4.5
0.25 4.2 4.9 5.2 4.6 6.0 5.5 4.8 4.5 4.7
0.5 5.0 5.7 5.3 5.5 6.4 5.4 4.9 4.9 4.9
0.75 4.4 5.4 5.0 5.5 6.0 5.0 4.9 4.5 4.7
0.95 4.8 5.3 4.9 5.0 5.6 4.9 4.7 4.4 4.8
Note: For different i, the null hypothesis for marginal test is H0 : βiτ = 0.
Table B.4: Power and Size Performance of Two-sided Test Results (%)
of Qm−qr in Multivariate Model with K=8.
8
Table B.4 – continued from previous page
Null Hypothesis β̃τ 0 0.02 0.04 0.06 0.08 0.1
τ =0.95 4.7 9.6 26.3 51.3 76.2 90.6
τ =0.05 4.8 7.3 17.1 32.8 52.1 70.6
τ =0.25 4.5 11.3 33.5 61.9 85.1 96.0
H0 : β7τ = 0 τ =0.5 4.9 12.8 37.7 68.7 90.3 98.0
τ =0.75 4.5 11.1 32.4 61.7 85.1 95.8
τ =0.95 4.4 7.8 17.0 32.4 51.5 70.5
τ =0.05 4.5 8.1 18.1 35.5 56.9 76.0
τ =0.25 4.7 12.3 36.9 68.2 89.5 97.7
H0 : β8τ = 0 τ =0.5 4.9 13.9 42.7 74.7 93.1 98.9
τ =0.75 4.7 12.0 36.8 68.8 89.5 97.6
τ =0.95 4.8 8.3 17.5 36.0 56.2 76.8
β1τ = β2τ = · · · = βKτ = β̃τ .
Table B.5: Size Performance of L2016, CCL2023 and Qm−qr for Joint
Test H0 : βτ = 0
Panel A: L2016
τ K=2 K=3 K=4 K=5 K=6 K=7 K=8
0.05 8.0 10.6 18.6 22.6 25.6 29.4 34.2
0.25 5.3 5.9 13.0 13.2 15.3 15.1 16.8
0.5 4.6 5.3 12.2 13.1 14.0 14.0 14.6
0.75 4.8 6.2 12.7 13.5 15.3 15.7 15.1
0.95 8.4 10.7 18.5 22.9 26.1 28.7 33.1
Panel B: CCL2023
τ K=2 K=3 K=4 K=5 K=6 K=7 K=8
0.05 8.1 11.3 11.8 15.8 18.6 21.3 25.7
0.25 5.0 7.2 6.8 6.8 8.1 8.1 10.2
0.5 4.3 6.2 5.2 7.0 7.4 7.4 8.8
0.75 5.2 6.9 4.4 8.5 7.0 9.9 9.9
0.95 8.9 10.2 12.2 16.1 19.6 23.5 26.9
Panel C: Qm−qr
τ K=2 K=3 K=4 K=5 K=6 K=7 K=8
0.05 4.3 4.6 4.6 5.4 5.2 4.9 4.6
0.25 4.7 5.1 5.3 5.2 5.0 4.0 5.2
0.5 4.9 4.8 5.1 5.5 5.3 5.0 5.1
0.75 4.7 4.8 4.4 4.5 4.5 5.0 4.9
0.95 4.8 5.0 4.8 4.5 4.2 4.8 4.7
9
Table B.6: Right-sided Test Results (%) of Q̌m−qr in Multivariate Model
with K=8.
10
C Additional Theoretical Results
We eliminate the first one of two higher-order terms by applying a new instrumental
variable which is constructed by the sample splitting method and the instrumental variable
zt−1 . The size distortion induced by the higher-order terms CT in Proposition 3 could be
eliminated by removing the term T1 Tt=1 zt−1 in the test statistics. We split the full sample
P
The fitted value of equation (C.1) is x̃at−1 = µ̂ax + θ̂a zt−1 , and the residual is ṽt−1
a
= xt−1 − x̃at−1
and t = 1, 2, · · · , T0 . The second step based on the first sub-sample is as follows.
h i⊤ T0
X
a ⊤ a ⊤
a
ρτ yt − µτ − βτ⊤ x̃at−1 − γτ⊤ ṽt−1
a
µ̂τ , (β̂τ ) , (γ̂τ ) = arg min . (C.2)
µτ ,βτ ,γτ
t=1
The fitted value of equation (C.3) is x̃bt−1 = µ̂bx + θ̂b zt−1 , and the residual is ṽt−1
b
= xt−1 − x̃bt−1
and t = 1, 2, · · · , T0 . The second step based on the second sub-sample is as follows.
h i⊤ T
X
µ̂bτ , (β̂τb )⊤ , (γ̂τb )⊤ ρτ yt − µτ − βτ⊤ x̃bt−1 − γτ⊤ ṽt−1
b
= arg min . (C.4)
µτ ,βτ ,γτ
t=T0 +1
11
a 1
PT0 b 1
PT
where z̄t−1 = zt−1 − T0 t=1 zt−1 and z̄t−1 = zt−1 − T −T0 t=T0 +1 zt−1 . Second, we apply
equations (1), (2), (C.5) and (C.6) to obtain the following equations.
T
1
P
T T T
zt−1 X
T
X 1 X t=1
z̄t−1 x⊤
t−1 (β̂τ − β) = zt−1 ψτ (utτ ) − ψτ (utτ ) + op (DT−1 ), (C.7)
t=1
fuτ (0) t=1
fuτ (0) t=1
T0
1
P
T0 T0 T0
zt−1 X
T0
X 1 X t=1
a
z̄t−1 x⊤ a
t−1 (β̂τ − β) = zt−1 ψτ (utτ ) − ψτ (utτ ) + op (DT−1 ), (C.8)
t=1
fuτ (0) t=1
fuτ (0) t=1
T T
P
P zt−1
T zt−1 ψτ (utτ ) t=T0 +1 T
t=T0 +1
X T −T0
X
b
z̄t−1 x⊤ b
t−1 (β̂τ − β) = − ψτ (utτ ) + op (DT−1 ). (C.9)
t=T0 +1
fuτ (0) fuτ (0) t=T0 +1
(C.11)
By subtracting the sum of equations (C.10) and (C.11) from equation (C.7), the term
PT
t=1 ψτ (utτ ) of β̂τ is removed as follows.
PT ⊤
PT0 a ⊤
PT b ⊤
where W1 = t=1 z̄t−1 xt−1 , W2 = Sa t=1 z̄t−1 xt−1 , W3 = Sb t=T0 +1 z̄t−1 xt−1 . Third,
we could utilize another instrumental variable estimator with the IV z̃t−1 . Define the IV
estimator β̂τl0 as follows.
12
T
P
By equations (C.12) and (C.13) and the fact W1 − W2 − W3 = z̃t−1 xt−1 , it follows that
t=1
T
!−1 T
1 X X
β̂τl0 − β = z̃t−1 x⊤
t−1 z̃t−1 ψτ (utτ ) + op (DT−1 ). (C.14)
fuτ (0) t=1 t=1
PT
The key and desirable property of the new instrumental variable z̃t−1 is t=1 z̃t−1 = 0.
Thus it follows that
β̂τl0 − β
" T T
! #−1 T T
!
1 X 1X X 1X
= z̃t−1 − z̃t−1 xt−1⊤
z̃t−1 − z̃t−1 ψτ (utτ ) + op (DT−1 )
fuτ (0) t=1 T t=1 t=1
T t=1
T
! −1 T
1 X X
= ⊤
z̃t−1 xt−1 z̃t−1 ψτ (utτ ) + op (DT−1 ).
fuτ (0) t=1 t=1
Although the higher-order term CT disappear in Q̌l−qr and Ql−qr , the higher-order term
BTl which arise by the same reason of BT still exists. To analyze these distortion effects
of Q̌l−qr and Ql−qr intuitively, we show their higher-order terms with J = K = 1 here.
Following the Proposition 2 of Liao et al. (2024), we have the proposition as follows.
Proposition C.1. Under Assumptions 1 and 2 and the null hypothesis H0 : βτ = 0 and
J = K = 1, for SD predictors, the following equation holds as T → ∞ that
PT d
where ZTl = (Σzz )−1/2 → N(0, 1) with BTl → 0 and
t=1 z̃t−1 ψτ (utτ )
l l 1n XT
∗ ∗ ⊤
o
BT = ϖl ZT , ϖl = − τ (1 − τ ) z̃t−1 (z̃t−1 ) − 1 (C.15)
2 t=1
∗
where z̃t−1 = (Σzz )−1/2 z̃t−1 and
√
T (1−δ)/2 BTl = RTl − Wl ρvψ / −2cz (C.16)
√ √
and RTl = Wa R1,T l l
+ Wb R2,T , E R1,T l
= E R2,Tl
= 0 and Wl = Wa / λ + Wb / 1 − λ and
Wa = [T −(1+δ) Tt=1 z̃t−1 z̃t−1 ⊤
û2t ]−1/2 (1 − Sa )[T −(1+δ) Tt=1 ⊤
û2t ]1/2 and Wb = [T −(1+δ)
P P 0
zt−1 zt−1
PT ⊤ 2 −1/2
(1 − Sb )[T −(1+δ) Tt=T0 +1 zt−1 zt−1⊤
û2t ]1/2 . More details about R1,T l
P
t=1 z̃t−1 z̃t−1 ût ] and
l
R2,T are discussed in Proposition 2 of Liao et al. (2024).
13
Remark 3. As specified by Liao et al. (2024), the term RTl could be regarded as the “residual
√
term” of BTl and thus only −Wl ρvψ / −2cz is the main source of size distortion by BTl for
Q̌l−qr . So it is evident that the bigger |cz | = −cz is, the smaller the size distortion by BTl
is.
14
CP
4
2
0
−2
−4
LN1
4
2
0
−2
Transformed Test Statistics
−4
LN2
4
2
0
−2
−4
TDH
4
2
0
−2
−4
TI
4
2
0
−2
−4
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
τj
In Figure D.1, the test statistic Q̌tm−qr is significantly positive at the 0.01 level in the
light green area, while it is significantly negative at level 0.01 in the gray area. These
results are consistent with those of Qm−qr and the economic explanations in subsection 4.2
of the main text.
15
D.2 Tail Risk Indicator for U.S. Treasury Bonds
ˆτj (T̃ ), β̇ˆτj (T̃ )⊤ ⊤ = arg minµτ ,βτ
PT̃ −1
to obtain the real time coefficient estimator µ̇ j j t=1 ρτj (yt −
µτj − β̇τ⊤j ẋjt−1 ), where ẋjt−1 is the one used in (29). Then the out-of-sample predicted value
of bond risk premia at quantiles τ is ŷ o (T̃ ) = µ̇
j
ˆ (T̃ ) + β̇ˆ (T̃ )⊤ ẋj and the out-of-sample
τj τj τj t−1
quantile regression error is ûoτj (T̃ ) = yT̃ − ŷτoj (T̃ ). The out-of-sample performance in the
period from T0 to T is evaluated by qw-CRPS. Using qw-CRPS and the out-of-sample pre-
diction results, we construct the tail-risk indicator by the following procedure. The basic
idea to construct tail risk indicator is as follows. Since W (τj ) > 0, then qw-CRPS could
be written as
J−1 J−1
2 X 2 X
qw − CRPST̃ = W (τj ) ρτj [ûoτj (T̃ )] = ρτj [ûw
τj (T̃ )], (D.1)
J − 1 j=1 J − 1 j=1
where ûw o o o w w
τj (T̃ ) = W (τj ) ûτj (T̃ ). Moreover, since yT̃ = ŷτj (T̃ )+ ûτj (T̃ ), we have yT̃ = ŷτj (T̃ )+
ûw w w o
τj (T̃ ), where j = 1, 2, · · · J, yT̃ = W (τj ) yT̃ and ŷτj (T̃ ) = W (τj ) ŷτj (T̃ ). As per equation
(D.1), qw − CRPST̃ could be viewed as the quantile loss function for the aforementioned
weighted quantile regression at quantile τj , where j = 1, 2, · · · , J. Therefore, it is reasonable
to use the out-of-sample prediction value of weighted quantile regression to construct the
left and the right tail risk indicators for U.S. treasury bonds as
J
X J
X
R(tail)T̃ = ŷτwj (T̃ ) = o
W (τj ) ŷtτj
, T̃ = Tm , Tm + 1, · · · , T.
t=1 t=1
To model the right tail risk indicator, we set W (τj ) = Wr ≡ (τj ) = τj2 (τ12 + τ22 + · · · + τJ2 )−1
and W (τj ) = Wl (τj ) = −(1 − τj )2 [(1 − τ1 )2 + (1 − τ2 )2 + · · · + (1 − τJ )2 ]−1 to model the
20
left tail risk indicator. Therefore, the larger the values of the left and right tail risk
indicators, R(tail)t , the greater the risk associated with 2- and 3-year bonds and 4- and
5-year bonds, respectively.
20
Wr (τj ) and Wl (τj ) are the unitized aforementioned weight τj2 evaluating the performance at right tail
quantiles and (1 − τj )2 evaluating the performance at left tail quantiles respectively such that Wr (τ1 ) +
Wr (τ2 ) + · · · + Wr (τJ ) = 1 and Wl (τ1 ) + Wl (τ2 ) + · · · + Wl (τJ ) = −1.
16
Left Tail Risk Indicator of 2−year Bond
0.00
−0.25
−0.50
−0.75
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
0.0
Tail Risk Indicator R(tail)t
−0.5
−1.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
0.0
−0.5
−1.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
0.5
0.0
−0.5
−1.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
Time
17
Right Tail Risk Indicator of 2−year Bond
0.4
0.2
0.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
0.2
0.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
0.8
0.4
0.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
0.8
0.4
0.0
1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024
Time
In the following part, the left and the right tail indicator constructed based on Qm−qr
and CCL2023 in 1992–2022 is shown in Figures D.2 and D.3. In Figure D.2, the left tail
risk indicator successfully highlights two peaks during the 2008 global financial crisis and
the COVID-19 pandemic, indicating increased risk for 2- and 3-year bonds during these
periods. Meanwhile, in Figure D.3, the right tail risk indicator successfully highlights four
peaks during the periods of 1990-1992 (third oil crisis), 2002-2003 (dotcom bubble crisis),
and twice for the 2008 global financial crisis, indicating elevated risk for 4- and 5-year bonds
18
during these times. Although the left and right tail risk indicators based on CCL2023 have
a similar shape to those based on Qm−qr and also capture these peaks, our approach is
ahead of CCL2023. As a result, our tail risk indicators are capable of detecting bond risk
more accurately and at an earlier stage than those derived from CCL2023. This superiority
is attributed to several key factors: the test conducted in CCL2023 (and thus its tail risk
indicators) did not find significant predictive power for TDH across all bond maturities at
both tails quantiles, nor for TI at the right tail quantiles of 4- and 5-year bonds. This
is in contrast to the findings associated with our Qm−qr test, which identified significant
predictive power of those predictors under the same conditions.
E Selected Proof
Theorems with WD predictors are easy to show using conventional technics, so we focus
on the proof with SD predictors in this section. Equation numbers without the appendix
⊤
section prefix refer to those in the main text. Define Xt−1 = T −1/2 , DT−1 x̃⊤
t−1 , ḊT
−1 ⊤
ṽt−1 .
Lemma 1. Under Assumption 1 and 2, the following equation holds for SD predictors.
X T Z 1
⊤ ⊤ −1
Xt−1 Xt−1 ⇒ diag[1, Ωzx Ωzz Ωzx , J¯xc (r)J¯xc (r)⊤ dr]
t=1 0
By equation (5) and (E.1) and continuous mapping theorem, it follows that
T T
!−1
X X
θ̂ = DT−2 ⊤
xt−1 z̄t−1 DT−2 ⊤
z̄t−1 z̄t−1 ⇒ Ω⊤ −1
zx Ωzz . (E.2)
t=1 t=1
By equation (E.1) and (E.2), definition of x̃t−1 and continuous mapping theorem,
T
X T
X
T −1/2
DT−1 x̃t−1 = θ̂DT−1 T −1/2 zt−1 = Op (1)Op [T (δ−1)/2 ] = op (1). (E.3)
t=1 t=1
19
By the equation xt−1 = µ̂x + θ̂zt−1 + ṽt−1 and equation (E.4), we have
T T
1X 1X
µ̂x = xt−1 − θ̂ zt−1 . (E.5)
T t=1 T t=1
√
By equations (E.2) and (E.5) and the fact that x⌊rT ⌋ / T ⇒ Jxc (r) for 0 ≤ r ≤ 1, it follows
that
T Z 1
−1/2 1 X
T µ̂x = xt−1 + op (1) ⇒ Jxc (r)dr. (E.6)
T 3/2 t=1 0
By equation (E.1), (E.2), (E.3) and (E.6), definition of x̃t−1 and continuous mapping the-
orem,
T
X T
X
DT−1 x̃t−1 x̃⊤ −1
t−1 DT = θ̂DT−1 ⊤
zt−1 zt−1 DT−1 θ̂⊤ + op (1) ⇒ Ω⊤ −1
zx Ωzz Ωzx . (E.7)
t=1 t=1
As per the equation xt−1 = µ̂x + θ̂zt−1 + ṽt−1 and equation (E.6), it becomes evident that
T −1/2 ṽ⌊rT ⌋ = T −1/2 x⌊rT ⌋ − T −1/2 µ̂x + op (1) ⇒ J¯xc (r), 0 ≤ r ≤ 1. (E.8)
Lemma 2. Under Assumption 1 and 2, the following equation holds for SD predictors.
T Z 1 ⊤
X ⊤
⊤ −1
J¯xc (r)⊤ dBψτ (r)
Xt−1 ψτ (utτ ) ⇒ MN 0, diag(1, Ωzx Ωzz Ωzx ) , . (E.10)
t=1 0
20
holds by Lee (2016). In light of equation (E.2) and (E.11) and the continuous mapping
theorem, it is clear that
T T
X ⊤ X ⊤
T −1/2
, DT−1 x̃⊤
t−1 ψτ (utτ ) = diag(1, θ̂) T −1/2 , DT−1 zt−1
⊤
ψτ (utτ ) (E.12)
t=1 t=1
⇒ MN 0, diag(1, τ (1 − τ )Ω⊤ −1
zx Ωzz Ωzx ) .
√
FCLT shown in equation (3) and the fact x⌊rT ⌋ / T ⇒ Jxc (r) for 0 ≤ r ≤ 1 imply the
following equation by the continuous mapping theorem.
T Z 1
1X
ṽt−1 ψτ (utτ ) ⇒ J¯xc (r)dBψτ (r). (E.13)
T t=1 0
Using equations (E.12) and (E.13) and the joint convergence shown by equation (3), we
can determine that
T Z 1 ⊤
X ⊤
⊤ −1 ¯c ⊤
Xt−1 ψτ (utτ ) ⇒ MN 0, diag(1, τ (1 − τ )Ωzx Ωzz Ωzx ) , Jx (r) dBψτ (r) .
t=1 0
(E.14)
Lemma 3 (Lemma A.1 of online appendix of CCL2023). Let VT (v) be a vector function that
satisfies: (i) −v ⊤ VT (λv) ≥ −v ⊤ VT (v) for any λ ≥ 1. (ii) sup ∥VT (v) + N v − AT ∥ = op (1)
∥v∥<M
where ∥AT ∥ = Op (1), 0 < M < ∞, and N is a positive-definite random matrix. Suppose
that vT is a vector such that ∥VT (vT )∥ = op (1), then ∥vT ∥ = Op (1) and vT = N −1 AT +op (1).
Proof of Proposition 1. Following the same procedure of the proof of Theorem 1 in online
appendix of CCL2023, Proposition 1 holds by Lemma 1, 2 and 3.
Proof of Theorem 1. From Proposition 1 and equation (E.2) and the orthogonality shown
21
in equation (E.4), we can deduce the following equation.
T
!−1 T
1 X ⊤
X
DT β̂τ − βτ = DT−2 x̃t−1 x̃t−1 DT−1 x̃t−1 ψτ (utτ ) + op (1) (E.15)
fuτ (0) t=1 t=1
" T
#−1 T
1 X X
= DT−2 θ̂z t−1 (θ̂z t−1 )⊤ DT−1 x̃t−1 ψτ (utτ ) + op (1)
fuτ (0) t=1 t=1
T
!−1 T
1 −2
X
⊤ ⊤ −1
X
= θ̂DT z t−1 z t−1 θ̂ θ̂DT z t−1 ψτ (utτ ) + op (1)
fuτ (0) t=1 t=1
T
!−1 T
1 ⊤ −1 −2
X
⊤ −1 −1
X
= (θ̂ ) DT z t−1 z t−1 θ̂ θ̂DT z t−1 ψτ (utτ ) + op (1)
fuτ (0) t=1 t=1
T
!−1 T
1 X X
= (θ̂⊤ )−1 DT−2 z t−1 z ⊤t−1 D −1
T z t−1 ψτ (utτ ) + op (1)
fuτ (0) t=1 t=1
T
!−1 T
1 X X
= DT−2 z̄t−1 x⊤t−1 D T
−1
z̄t−1 ψτ (utτ ) + op (1)
fuτ (0) t=1 t=1
Equation (E.1), (E.11) and (E.15) and Slutsky’s Theorem imply that
h i
MN 0K , 1 2 Ω−1 −1 ⊤
zx Ωzz (Ωzx ) , SD;
d fuτ (0)
DT β̂τ − βτ → − h i . (E.16)
N 0K , 1 2 Ω−1 Ωzz (Ω−1 )⊤ , W D;
fu (0) zx τ
zx
Proof of Proposition 2. Using equations (9), (10), (14) and (E.1), Theorem 1 and Slutsky’s
theorem, Proposition 2 is proved.
The proof of equation (C.5) and (C.6) closely mirrors the proof of Theorem 1, therefore,
we omit it for brevity.
−−→
Proof of Equation (13). Equation (12) and the fact that ξt−1 follows i.i.d. N(0M , IM ) imply
T T
(i)
that √1T ξj,t−1 ψτ (utτ ) FT ∼ N 0, T1 ψτ (utτ )2 and thus
P P
t=1 t=1
" T
#−1/2 T
1X 1 X (i)
ψτ (utτ )2 √ ξj,t−1 ψτ (utτ ) FT ∼ N(0, 1), (E.17)
T t=1 T t=1
22
h i⊤
(i) (i) (i) (i)
where i = 1, · · · , M1 and j = 1, · · · , M2 and ˆlτ = ˆl1,τ , ˆl2,τ , · · · , ˆlM2 ,τ . By the indepen-
(i)
dence of ξj,t−1 for j = 1, · · · , M2 and i = 1, · · · , M1 , it can be verified that conditional
h P i−1/2 PT (i)
1 T
on the whole sample FT , T t=1 ψτ (utτ ) 2 √1
T t=1 ξj,t−1 ψτ (utτ ) are uncorrelated for
all i and j. Therefore, equation (E.17) and the uncorrelation between the (i,j)th items of
h P i−1/2 PT (i)
1 T
ψ (u ) 2 √1
T t=1 τ tτ T t=1 ξj,t−1 ψτ (utτ ) lead to the following joint convergence.
" T
#−1/2 ( T T
1X 1 X (1) ⊤ 1 X (2) ⊤
ψτ (utτ )2 √ (ξt−1 ) ψτ (utτ ), √ (ξ2,t−1 ) ψτ (utτ ), · · · , (E.18)
T t=1 T t=1 T t=1
T
)
1 X (M1 ) ⊤
√ (ξt−1 ) ψτ (utτ ) FT ∼ N [0M , IM ] ,
T t=1
where M = M1 M2 . Hence,
" T
#−1 M M " T
#2
1 X 2
1X X 1 X (i)
ψτ (utτ )2 √ ξj,t−1 ψτ (utτ ) FT ∼ χ2M , (E.19)
T t=1 i=1 j=1
T t=1
[τ (1 − τ )]2 1
= 2M + op (1)
fuτ (0)4 M 2
[τ (1 − τ )]2 2
= + op (1) → 0,
fuτ (0)4 M
23
as M → ∞ and T → ∞. Similarly, equations (12) and (E.20) imply that
( M1 XM2 h
)
1 X √ (i) i2
E T ˆlj,τ FT (E.23)
M i=1 j=1
" #2
M1 X M2 T
1 1 X 1 X (i)
= E √ ξj,t−1 ψτ (utτ ) FT + op (1)
fuτ (0)2 M i=1 j=1 T t=1
" #−1 M M " #2
T 1 2 T
τ (1 − τ ) 1 1X 2
X X 1 X (i)
= E ψτ (utτ ) √ ξj,t−1 ψτ (utτ ) + op (1) FT + op (1)
fuτ (0)2 M T t=1 i=1 j=1
T t=1
" #−1 M M " #2
T 1 X 2 T
τ (1 − τ ) 1 1 X
2
X 1 X (i)
= E ψτ (u tτ ) √ ξ ψ
j,t−1 τ (utτ ) F T + op (1)
fuτ (0)2 M T t=1 i=1 j=1
T t=1
τ (1 − τ ) 1 τ (1 − τ )
= 2
M + op (1) = + op (1),
fuτ (0) M fuτ (0)2
Equations (E.22) and (E.23) and the continuous mapping theorem imply that
M1 XM2 h
1 X √ (i) i2 P τ (1 − τ ) 1 1
T ˆlj,τ FT −
→ 2
M= τ (1 − τ ), (E.24)
M i=1 j=1 fuτ (0) M fuτ (0)2
as M → ∞ and T → ∞.
24
Proof of Lemma 4, Theorem 2. The proofs are very similar to the Theorem 2 of Liao et al.
(2024), and thus are omitted for brevity.
Proof of Proposition 4. The proof closely follows the proof of Theorem 3 of Liao et al.
(2024), and is therefore omitted.
Proof of Proposition C.1. The proof directly follows from the proof of Proposition 2 in Liao
et al. (2024).
R(β̂τl − βτ ) Rβτ − rτ
=h i1/2 + h i1/2
[ β̂τl )R⊤
R Avar( [ β̂τl )R⊤
R Avar(
As per Lemma 5 and Theorem 2 and the continuous mapping theorem, it become evident
that
Rβ̂τl − βτ d
h i1/2 →
− N(0, 1). (E.30)
[ β̂τl )R⊤
R Avar(
Rβτ − rτ h
l ⊤
i−1/2
h i1/2 ⇒ R Avar(β̂τ )R bτ (E.31)
[ β̂τl )R⊤
R Avar(
Using equations (E.30) and (E.31) and the continuous mapping theorem, we have the
following result that
h i−1/2
d d
− Q̌∗l−qr = N(0, 1) + R Avar(β̂τl )R⊤
Q̌l−qr → bτ (E.32)
25
when J=1. Moreover, by equation (19), it follows that
⊤ n o−1
Ql−qr = Rβ̂τl − rτ [ β̂τl )R⊤
R Avar( Rβ̂τl − rτ (E.33)
⊤ n o−1
= Rβ̂τl − Rβτ + Rβτ − rτ [ β̂τl )R⊤
R Avar( Rβ̂τl − Rβτ + Rβτ − rτ
⊤ n o−1
= Rβ̂τl − Rβτ [ β̂τl )R⊤
R Avar( Rβ̂τl − Rβτ
⊤ n o−1
+ Rβ̂τl − Rβτ [ β̂τl )R⊤
R Avar( (Rβτ − rτ )
n o−1
+ (Rβτ − rτ )⊤ R Avar(
[ β̂τl )R⊤ Rβ̂τl − Rβτ
n o−1
+ (Rβτ − rτ )⊤ R Avar(
[ β̂τl )R⊤ (Rβτ − rτ )
As per Lemma 5 and Theorem 2 and the continuous mapping theorem, it becomes evident
that
⊤ n o−1
Rβ̂τl − Rβτ [ β̂τl )R⊤
R Avar( Rβ̂τl − Rβτ ⇒ χ2J , (E.34)
and
⊤ n o−1 h i−1
d
Rβ̂τl − Rβτ [ β̂τl )R⊤
R Avar( − ZQ = Zβ⊤ R Avar(β̂τl )R⊤
(Rβτ − rτ ) → bτ , (E.35)
and
n o−1
d
(Rβτ − rτ )⊤ R Avar(
[ β̂τl )R⊤ Rβ̂τl − Rβτ →− ZQ⊤ = ZQ . (E.36)
By equations (E.33), (E.34), (E.35), (E.36) and (E.37) and Slutsky’s Theorem, Proposition
5 holds.
The proof of equation (24) is very similar to the proof of Proposition 5, so we omit it.
Proof of Theorem 3. By Proposition 5 and equations (23) and (24) and the continuous
mapping theorem, Theorem 3 holds.
Proof of Theorem 4. By Proposition 5 and equations (26) and (27) and the continuous
mapping theorem, Theorem 4 holds.
26