QuantLib 金融计算——收益率曲线之构建曲线(3)

载入 QuantLib 和其他包:

import QuantLib as ql
import seaborn as sb
import numpy as np
import pandas as pd print(ql.__version__)
1.15

概述

本文展示利用 quantlib-python 根据样本券的交易数据估算出即期利率的期限结构的完整流程,并指出当前实现所存在的问题。

示例所用的样本券交易数据来自专门进行期限结构分析的 R 包——termstrc。具体来说是数据集 govbonds 中的 GERMANY 部分,包含 2008-01-30 这一天德国市场上 52 只固息债的成交数据。

估算期限结构的步骤

QuantLib 中估算期限结构的核心流程有两步:

  1. 配置 *Helper 对象,描述样本券信息,包括付息时间表(schedule)、价格(默认用净价)、票息等;
  2. 配置期限结构模型,可以额外提供样本券权重、优化方法、参数正则化条件等参数辅助计算。

读取样本券数据

govBond = pd.read_csv(
'GERMANY_INFO.csv',
parse_dates=['MATURITYDATE', 'ISSUEDATE']) numberOfBonds = govBond.shape[0] PRICE = [
ql.QuoteHandle(ql.SimpleQuote(p)) for p in govBond['PRICE']]
MATURITYDATE = [
ql.Date(m.day, m.month, m.year) for m in govBond['MATURITYDATE']]
ISSUEDATE = [
ql.Date(i.day, i.month, i.year) for i in govBond['ISSUEDATE']]
COUPONRATE = govBond['COUPONRATE'].values

一些基本配置

# 查看 govbonds 数据集可知样本券均为每年付息一次
frequency = ql.Annual
# termstrc 的日期计算并不如 QuantLib 精细,
# 为了和 termstrc 的算法保持一致,示例使用如下天数计算规则
dc = ql.Actual365Fixed(ql.Actual365Fixed.Standard)
paymentConv = ql.Unadjusted
terminationDateConvention = ql.Unadjusted
convention = ql.Unadjusted
redemption = 100.0
faceAmount = 100.0
# 其实我不知道样本券所在的交易所,
# 所以不确定是不是该用这个日历 :)
calendar = ql.Germany(ql.Germany.Eurex) # 估值日期 2008-01-30
today = calendar.adjust(ql.Date(30, 1, 2008))
ql.Settings.instance().evaluationDate = today # 为了和 termstrc 的算法保持一致,示例采用 T+0 结算
bondSettlementDays = 0
bondSettlementDate = calendar.advance(
today,
ql.Period(bondSettlementDays, ql.Days))

配置 *Helper 对象

instruments = []

for j in range(numberOfBonds):
# 配置付息时间表
schedule = ql.Schedule(
ISSUEDATE[j],
MATURITYDATE[j],
ql.Period(frequency),
calendar,
convention,
terminationDateConvention,
ql.DateGeneration.Backward,
False) # 配置 Helper 对象
# 因为样本券均为固息债,所以采用 FixedRateBondHelper 类
# 对于其他金融工具,需要使用对应的 Helper 类
helper = ql.FixedRateBondHelper(
PRICE[j],
bondSettlementDays,
faceAmount,
schedule,
[COUPONRATE[j]],
dc,
paymentConv,
redemption) instruments.append(helper)

配置期限结构

tolerance = 1.0e-6
max = 5000 # 即期利率的 Svensson 模型
sf = ql.SvenssonFitting()
# 即期利率的 Nelson Siegel 模型
nsf = ql.NelsonSiegelFitting()
# 用指数样条函数拟合贴现因子
esf = ql.ExponentialSplinesFitting()
# 用简单多项式函数拟合贴现因子
spf = ql.SimplePolynomialFitting(8)
# 用三次 B-样条函数拟合贴现因子
knots = [-20.0, -10.0, 0.0, 0.25, 0.5, 1, 3, 5, 10, 20, 30, 40, 50]
cbsf = ql.CubicBSplinesFitting(knots) tsSvensson = ql.FittedBondDiscountCurve(
bondSettlementDate, instruments, dc,
sf,
tolerance, max) tsNelsonSiegel = ql.FittedBondDiscountCurve(
bondSettlementDate, instruments, dc,
nsf,
tolerance, max) tsExponentialSplines = ql.FittedBondDiscountCurve(
bondSettlementDate, instruments, dc,
esf,
tolerance, max) tsSimplePolynomial = ql.FittedBondDiscountCurve(
bondSettlementDate, instruments, dc,
spf,
tolerance, max) tsCubicBSplines = ql.FittedBondDiscountCurve(
bondSettlementDate, instruments, dc,
cbsf,
tolerance, max)

估算期限结构

