Excess Returns : Empirical Investigations

We review the classical work of Campbell and Shiller [1] on excess returns from trading long-maturity bonds and funding it with a short-maturity bond. Subsequently we also look at the modern developments from the work of Cochrane and Piazzesi [2].

Definitions

We define general expression for excess returns that involves taking long and short positions in two bonds of different maturities.

Let $P_t^k$ be the price of the $k$-maturity discount bound at time $t$.

At time $t=0$, we buy a $\frac{1}{P_0^N}$ units of $N$-maturity bond for \$1. To fund the \$1, we short $\frac{1}{P_0^n}$ units of $n$-maturity bond to get \$1.

So our portfolio at time $t=0$ is, $$ \mathbb{\Pi}_0(\tau, N, n) = \left(\frac{1}{P_0^N}\right)P_0^N -\left(\frac{1}{P_0^n}\right)P_0^n = 0 $$ where the quantities in the brackets are the notionals.

We then wait for a period of time $\tau$ to elapse, and at which time we unwind the positions. We have to sell the $N$-maturity bond we bought at time $t=0$ for price $P_{\tau}^{N-\tau}$, and buy back the short position for $P_{\tau}^{n-\tau}$. The time-$\tau$ portfolio will be worth,

$$ \mathbb{\Pi}_\tau(\tau, N, n) = \left(\frac{1}{P_0^N}\right)P_\tau^{N-\tau} -\left(\frac{1}{P_0^n}\right)P_\tau^{n-\tau} $$

We therefore define the quantity $\mathbb{\Pi}_\tau-\mathbb{\Pi}_0 = \mathbb{\Pi}_\tau$ to be the $\tau$-period excess return from investing in the $N$-period bond and financing it with the $n$-period bond,

$$ \textit{xret}(\tau, N, n) =\mathbb{\Pi}_\tau (\tau, N, n) $$

Case $\tau=n=1$

In this case we short a one year bond, use it to finance a $N$-maturity bond, and unwind it after 1 year. The excess return is,

$$ \begin{aligned} \textit{xret}(1, N, 1) &= \left(\frac{1}{P_0^N}\right)P_1^{N-1} -\left(\frac{1}{P_0^n}\right)P_1^0 = \left(\frac{1}{P_0^N}\right)P_1^{N-1} -\left(\frac{1}{P_0^n}\right)\\ &=\textit{xret}_{t+1} \end{aligned} $$

We know that by definition, $P_t^T = e^{-y_t^T (T-t)} \approx 1 - y_t^T(T-t) + \ldots$, using this, we have,

$$ \begin{aligned} \textit{xret}_{t+1} &= \left(\frac{1}{P_0^N}\right)P_1^{N-1} -\left(\frac{1}{P_0^n}\right)\\ &=(1+ y_0^N(N-0))(1 - y_1^{N-1} ( N-1)) - (1+y_0^1 ( 1-0))\\ &=(1+ y_0^NN)(1 - y_1^{N-1} ( N-1)) - (1+y_0^1)\\ &=y_0^NN - y_1^{N-1} ( N-1) -y_0^1 \qquad \text{(retaining 1st order terms)} \end{aligned} $$
In [2]:
import datetime
import quandl
import pandas as pd
import matplotlib.pyplot as plt
from   mpl_toolkits.mplot3d import Axes3D
from   matplotlib import cm
from   matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
from   numpy.linalg import eig
from   scipy.stats import linregress
from   statsmodels.formula.api import ols
import warnings
from   IPython.core.display import HTML, display

warnings.filterwarnings('ignore')
%matplotlib inline
plt.style.use(('bmh'))
In [3]:
# data = quandl.get("FED/SVENY") # get from quandl if we don't have a local copy 
dataLoc = './FED-SVENY.csv'
data = pd.DataFrame.from_csv(dataLoc)
data.sort_index(inplace=True)
# find the index onwards from which we have data for all tenors
idx_mat = np.isnan(data).apply( np.nonzero, axis=0).apply(np.ravel)
first_idx = np.max( [ np.max(a) for a in idx_mat.values if len(a) > 0 ])
data_subset = data[first_idx+1:]
data_subset = data_subset.applymap( lambda x : x * 0.01 ) # convert from percentage to decimals
data_subset = data_subset.resample('M').copy() # down sample to monthly frequency
In [16]:
fig = plt.figure( figsize = ( 12, 6))
for n, col in enumerate( data_subset.columns ):
    if n % 3 == 0:
        # plot last 200 points
        plt.plot( data_subset.index[-200:], data_subset[col][-200:], lw=2. )

