QuantLibPython/testHullWhiteBermudan.py at master · sschlenkrich/QuantLibPython · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import LinearLocator, FormatStrFormatter

import pandas

import QuantLib as ql
import QuantLibWrapper.YieldCurve as yc
from QuantLibWrapper.HullWhiteModel import HullWhiteModel, HullWhiteModelWithDiscreteNumeraire
from QuantLibWrapper.MCSimulation import MCSimulation
from QuantLibWrapper import Payoffs
from QuantLibWrapper.BermudanOption import BermudanOption
from QuantLibWrapper.DensityIntegrations import DensityIntegrationWithBreakEven, \
    SimpsonIntegration, HermiteIntegration, CubicSplineExactIntegration
from QuantLibWrapper.PDESolver import PDESolver
from QuantLibWrapper.AMCSolver import AMCSolver, AMCSolverOnlyExerciseRegression

# yield curves

flatCurve = yc.YieldCurve(['30y'],[0.03])

terms = [    '1y',    '2y',    '3y',    '4y',    '5y',    '6y',    '7y',    '8y',    '9y',   '10y',   '12y',   '15y',   '20y',   '25y',   '30y' ]
rates = [ 2.70e-2, 2.75e-2, 2.80e-2, 3.00e-2, 3.36e-2, 3.68e-2, 3.97e-2, 4.24e-2, 4.50e-2, 4.75e-2, 4.75e-2, 4.70e-2, 4.50e-2, 4.30e-2, 4.30e-2 ]
fwdRateYC = yc.YieldCurve(terms,rates)

# Hull-White model

meanReversion    = 0.05
volatilityTimes  = np.array([  1.0 , 2.0 , 5.0 , 10.0  ])
volatilityValues = np.array([  0.01, 0.01, 0.01,  0.01 ])

#hwModel          = HullWhiteModel(fwdRateYC,meanReversion,volatilityTimes,volatilityValues)
hwModel          = HullWhiteModelWithDiscreteNumeraire(fwdRateYC,meanReversion,volatilityTimes,volatilityValues)
#method           = SimpsonIntegration(hwModel,101,5)
#method           = DensityIntegrationWithBreakEven(SimpsonIntegration(hwModel,101,5))
#method           = HermiteIntegration(hwModel,10,101,5)
#method           = CubicSplineExactIntegration(hwModel,101,5)
#method           = DensityIntegrationWithBreakEven(CubicSplineExactIntegration(hwModel,101,5))
#method           = PDESolver(hwModel,11,2.0,0.5,1.0/12.0,0.0)
#method           = PDESolver(hwModel,11,2.0,0.5,1.0/12.0)

times  = np.linspace(0.0,19.0,20)
times  = np.array([0.0, 2.0, 6.0] + [ 12.0+k for k in range(8)])
nPaths = 10000
mcSim  = MCSimulation(hwModel,times,nPaths)
method = AMCSolver(mcSim,3,0.25)
#method = AMCSolverOnlyExerciseRegression(mcSim,3,0.25)

# now we test coupon bond opition pricing

# coupon bond option details

exercise  = 12.0
payTimes  = [ 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 20.0 ]
cashFlows = [ -1.0, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03,  1.0 ]

print('Bond option with zero strike: ' + str(hwModel.couponBondOption(exercise,payTimes,cashFlows,0.0,1.0)))
print('Bond option with unit strike: ' + str(hwModel.couponBondOption(exercise,payTimes[1:],cashFlows[1:],1.0,1.0)))

expiryTimes = np.array([ 2.0, 6.0, 12.0])
underlyings = [ Payoffs.Zero(), Payoffs.Zero(), Payoffs.CouponBond(hwModel,12.0,payTimes,cashFlows) ]
berm = BermudanOption(expiryTimes, underlyings, method)
print('Pseudo-Bermudan opt. (3 ex.): ' + str(berm.npv()))
payoff = Payoffs.Pay(Payoffs.VanillaOption(Payoffs.CouponBond(hwModel,12.0,payTimes,cashFlows),0.0,1.0),12.0)
print('MC bond option              : ' + str(mcSim.npv(payoff)))

#exit()
# test Bermudan bond option pricing

expiryTimes = np.array([ 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0])
underlyings = []
print('European bond option prices 12y - 19y...')
for k in range(8):
    pTimes = [payTimes[k]] + payTimes[(k+1):]
    cFlows = [-1.0 ] + cashFlows[(k+1):]
    print(hwModel.couponBondOption(payTimes[k],pTimes,cFlows,0.0,1.0))
    underlyings.append(Payoffs.CouponBond(hwModel,payTimes[k],pTimes,cFlows))
berm = BermudanOption(expiryTimes, underlyings, method)
print(berm.npv())