# Applying Differential Evolution to the Nelson Siegel Model in Python

I am trying to create a zero curve from a series of government bonds by minimizing the differences in the dirty prices.

The problem is that there is something wrong with my differential evolution algorithm.

My process is as follows:

1. I have utilized quantlib to generate the cash flow schedules and year fractions
2. I have then moved on to a function that attempts to discount the cash flows from from the zero rate approximated by the nelson siegel function.
3. Then use scipy's differential evolution algorithm to so solve for the 4 parameters

The outputs that I get is far from the ones I would expect. Given, the below data : I'd expect:

• b1 to be closer to the long term rate ~= (0.0267)
• B2 to be closer to the negative of the slope ~= (-0.0255)

Would anyone know where I am going wrong?

import QuantLib as ql
import pandas as pd
import numpy as np
import datetime
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution
import functools

def text_to_pd(text_input,isfut=0):
parsed_data = [tmp_list.split('\t') for tmp_list in text_input.splitlines() if len(tmp_list)>0]
columns = parsed_data.pop(0)
ds = pd.DataFrame(parsed_data,columns=columns)
return ds

def gen_cf_schedule(bond,name=None,effective_date=pd.to_datetime("today")):

if name == None: name = 'tmp_bond'

cashflows = pd.DataFrame({
'date': pd.to_datetime(cf.date().to_date()),
name: cf.amount()
} for cf in bond.cashflows()
)
cashflows =cashflows[cashflows.date >= effective_date]
cashflows = cashflows.set_index('date').groupby('date').sum()
return cashflows

def gen_cf_schedule_many_bonds(list_bonds,list_bond_name):

if len(list_bonds) != len(list_bond_name):
print('Error : Bond list and name list are not of same size')

cf_list = [
gen_cf_schedule(tmp_bond,tmp_name)
for tmp_bond,tmp_name in zip(list_bonds,list_bond_name)
]

# below starts from first 2 pairs in list and does an outter join of call cash flows
df_final = functools.reduce(
lambda left,right: pd.merge(left,right,how = 'outer',left_index=True, right_index=True), cf_list)

return df_final

data = """
MATURITY    CPN CPN_FREQ    DIRTY_PRICE YLD
15/07/2022  5.75    2   102.8251133 0.114989457
21/11/2022  2.25    2   101.9755    0.350344935
21/04/2023  5.5 2   107.426 0.664767214
21/04/2024  2.75    2   104.4845    1.06749642
21/11/2024  0.25    2   96.9755 1.386223437
21/04/2025  3.25    2   106.31  1.539978867
21/11/2025  0.25    2   94.8665 1.675063719
21/04/2026  4.25    2   111.523 1.715012175
21/09/2026  0.5 2   94.583  1.775096919
21/04/2027  4.75    2   115.9925    1.806204173
21/11/2027  2.75    2   105.409 1.871278842
21/05/2028  2.25    2   102.4795    1.915042694
21/11/2028  2.75    2   105.707 1.94498807
21/04/2029  3.25    2   109.5315    1.972445794
21/11/2029  2.75    2   105.9405    2.007513096
21/05/2030  2.5 2   104.108 2.032546486
21/12/2030  1   2   91.543  2.06749518
21/06/2031  1.5 2   95.3565 2.07373467
21/11/2031  1   2   90.6355 2.089994302
21/05/2032  1.25    2   92.42   2.105042334
21/11/2032  1.75    2   96.825  2.122473219
21/04/2033  4.5 2   125.034 2.116228619
21/06/2035  2.75    2   106.0135    2.2600066
21/04/2037  3.75    2   118.854 2.357477112
21/06/2039  3.25    2   111.195 2.484967213
21/05/2041  2.75    2   103.3825    2.567442991
21/03/2047  3   2   107.288 2.664936079
21/06/2051  1.75    2   81.654  2.667453281

"""

bonds_data = text_to_pd(data)
#Convert string to quantlib date
date_format = '%d/%m/%Y'
lambda_str_to_qlDate = lambda x : ql.DateParser.parseFormatted(x, date_format)
bonds_data = bonds_data.assign(MATURITY = lambda df_: df_['MATURITY'].apply(lambda_str_to_qlDate))