We calculate the excess returns from long a 10-yr bond while funding it with a 1-yr bond, then unwinding the positions after 1yr.

In [17]:
# compute the excess returns. first we truncate the data and drop the last 1-yr (365days)
# since we won't be able to calculate returns from them.
data_1_10 = data_subset[['SVENY01','SVENY09','SVENY10']]

def excess_returns( data_set, N ):
    xrets = []
    time_points = []
    
    for ts, series in data_set.iterrows():
    
        #find the yield of the 9-yr bond, one yr from current date
        one_yr = ts.date() + datetime.timedelta(days=365)
        idx = data_set.index.searchsorted( one_yr )
    
        if idx == len( data_set ):
            break
        
        if data_set.index[idx] == one_yr:
            next_yr_bond = data_set.SVENY09.iloc[idx]
        else:
            next_yr_bond = data_set.SVENY09.iloc[idx-1]
        
        xrets.append( series.SVENY10 * N - next_yr_bond * (N-1) - series.SVENY01 )
        time_points.append( ts.date() )
        
    return ( np.array(xrets), time_points )

ret, xs = excess_returns( data_1_10, 10 )
In [19]:
fig = plt.figure( figsize = ( 12, 6 ) )
plt.plot( xs, ret, lw='2.0')
plt.title('Excess returns from long 10-y and short 1-yr bond strategy. Sharpe ratio = %s'
           % (ret.mean()/ret.std(), ) )
plt.show()

Slope as a Predicting Factor

We now test the claim that the slope of the yield curve predicts excess returns. In the work of Campbell and Shiller (1991) [1], they regressed the differences between the time-$(t+1)$, $(n-1)$-maturity yield and the time-$t$, $n$-maturity yield (the right-hand variable) against the time-$t$ differences between the $n$-maturity yield and the funding rate, $y_t^1$, (the right-hand variables):

$$ y_{t+1}^{n-1} - y_t^n= \alpha_n + \beta_n \frac{y_t^n - y_t^1}{n-1}\qquad \text{(*)} $$

What's the intuition behind the above?

Recall that previously we defined the excess return as,

$$ \begin{aligned} \textit{xret}_{t+1}&=y_0^NN - y_1^{N-1} ( N-1) -y_0^1 \\ &= \underbrace{(y_0^N - y_1^{N-1})}_{\text{roll-down term}} ( N-1) + \underbrace{y_0^N-y_0^1}_{\text{slope or carry}} \end{aligned} $$

The above expression after rearrangement looks very similar to that in $(*)$. Suppose if we regress the carry ($\frac{y_0^N-y_0^1}{N-1}$) against the roll-down $y_0^N - y_1^{N-1}$ and if we find a relationship between the two, then we can use it to predict the excess return at any point in time.

As long as the yield curve is upward sloping, we get a boost from the positive carry term. The other component of the excess returns, the roll-down, it dependent on the behavior of the $(n − 1)$-year yield at time-$(t + 1)$ : $y_{t+1}^{n-1}$.

We shall see that indeed there is a relationship which is of statistical significance.

In [20]:
LHS, RHS = [], []
N=10
for ts, series in data_1_10.iterrows():
    
    #find the yield of the 9-yr bond, one yr from current date
    one_yr = ts.date() + datetime.timedelta(days=365)
    idx = data_1_10.index.searchsorted( one_yr )
    
    if idx == len( data_1_10):
        break
        
    if data_1_10.index[idx] == one_yr:
        next_yr_bond = data_1_10.SVENY09.iloc[idx]
    else:
        next_yr_bond = data_1_10.SVENY09.iloc[idx-1]
    
    LHS.append( next_yr_bond - series.SVENY10 )
    RHS.append( ( series.SVENY10 - series.SVENY01 ) / ( N - 1 ) )
    
fig = plt.figure( figsize = ( 12, 6 ) )

