Empirical duration and convexity for bonds using linear regression

Question

I have a given time series of bond yields from Quandl. From the time series, I have taken a sample to simulate a path of bond yields by Monte Carlo in Python.

I have to do the following task:

"Calculate the average empirical duration and convexity of any one of the bond indices for the sample period.

The following approximate relationship exists between yields changes and total returns (log changes in total return index):

$$\ln \left( \frac{TRI_{t}}{TRI_{t-1}} \right)≈Q -D(y_{t}-y_{t-1} )+\frac{C}{2} (y_{t}-y_{t-1} )^2$$

Where $TRI$ is the total return index, $y$ is the yield, $Q$ is the empirical carry, $D$ is the empirical duration, and $C$ is the empirical convexity. This can be expressed as a simple linear regression of the form $y = a + bx + cz$. Estimating the parameters will help you determine the required quantities. Note: yields are in percentage terms so converting to decimal, or alternatively scaling the left-hand side by 100, will give sensible duration numbers."

I don't know anything about finance, or regression. A simple google search did not show the meaning of "empirical carry", so I don't really understand how to approach the question. Since the given relationship above has a quadratic term (in bond yields $y$), how can it fit the linear form $y = a +bx + cz$?

A general explanation of how I should approach this would be appreciated. I need to do this in Python (using pandas, numpy, all the standard packages). Thanks!

Answer 1

It really doesn't matter if the term $(y_t - y_{t-1})^2$ is present, as long as there are no terms in the model that involve non-linear operations of the weights, you are good to go. As an example if the model is something like $u_t = a + b e^a x_t$ then you're in trouble, but fortunately this is not your case.

Now, to take care of the non-linearity in the features, just define a new set of variable, e.g.

\begin{eqnarray} u_t &=& \ln \frac{{\rm TRI}_t}{{\rm TRI}_{t-1}} \\ x_t^{(1)} &=& y_t - y_{t-1} \\ x_t^{(2)} &=& (y_t - y_{t-1})^2 \tag{1} \end{eqnarray}

So that the model becomes

$$ u_t = \beta_0 + \beta_{1} x_t^{(1)} + \beta_2 x_t^{(2)} \tag{2} $$

with

$$ \beta_0 = Q, ~~~~~ \beta_1 = -D ~~~\mbox{and}~~~ \beta_2 = C/2 \tag{3} $$

The problem now is how to estimate the $\beta$s. I assume you have a pandas dataframe df with columns called $y$ and ${\rm TRI}$. This is the what you should do

import pandas
import numpy as np
import statsmodels.formula.api as sm

# (read data)

# prepare features
u = np.log(df['TRI']).diff(1)[1 : ]
x1 = df['y'].diff(1)[1 : ]
x2 = (df['y'].diff(1)[1 : ]) ** 2
dfnew = pandas.DataFrame({'u' : u, 'x1' : x1, 'x2' : x2})

# fit
result = sm.ols(formula="u ~ x1 + x2", data = dfnew).fit()

# weights
Q = result.params['Intercept']
D = -result.params['x1']
C = 2 * result.params['x2']
print 'Q = {}, D = {}, C = {}'.format(Q, D, C)