#parse to numeric
cols_to_numeric = [i for i in bonds_data.columns  if i!='MATURITY'  ]
bonds_data[cols_to_numeric] = bonds_data[cols_to_numeric].apply(pd.to_numeric)
coupon_frequency

# bond conventions
# Conventions used for Aus sov bonds
calendar = ql.Australia()
day_count =  ql.ActualActual(ql.ActualActual.ISMA)
settlement_days = 2
coupon_frequency   = ql.Period(ql.Semiannual)
​effective_date = ql.Date(9,2,2022)
ql.Settings.instance().evaluationDate=settlement_date

# to generate bonds
bonds = [    # creating the bond
ql.FixedRateBond(
settlement_days,
calendar,
100.0,
calendar.advance(effective_date,ql.Period(-6,ql.Months)), # ensure a previous coupon to calc accrued amount
maturity,
coupon_frequency,
[coupon/100],
day_count,
100,
ql.Date(), ql.Date(), # default dates,
ql.DateGeneration.Backward,
False, # default, end of month
calendar,
ql.Period(7,ql.Days),
calendar,
)

for maturity,coupon in zip(bonds_data.MATURITY,bonds_data.CPN)
]

#cashflows:
pd_cf=gen_cf_schedule_many_bonds(bonds,bonds_data.index)

#Now to simply generate t (yearfractions)
tmp_format = '%d-%m-%Y'
list_cf_dates = [
ql.DateParser.parseFormatted(x.strftime('%d-%m-%Y'), '%d-%m-%Y')
for x in pd_cf.index
]

year_frac = [day_count.yearFraction(settlement_date,cf_date) for cf_date in list_cf_dates]

dirty_prices = bonds_data.DIRTY_PRICE.to_numpy()



The task next is to create the objective function, which is defined below:

def NS(x):
x1,x2,x3,lam1 = x

tmp_ns = lambda x1,x2,x3,lam1,t : x1 + x2*(lam1/t)*((1-np.exp(-t/lam1))) + x3*((lam1/t)*(1-np.exp(-t/lam1)) -np.exp(-t/lam1))

#exp(-t*r)
neg_t_times_r  = [-t*(tmp_ns(x1,x2,x3,lam1,t)) for t in year_frac ]

# CF.exp(-t*r)
cf_discounted = pd_cf.apply(lambda x : x*np.exp(neg_t_times_r),axis=0)
test_prices = cf_discounted.sum(axis=0)
test_prices = test_prices.to_numpy()

squared_dif = [(a-b)**2 for a,b in zip(dirty_prices,test_prices)]
sum_squared_diff = sum(squared_dif)

return sum_squared_diff



At this stage, I have set the bonds and attempted to run the code as per below, but the fit of the curve is particularly bad and the parameters are far from what I expect them to be.

Note that I have made the constraints exceptionally wide, but even when I narrow them, the results seem to be off.

bounds_ns=[(0,0.5),(-0.5,0.5),(-1,1),(0,100)]
result = differential_evolution(NS_Proper,bounds=bounds_ns,maxiter=1000_000).x
print('b_1: {} '.format(result[0]))
print('b_2: {}'.format(result[1]))
print('b_3: {}'.format(result[2]))
print('b_4: {}'.format(result[3]))


The result is as below :

b_1: 0.0

b_2: 0.01145940645447193

b_3: 0.07370270399890494

b_4: 22.21196983881775

Even when plotted it is obviously wrong:

year_frac_bond_maturity = [day_count.yearFraction(settlement_date,tmp_date) for tmp_date in bonds_data.MATURITY]
tmp_ns = lambda x1,x2,x3,lam1,t : x1 + x2*(lam1/t)*((1-np.exp(-t/lam1))) + x3*((lam1/t)*(1-np.exp(-t/lam1)) -np.exp(-t/lam1))
zero_ns = [tmp_ns(*[2.53546089e-04, 1.13003488e-02, 7.32022467e-02, 2.23133786e+01],t) for t in year_frac_bond_maturity]
Fig,ax = plt.subplots()
Fig.set_size_inches((15,5))
plt.plot(year_frac_bond_maturity,zero_ns,label='zero_ns')
plt.plot(year_frac_bond_maturity,bonds_data.YLD/100,label='yield',marker='o',linestyle="None")
plt.legend()
plt.show()