plt.scatter(RHS, LHS, color='#EC7063', edgecolors='black')
plt.show()
In [21]:
res = linregress( RHS, LHS )
regression_data = pd.DataFrame( {'y':RHS, 'x':LHS} )

model = ols( 'y ~ x', regression_data ).fit()
model.summary()
Out[21]:
OLS Regression Results
Dep. Variable: y R-squared: 0.034
Model: OLS Adj. R-squared: 0.031
Method: Least Squares F-statistic: 12.56
Date: Sat, 22 Oct 2016 Prob (F-statistic): 0.000447
Time: 12:05:00 Log-Likelihood: 1895.7
No. Observations: 360 AIC: -3787.
Df Residuals: 358 BIC: -3780.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [95.0% Conf. Int.]
Intercept 0.0017 7.08e-05 24.671 0.000 0.002 0.002
x -0.0294 0.008 -3.543 0.000 -0.046 -0.013
Omnibus: 158.509 Durbin-Watson: 0.028
Prob(Omnibus): 0.000 Jarque-Bera (JB): 21.000
Skew: -0.125 Prob(JB): 2.75e-05
Kurtosis: 1.843 Cond. No. 125.

As we can see from the above, the relationship is statistically significant.

Excess Returns in terms of Forward Rates

We now explore the formulation of excess return in terms of forward rates.

Definitions

We work in time intervals of 1-yr and consider the construction of the following portfolio.

At time $t=0$,

  • Buy a $N$-maturity bond, paying $-P^N$ dollars
  • Sell $\frac{P^{N}}{P^{N+1}} \times P^{N+1}$ for $P^N$ dollars

So the initila value of our portfolio is $0$.

At time $t=N$,

  • We get $\$1$ from our $N$-maturity bond.

At time $t=N+1$,

  • We pay $\$1$ dollar for every unit notional of the $P^{N+1}$ that we sold. Since we sold $\frac{P^{N}}{P^{N+1}}$ units, we pay $\frac{P^{N}}{P^{N+1}} > 1$ at time $t=N+1$.

So what we have done is that we have synthesized a contract today for borrowing $\$1$ at $t=N$ and paying back at $t=N+1$. This is essentially, a forward rate.

Let us now define log-prices and forward rates,

  • log-price of a bond : $p_t^T=\log P_t^T$
  • the forward rate : we saw earlier that we borrowed $\$1$ at time $t=N$ and had to pay back $\frac{P^{N}}{P^{N+1}}$ at $t=N+1$. Assuming continuous compounding over 1 period, it means that the forward rate $f_t^N$ satisfies $e^{f_t^N \times1}=\frac{P^{N}}{P^{N+1}}$, which gives the forward rate : $f_t^N = p_t^N - p_t^{N+1}$.

By convention, we use the forward rate definition whereby we borrow at time-$t+n-1$ and pay back at time-$t+n$,

$$ f_t^n = p_t^{t+n-1} - p_t^{t+n} $$

Assume $p_t^{t+0}=0$, then we can write the following,

$$ \begin{aligned} f_t^1 &= p_t^{t+0} - p_t^{t+1}\\ &= 0 - p_t^{t+1}\\ f_t^2 &= p_t^{t+1} - p_t^{t+2}\\ \ldots\\ f_t^n &= p_t^{t+n-1} - p_t^{t+n} \end{aligned} $$

Clearly, from the above we see that we have a telescopic sum. Adding up $f_t^1, f_t^2, \ldots f_t^n$, we get,

$$ p_t^{t+n} = -\left(\sum_{i=1}^n f_t^i \right) $$

Given this, we can express the excess returns formula as,

$$ \begin{aligned} \textit{xret}_{t+1} &= y_t^nn - y_{t+1}^{n-1} ( n-1) -y_t^1 \\ &=p_{t+1}^{n-1} - p_t^n -y_t^1\\ &=\sum_{i=1}^{n} f_t^i - \sum_{i=1}^{n-1} f_{t+1 }^i + p_t^1\\ &=\sum_{i=2}^{n} f_t^i - \sum_{i=2}^{n} f_{t+1 }^{i-1} \qquad \text{(adjusting the index in 2nd sum and subtract $p_t^1$ from the 1st)} \\ &=\sum_{i=2}^{n} \left(f_t^i - f_{t+1 }^{i-1} \right) \end{aligned} $$

