Convexity for historical bond data

Question

I'm trying to write a program to calculate the convexity of a bond. The bigger idea is, that if I have access to the actual price for each point in time, I should be able to calculate various features of the bond given that date. The convexity program is written as follows:

""" Calculate convexity of a bond """
    def bond_convexity(price, par, T, coup, freq, dy=0.0001):
        ytm = bond_ytm(price, par, T, coup, freq)

        ytm_minus = ytm - dy
        price_minus = bond_price(par, T, ytm_minus, coup, freq)

        ytm_plus = ytm + dy
        price_plus = bond_price(par, T, ytm_plus, coup, freq)

        convexity = (price_minus+price_plus-2*price)/(price*dy**2)
        return convexity

Where the yield to maturity and bond price functions are as follows:

"""Interpolate yield to maturity"""
def bond_ytm(price, par, coupdates, coup, freq=2, guess=0.05):
     freq = float(freq)
     MAT = coupdates[len(coupdates) - 1] #maturity date
     timedif = list() #init
     periods = list() #init
     for i in range(0,len(coupdates)-1):
         timedif.append( round((MAT - coupdates[i]).days/365,1) )     #collection of time to maturity - coupon date
         periods.append(timedif[i]*freq)
     coupon = coup/100.*par/freq
     dt = [(x)/freq for x in periods]
     ytm_func = ytm_func = lambda y : \
         sum([coupon/(1+y/freq)**(freq*t) for t in dt]) + par/(1+y/freq)**\
(freq*max(dt)) - price

 return optimize.newton(ytm_func, guess)

 """Price bond"""
 def bond_price(par, coupdates, ytm, coup, freq=2):
      freq = float(freq)
      MAT = coupdates[len(coupdates) - 1] #maturity date
      timedif = list() #init
      periods = list() #init
      for i in range(0,len(coupdates)-1):
          timedif.append( round((MAT - coupdates[i]).days/365,1) ) #collection of time to maturity - coupon date
          periods.append(timedif[i]*freq)
      coupon = coup/100.*par/freq
      dt = [(x)/freq for x in periods]
      price = sum([coupon/(1+ytm/freq)**(freq*t) for t in dt]) + par/(1+ytm/freq)**(freq*max(dt))
 return price

Now, the following data shows the upcoming coupon dates for the bond and today's date at the top. I've saved this variable as cdd:

cdd
Out[28]: 
[datetime.datetime(2017, 10, 31, 15, 20, 0, 480212),
 datetime.datetime(2018, 2, 1, 0, 0),
 datetime.datetime(2018, 8, 1, 0, 0),
 datetime.datetime(2019, 2, 1, 0, 0),
 datetime.datetime(2019, 8, 1, 0, 0),
 datetime.datetime(2020, 2, 1, 0, 0),
 datetime.datetime(2020, 8, 1, 0, 0),
 datetime.datetime(2021, 2, 1, 0, 0),
 datetime.datetime(2021, 8, 1, 0, 0),
 datetime.datetime(2022, 2, 1, 0, 0),
 datetime.datetime(2022, 8, 1, 0, 0),
 datetime.datetime(2023, 2, 1, 0, 0),
 datetime.datetime(2023, 8, 1, 0, 0),
 datetime.datetime(2024, 2, 1, 0, 0),
 datetime.datetime(2024, 8, 1, 0, 0),
 datetime.datetime(2025, 2, 1, 0, 0),
 datetime.datetime(2025, 8, 1, 0, 0),
 datetime.datetime(2026, 2, 1, 0, 0),
 datetime.datetime(2026, 8, 1, 0, 0),
 datetime.datetime(2027, 2, 1, 0, 0),
 datetime.datetime(2027, 8, 1, 0, 0),
 datetime.datetime(2028, 2, 1, 0, 0),
 datetime.datetime(2028, 8, 1, 0, 0),
 datetime.datetime(2029, 2, 1, 0, 0),
 datetime.datetime(2029, 8, 1, 0, 0),
 datetime.datetime(2030, 2, 1, 0, 0),
 datetime.datetime(2030, 8, 1, 0, 0),
 datetime.datetime(2031, 2, 1, 0, 0),
 datetime.datetime(2031, 8, 1, 0, 0),
 datetime.datetime(2032, 2, 1, 0, 0),
 datetime.datetime(2032, 8, 1, 0, 0),
 datetime.datetime(2033, 2, 1, 0, 0),
 datetime.datetime(2033, 8, 1, 0, 0),
 datetime.datetime(2034, 2, 1, 0, 0),
 datetime.datetime(2034, 8, 1, 0, 0),
 datetime.datetime(2035, 2, 1, 0, 0),
 datetime.datetime(2035, 8, 1, 0, 0),
 datetime.datetime(2036, 2, 1, 0, 0),
 datetime.datetime(2036, 8, 1, 0, 0),
 datetime.datetime(2037, 2, 1, 0, 0),
 datetime.datetime(2037, 8, 1, 0, 0),
 datetime.datetime(2038, 2, 1, 0, 0),
 datetime.datetime(2038, 8, 1, 0, 0),
 datetime.datetime(2039, 2, 1, 0, 0),
 datetime.datetime(2039, 8, 1, 0, 0),
 datetime.datetime(2040, 2, 1, 0, 0),
 datetime.datetime(2040, 8, 1, 0, 0),
 datetime.datetime(2041, 2, 1, 0, 0),
 datetime.datetime(2041, 8, 1, 0, 0),
 datetime.datetime(2042, 2, 1, 0, 0),
 datetime.datetime(2042, 8, 1, 0, 0),
 datetime.datetime(2043, 2, 1, 0, 0),
 datetime.datetime(2043, 8, 1, 0, 0),
 datetime.datetime(2044, 2, 1, 0, 0),
 datetime.datetime(2044, 8, 1, 0, 0),
 datetime.datetime(2045, 2, 1, 0, 0),
 datetime.datetime(2045, 8, 1, 0, 0)]

However, if I run the program to calculate convexity for a bond priced at 112.057, par = 100, dates as cdd, coupon of 4.9, and semi annual frequency, I get the following solution:

bond_convexity(112.057, 100, cdd, 4.9, 2)
Out[32]: 351.98162487756656

Bloomberg is showing that the convexity should be approximately 3.467. What could I possibly be doing wrong? Any help is appreciated.

Answer 1

A couple of thoughts:

1. There are different ways to express convexity, depending how you express yields. An example might help: If yields are expressed in percent (e.g., 5 for 5%), then the convexity of a zero coupon bond with a duration of 15 is roughly $15\times15/100 = 2.25$. Alternatively, if you expressed yields in decimals (0.05 for 5%), then the corresponding convexity is expressed as $15\times 15 = 225$. The first approach is more common and is used by BBG. All you have to do is to divide your number by 100 to be consistent.

2. You are using an approximation formula for convexity, but yield-based convexity for vanilla bonds has an exact closed-form solution (simply take the second derivative of the price/yield formula relative to yield).

3. Your implementation does not conform to market conventions (stuff like Actual/Actual day count). There are a lot of subtleties involved in bond pricing and I recommend using a professional implemented solution if precision is important. QuantLib, for example, has a very easy to use Python binding that you can use.