• Instead of fitting to price errors, have you tried fitting to yield errors (or price errors divided by duration^2 which might be faster)? Feb 10, 2022 at 6:14
• Hi Helin. Yes, I had extended it to convert the prices into ytm and then taken the yield errors as suggested. I ran into a similar issue. The reason I put price errors here was because I wanted to keep the code a bit short. Let me try with the price errors divided by duration^2 and I'll revert. Thanks Feb 10, 2022 at 6:21
• Helin, I rebuilt the function that was minimizing the ytm. It finally worked. I have a suspicion that the prior error was coming from not converting to the clean price before calculating the ytm. Thanks for this. Feb 12, 2022 at 11:02
• The reason for the poor fit in price space is that shorter term bonds have low duration, so even a large yield error might not translate into a sufficiently large price error for the optimizer to correct for. Minimizing yield error is more direct. Minimizing price / duration^2 is roughly the same as minimizing yield error (first order), but it's faster because you don't need to calculate yields at each optimization step, which can be slow. Of course, if you don't care about speed and don't have a ton of bonds, minimizing yield error probably produces better results =) Feb 14, 2022 at 6:38

If we want to create a zero-coupon curve we can take the market yield curve and bootstrap zero-coupon curve. There are a few complexities which make this difficult like non-annual coupons and tenors not falling on exact years but I think I figured a way around it.

I had great difficulty finding the answer to this online after a great deal of looking so I hope this will be useful to others.

Note the difficulty here is dealing with non-annual coupon payments and real-world dates that do not match round years.

What we are essentially doing when creating a zero-coupon curve is creating a curve that gives us a discounting rate for various tenors (maturities) that is not 'dirtied' by the differing coupons of each of the bonds that have been issued. Given that bonds with a higher coupon will have a lower duration (you get your money back quicker so there is less risk) we need to make an adjustment. Zero-coupon curves allow us to assess relative value without the distortion of the coupon's impact on duration.

So in simple terms the steps to take are:

1. Get the yield to maturity and tenor (in years) for each bond for the issuer. Interpolate to fit a curve to the points (e.g. Nelson Siegel OLS regression) which will give you the parameters to get the corresponding yield for any/all maturity/tenors.
2. Build a full curve of say 30 years at semi-annual increments using the Nelson Siegel formula and the parameters that were optimized above by inputting each time value (0.5, 1, 1.5, ..., 29.5, 30) into the NS formula with the parameters. You now have the 'par curve' for semi-annual maturities going out 30 years.
3. The yields at a tenor of 0.5 years calculated above is a zero-coupon rate and your starting point for bootstrapping the zero-coupon curve.
4. We then use bootstrapping to construct the zero/spot curve. We use the interpolated yield for each tenor as the ANNUAL COUPON which defines the cash flows before maturity. The price to which we discount each set of cashflows is $100 (or par value). There are plenty of explanations of bootstrapping online but can post here if necessary. Why is the coupon and price determined like this you ask? because we can say that if the issuer was to issue a new bond today for that tenor they would be issuing at the interpolated yield to maturity and an issue price of$100 (par) (therefore the coupon is equal to the yield).

I hope this makes sense - please ask any questions or let me know where I have gone wrong.

• I'm sorry but this is not how you use the NS model: 1) there are no issues with different bonds' cashflows landing on the same dates; 2) it is completely plausible for two bonds with different coupons to both trade at fair value, in spite of different durations, so the adjustment is not needed (unless there are legitimate market segmentation/liquidity concerns); 3) there really should only be one step – strip out the cash flows of every bond, discount them using the NS discount curve, optimize the parameters so that all the bonds are reasonably priced by the single discount curve. Feb 14, 2022 at 6:32
• @powpow's code is actually the correct way; just needed a small change in the optimization target. Feb 14, 2022 at 6:48
• I think you might have misunderstood what was being said here @helin. The above answer is for those looking to not use the quantlib library.
– JPI
Feb 15, 2022 at 7:22
• Hi @JPI I didn’t. You don’t need Quantlib for curve fitting. Also powpow didn’t use quantlib for the NS part; only for cash flow generation. Feb 15, 2022 at 11:39
• Here's a good paper from the Fed for reference and illustrate how the model is used in real life federalreserve.gov/pubs/feds/2006/200628/200628pap.pdf Feb 16, 2022 at 5:49