This expression shows that the excess returns reaped by engaging in the 'carry' trade with the $n$-maturity bond can be expressed as the sum of $n-1$ terms. Each term is the difference between today's forward rate for maturity $i$ $\left(f_t^i\right)$ and the forward rate in one year's time with maturity decreased by 1 yr $\left(f_{t+1}^{i-1}\right)$.

Cochrane and Piazzesi - Forward Rates as Returning Predicting Factors

In their work in [2], Cochrane and Piazzesi found that five forward rates predict excess returns with $R^2$ as high as $44\%$.

They define,

  • excess return for holding an $n$-maturity bond : $rx_{t+1}^{(n)} = p_{t+1}^{(n-1)} - p_t^{(n)} - y_t^{(1)} $.
  • the forward rate : $f_t^{(n)}= p_t^{(n-1)}-p_t^{(n)}$

which we derived above.

They then regressed the excess return against the five forward rates : $f_t^{(1)}=y_t^{(1)}$, $f_t^{(2)}$, $f_t^{(3)}$, $f_t^{(4)}$, $f_t^{(5)}$.

We validate their claims below using a Python implementation and investigate how the excess returns look like as well as the shape of the return predicting factors.

Note: In their analysis, Cochrane and Piazzesi used CRSP data from Jan-1964 to Dec-2003.

In [22]:
st = datetime.date(1964, 1, 1)
et = datetime.date(2003, 12, 31)
CP_data = pd.DataFrame.from_csv(dataLoc)
CP_data.sort_index(inplace=True)
CP_data = CP_data.ix[st:et]
CP_data = CP_data.applymap( lambda x : x * 0.01 ) # convert from percentage to decimals
CP_data = CP_data.resample('M').copy() # down sample to monthly frequency

def compute_forward_rates( df ):

    # compute the new columns for the forward rates
    df['p0'] = 0.0 # for convenience later when computing the forward rates.
    for n in xrange(1,6):
        ns = str(n)
        colname = 'SVENY0' + ns
        df['p' + ns] = - df[colname] * float(n) # this is the log-price p(n)

    for n in xrange(1,6):
        ns   = str(n)
        nm1s = str(n-1)
        df['f' + ns] = df[ 'p' + nm1s] - df[ 'p' + ns ]
        
compute_forward_rates( CP_data )
In [23]:
def excess_returns( df ):
    xrets = {}
    time_points = []
    holding_period = 365
    
    for n, (ts, series) in enumerate(df.iterrows()):
        one_yr = ts.date() + datetime.timedelta( days = holding_period )
    
        idx = df.index.searchsorted( one_yr )
        if idx == len( df ):
            break
            
        for n in xrange(2, 6):
            ns, nm1s = str(n), str(n - 1)
            
            if df.index[idx] == one_yr:
                next_yr_bond = df[ 'p' + nm1s ].iloc[ idx ]                
            else:
                next_yr_bond = df[ 'p' + nm1s ].iloc[ idx - 1 ]
            
            xrets.setdefault( n, [] ).append( next_yr_bond - series[ 'p' + ns ] - series['f1'] )
            
        time_points.append( ts.date() )
            
    return xrets, time_points

ret, xs = excess_returns( CP_data )
In [24]:
fig = plt.figure( figsize = ( 12, 20 ) )
for i in range(1, 5):
    cur_rets = np.array( ret[i+1] )
    sr = cur_rets.mean() / cur_rets.std()
    plt.subplot(410 + i)
    plt.plot( xs, ret[i+1], lw='2.0')
    plt.title( 'Excess returns from investing in the %s-yr bond. Sharpe ratio = %s' % ( i+1, sr ) )
In [25]:
def make_regression_data( df, ret, timepts ):
    d = dict( (x, df[x][timepts]) for x in ['f1', 'f2', 'f3', 'f4','f5'] )
    for x in range( 2, 6 ):
        d['rx' + str(x)] = ret[x]
    
    return pd.DataFrame(d, index=xs)
    