sv = []
ns = []
es = []
sp = []
cbs = []
matList = [] matDate = bondSettlementDate while matDate <= bondSettlementDate + ql.Period(31, ql.Years):
matDate = matDate + ql.Period(1, ql.Days)
matList.append(
dc.yearFraction(bondSettlementDate, matDate)) sv.append(
tsSvensson.zeroRate(matDate, dc, ql.Continuous, frequency).rate() * 100)
ns.append(
tsNelsonSiegel.zeroRate(matDate, dc, ql.Continuous, frequency).rate() * 100)
es.append(
tsExponentialSplines.zeroRate(matDate, dc, ql.Continuous, frequency).rate() * 100)
sp.append(
tsSimplePolynomial.zeroRate(matDate, dc, ql.Continuous, frequency).rate() * 100)
cbs.append(
tsCubicBSplines.zeroRate(matDate, dc, ql.Continuous, frequency).rate() * 100)

汇总结果

# 以 termstrc 的估算结果作为比较基准
beta0 = 5.017052
beta1 = -1.117214
beta2 = -3.173622
tau = 2.443936
termstrc = [
beta0 + \
beta1 * (1 - np.exp(-m / tau)) / (m / tau) + \
beta2 * ((1 - np.exp(-m / tau)) / (m / tau) - np.exp(-m / tau)) for m in matList] df = pd.DataFrame(
dict(
maturity=matList * 6,
rate=sv + ns + es + sp + cbs + termstrc,
type=np.repeat(
['Svensson', 'NelsonSiegel', 'ExponentialSplines',
'SimplePolynomial', 'CubicBSplines', 'termstrc'], len(matList)))) print(tsSvensson.fitResults().solution())
print(tsNelsonSiegel.fitResults().solution())
print(tsExponentialSplines.fitResults().solution())
print(tsSimplePolynomial.fitResults().solution())
print(tsCubicBSplines.fitResults().solution()) print(tsSvensson.fitResults().minimumCostValue())
print(tsNelsonSiegel.fitResults().minimumCostValue())
print(tsExponentialSplines.fitResults().minimumCostValue())
print(tsSimplePolynomial.fitResults().minimumCostValue())
print(tsCubicBSplines.fitResults().minimumCostValue()) sb.relplot(
x='maturity', y='rate', kind='line', hue='type',
size='type', sizes=[2, 2, 2, 2, 2, 4],
data=df, height=5, aspect=1.6)
[ -134.509; 134.548; 136.375; -0.0234826; 0.000905134; 0.422332 ]
[ -9.7222; 9.75752; 38.0528; 3.14624e-05 ]
[ -31238.6; -95904.4; 46747.5; 53652; 21350.9; 2438.62; -77566.2; 36236.2; 0.000302241 ]
[ -0.0373576; 0.00165862; -0.000203106; 9.83472e-06; -1.47693e-07 ]
[ 1.13572; 0.977843; 0.948002; 0.90216; 0.798218; 0.616569; 0.355499; 0.293158 ] 0.00032255014803192505
0.0020683101877706977
0.0019235041138322383
0.0008108667989534623
0.00037198366450668654

QuantLib 金融计算——收益率曲线之构建曲线(3)-LMLPHP

图 1:QuantLib 的结果

QuantLib 金融计算——收益率曲线之构建曲线(3)-LMLPHP

图 2:termstrc 的结果

NelsonSiegelSimplePolynomialExponentialSplines 的结果与基准相去甚远。SvenssonCubicBSplines 的结果在短端与基准非常接近,但在长端依然有明显差距,SvenssonCubicBSplines 的结果要略低于基准。

考虑到基准似乎在长端高估了真实利率水平,SvenssonCubicBSplines 的结果可能要好于基准。另外,CubicBSplines 甚至顾及到了 0 附近的两个“异常值”。

当前实现存在的问题与对策

粗看结果似乎还可以接受,但实际上经不起推敲。

  • 首先,根据 Nelson Siegel 模型和 Svensson 模型的经济意义(参考文献 1),估计结果绝对值的数量级应该和利率处于同一水平,通常是 10 以内的某个数。尤其是模型中的 \(\beta_0\),大致应该等于利率的平均值。
  • 其次,三次 B-样条的结果在远端出现了 S 形的扭曲,可能是 knot 选择不当的结果,最终结果对 knot 的选择其实非常敏感(参考文献 3)。另外,QuantLib 采用了 \(d(t) = \sum_{i=0}^{n} c_i \times N_{i,3}(t)\) 的格式(\(N_{i,3}(t)\) 是基本样条函数,\(d(t)\) 是贴现因子),因为 \(N_{i,3}(t)\) 在最外测的两个 knot 上取值非常小(MATLAB 的演示),这使得使用者必须提供而外的 knot 完全覆盖当前的期限范围才能有合理的估计,相当反人类的设计。
  • 第三,指数样条等方法的参数过于极端。