regression_data = make_regression_data( CP_data, ret, xs )
regression_data.head()
Out[25]:
f1 f2 f3 f4 f5 rx2 rx3 rx4 rx5
1964-01-31 0.035507 0.040546 0.041772 0.041893 0.041903 0.002228 0.003552 0.004293 0.004958
1964-02-29 0.035441 0.040580 0.041392 0.041493 0.041504 0.002871 0.004910 0.005513 0.005296
1964-03-31 0.036966 0.042144 0.042184 0.042184 0.042184 0.003533 0.005960 0.007210 0.007673
1964-04-30 0.037072 0.041786 0.042384 0.042474 0.042488 0.002843 0.004738 0.005921 0.006691
1964-05-31 0.035948 0.040780 0.041818 0.042026 0.042068 0.002384 0.003983 0.004729 0.005061
In [26]:
models = {}
for x in range( 2, 6 ):
    models[x] = ols( 'rx%s ~ f1 + f2 + f3 + f4 + f5' % x, regression_data ).fit()
models[2].summary()
Out[26]:
OLS Regression Results
Dep. Variable: rx2 R-squared: 0.212
Model: OLS Adj. R-squared: 0.203
Method: Least Squares F-statistic: 24.84
Date: Sat, 22 Oct 2016 Prob (F-statistic): 3.46e-22
Time: 12:05:32 Log-Likelihood: 1268.1
No. Observations: 468 AIC: -2524.
Df Residuals: 462 BIC: -2499.
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [95.0% Conf. Int.]
Intercept -0.0118 0.003 -4.120 0.000 -0.017 -0.006
f1 -1.2430 0.249 -4.988 0.000 -1.733 -0.753
f2 2.1554 1.716 1.256 0.210 -1.217 5.527
f3 -2.3948 5.038 -0.475 0.635 -12.295 7.506
f4 3.5693 6.234 0.573 0.567 -8.681 15.820
f5 -1.9102 2.718 -0.703 0.483 -7.252 3.431
Omnibus: 43.042 Durbin-Watson: 0.139
Prob(Omnibus): 0.000 Jarque-Bera (JB): 77.895
Skew: -0.572 Prob(JB): 1.22e-17
Kurtosis: 4.639 Cond. No. 1.15e+04

Note the large value of 643 for the condition number above. This shows that there is multicollinearity in the data.

Regression Results

Now we show the results of multiple linear regression for the excess returns on 2, 3, 4, 5 yr bonds.

In [27]:
fig = plt.figure( figsize = ( 8, 6 ) )
for x in range( 2, 6 ):
    plt.plot(range(1,6), models[x].params.values[1:], lw=2, label="%s-yr" % x )
    plt.title('Regression coefficients for excess returns.')
plt.legend()
plt.show()

Rather surprisingly, we didn't get the tent shape as described by Cochrane and Piazzesi!

Instead, we have a 'bat' shaped diagram for the regression coefficients for the five forward rates. However, this is also not unexpected. It has been chronicled in Dai and Singleton [3] that slight differences in yield curve interpolation methods will produce very different shapes for the regression coefficients. In fact, we reproduced exactly the "bat" shape as described in their work in [3].

The reason for this is due to the multicollinearity issue in the forward rates used for the multiple linear regression. See [3] and [4] for a detailed discussion.

Regression Results - Cochrane and Piazzesi

We now retry the above using CRSP data as that in Cochrane and Piazzesi.

In [28]:
CRSP_data_loc = 'CRSP_yield.csv'
CRSP_data = pd.DataFrame.from_csv(CRSP_data_loc)
CRSP_data.sort_index(inplace=True)
CRSP_data = CRSP_data.ix[st:et]
CRSP_data = CRSP_data.applymap( lambda x : x * 0.01 ) # convert from percentage to decimals
CRSP_data.head()
Out[28]:
SVENY01 SVENY02 SVENY03 SVENY04 SVENY05
Date
1964-01-31 0.03789 0.03893 0.03939 0.03943 0.03984
1964-02-28 0.03947 0.03981 0.04025 0.04010 0.04023
1964-03-31 0.03910 0.04076 0.04125 0.04145 0.04073
1964-04-30 0.03783 0.03978 0.04047 0.04082 0.04077
1964-05-28 0.03832 0.03918 0.03964 0.04021 0.03979
In [29]:
compute_forward_rates( CRSP_data )
ret, xs = excess_returns( CRSP_data )