所有问题的根源是相同的,因为估算期限结构本质上是一个优化问题。以 Nelson Siegel 模型和 Svensson 模型为例,参数估计本身是一个相当有挑战性的非凸优化问题(参考文献 2),可能需要借助一些特殊的技术手段(参考文献 2),而不是依赖于某个优化算法。但是 quantlib-python 在封装期限结构接口的时候只保留了样本券权重一个自由度,优化算法、正则化条件等选项均被忽略,特别是优化算法,统一使用比较原始的单纯形算法

若要改良当前的结果,一种方法是编写 C++ 程序使用其他优化算法,并配置正则化条件;另一种方法是自定义 swig 的接口文件,修改 quantlib-python 期限结构类的接口,使其能使用其他优化算法,并接受正则化条件。

如果无法实现以上两种方法,在当前有限的条件下推荐使用三次 B-样条函数估算期限结构,但要注意 knot 的选择。

参考文献

  1. 《收益率曲线的建模和预测——基于 DNS 方法创新》,东北财经大学出版社
  2. R. Ferstl and J. Hayden, "Zero-Coupon Yield Curve Estimation with the Package termstrc" Journal of Statistical Software, August 2010, Volume 36, Issue 1.
  3. James, J. and N. Webber, "Interest Rate Modelling" John Wiley, 2000.

附录

样本券数据。

DE00011414142008-02-152002-08-140.0425100.0024.087
DE00011371312008-03-142006-03-080.0399.922.6557
DE00011414222008-04-112003-04-110.0399.8052.4262
DE00011371492008-06-132006-05-300.032599.752.069
DE00011350772008-07-041998-07-040.0475100.3052.7514
DE00011371562008-09-122006-08-300.03599.761.3579
DE00011414302008-10-102003-09-250.03599.751.0902
DE00011371642008-12-122006-11-300.037599.9750.5225
DE00011351012009-01-041999-01-040.0375100.04160.2869
DE00011371722009-03-132007-02-280.0375100.05743.3299
DE00011414482009-04-172004-02-020.032599.50492.5751
DE00011371802009-06-122007-05-300.045101.09712.877
DE00011351272009-07-041999-07-040.045101.1372.6066
DE00011371982009-09-112007-08-240.04100.71991.5628
DE00011414552009-10-092004-08-250.03599.88831.0997
DE00011372062009-12-112007-09-210.04100.9080.5683
DE00011351352010-01-041999-10-220.05375103.35530.4112
DE00011414632010-04-092005-02-240.032599.50342.6462
DE00011351502010-07-042000-05-050.0525103.9133.041
DE00011414712010-10-082005-08-260.02597.42290.7923
DE00011351682011-01-042000-09-290.0525104.56360.4016
DE00011414892011-04-082006-02-260.03599.75272.8593
DE00011351842011-07-042001-05-230.05104.37082.8962
DE00011414972011-10-142006-08-300.03599.60511.0519
DE00011351922012-01-042001-12-280.05104.86030.3825
DE00011415052012-04-132007-02-280.04101.34153.3661
DE00011352002012-07-042002-06-260.05105.292.8962
DE00011415132012-10-122007-08-240.0425102.49691.4631
DE00011352182013-01-042002-12-310.045103.76020.3443
DE00011352342013-07-042003-06-240.0375100.28032.1721
DE00011352422014-01-042003-10-210.0425102.60460.3251
DE00011352592014-07-042004-04-250.0425102.52912.4617
DE00011352672015-01-042004-10-270.037599.47480.2869
DE00011352832015-07-042005-04-280.032595.97021.8825
DE00011352912016-01-042005-10-300.03597.18150.2678
DE00011344682016-06-201986-06-200.06114.28493.7049
DE00011353092016-07-042006-04-260.04100.28472.3169
DE00011344922016-09-201986-09-200.05625112.232.0594
DE00011353172017-01-042006-10-310.037598.3970.2869
DE00011353332017-07-042007-04-270.0425102.02352.9262
DE00011353412018-01-042007-09-210.0499.84830.8415
DE00011349222024-01-041993-12-290.0625121.27110.4781
DE00011350442027-07-041997-07-030.065125.91573.765
DE00011350692028-01-041998-01-040.05625114.57910.4303
DE00011350852028-07-041998-10-070.0475103.22022.7514
DE00011351432030-01-042000-01-040.0625123.46680.4781
DE00011351762031-01-042000-10-270.055113.46940.4208
DE00011352262034-07-042003-01-220.0475103.18732.7514
DE00011352752037-01-042004-12-240.0491.56030.306
DE00011353252039-07-042006-12-280.042595.44414.3081
05-23 07:58