fig = plt.figure( figsize = ( 12, 20 ) )
for i in range(1, 5):
    cur_rets = np.array( ret[i+1] )
    sr = cur_rets.mean() / cur_rets.std()    
    plt.subplot(410 + i)
    plt.plot( xs, ret[i+1], lw='2.0')
    plt.title( 'Excess returns from investing in the %s-yr bond. Sharpe ratio = %s' % ( i+1, sr ) )
In [30]:
regression_data = make_regression_data( CRSP_data, ret, xs )
display( HTML( regression_data.head().to_html() ) ) 

models = {}
for x in range( 2, 6 ):
    models[x] = ols( 'rx%s ~ f1 + f2 + f3 + f4 + f5' % x, regression_data ).fit()
models[2].summary()
f1 f2 f3 f4 f5 rx2 rx3 rx4 rx5
1964-01-31 0.03789 0.03997 0.04031 0.03955 0.04148 -0.00080 0.00172 0.00163 0.00235
1964-02-28 0.03947 0.04015 0.04113 0.03965 0.04075 -0.00189 0.00060 0.00102 0.00160
1964-03-31 0.03910 0.04242 0.04223 0.04205 0.03785 0.00038 0.00397 0.00679 0.00447
1964-04-30 0.03783 0.04173 0.04185 0.04187 0.04057 0.00033 0.00306 0.00473 0.00590
1964-05-28 0.03832 0.04004 0.04056 0.04192 0.03811 -0.00093 0.00010 0.00135 -0.00013
Out[30]:
OLS Regression Results
Dep. Variable: rx2 R-squared: 0.339
Model: OLS Adj. R-squared: 0.332
Method: Least Squares F-statistic: 47.34
Date: Sat, 22 Oct 2016 Prob (F-statistic): 1.67e-39
Time: 12:05:55 Log-Likelihood: 1294.2
No. Observations: 468 AIC: -2576.
Df Residuals: 462 BIC: -2552.
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [95.0% Conf. Int.]
Intercept -0.0159 0.003 -6.013 0.000 -0.021 -0.011
f1 -0.9807 0.121 -8.083 0.000 -1.219 -0.742
f2 0.5493 0.252 2.180 0.030 0.054 1.044
f3 1.2118 0.209 5.795 0.000 0.801 1.623
f4 0.3024 0.154 1.963 0.050 -0.000 0.605
f5 -0.8640 0.130 -6.642 0.000 -1.120 -0.608
Omnibus: 19.533 Durbin-Watson: 0.342
Prob(Omnibus): 0.000 Jarque-Bera (JB): 27.025
Skew: -0.358 Prob(JB): 1.35e-06
Kurtosis: 3.935 Cond. No. 414.
In [31]:
fig = plt.figure( figsize = ( 8, 6 ) )
for x in range( 2, 6 ):
    plt.plot(range(1,6), models[x].params.values[1:], lw=2, label="%s-yr" % x )
    plt.title('Regression coefficients for excess returns.')
plt.legend()
plt.show()

So indeed we have reproduced the 'tent' shape as in Cochrane and Piazzesi. This shows how sensitive the regression betas are due to the multicolinearity problem and the issue of the smoothing of the yield curve. In Cochrane and Paizzesi the yield curve data is unsmoothed, while the quandl SVENY data set uses the Nelson-Sigel smoothed yield curve construction.

References

[1] Campbell J, Shiller, R, (1991), Yield Spreads and Interest Rate Movements: A Bird’s Eye View, Review of Economic Studies, 58, 495-514.

[2] Cochrane J H, Piazzesi M, (2005), Bond Risk Premia, American Economic Review, Vol 95, no 1, 138-160 — see also https://www.aeaweb.org/aer/data/mar05_app_cochrane.pdf (last accessed 25 November, 2014) for technical details, and the regression coefficients.

[3] Dai, Q, Singhleton K J, YAng W, (2004), Predictability of Bond Risk Premia and Affine Term Structure Models, working paper, Stern School of Business.

[4] Rebonato Riccardo, Return-predicting factors for us treasuries : on the similarity of “tents” and “bats”, International Journal of Theoretical and Applied Finance.- River Edge, NJ [u.a.] : World Scientific, ISSN 0219-0249, ZDB-ID 14289829. - Vol. 18.2015, 4, p. 1-14