-
Notifications
You must be signed in to change notification settings - Fork 1
/
HW_run.py
146 lines (118 loc) · 6.21 KB
/
HW_run.py
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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
import numpy as np
import pandas as pd
from term_structure import *
def HW_theta(a: float, sigma: float, P0t: callable, eps: float) -> callable:
"""
Calculates the theta value for the Hull-White model
using a numeric approximation of the instantaneous forward rate
and the spot rate.
Args:
a (float): Mean reversion rate parameter a.
sigma (float): Volatility parameter sigma.
P0t (function handle): Function that calculates the price of a
zero-coupon bond as a function of time.
eps (float): Increment of time used in the numeric calculation of the
derivative of the instantaneous forward rate.
Returns:
theta (function): Function that returns the parameter theta of
Hull-White model at the time t implied by the calibration.
"""
def theta(t:float)->float:
insta_forward_term = (f0t(t+eps, P0t, eps)
-f0t(t-eps,P0t,eps))/(2.0*eps)
forward_term = a*f0t(t, P0t, eps)
variance_term = sigma**2/(2.0*a)*(1.0-np.exp(-2.0*a*t))
return insta_forward_term + forward_term + variance_term
return theta
def Paths(NoOfPaths, NoOfSteps, T, P0t, a, sigma, epsilon):
"""
Simulates a series of stochastic interest rate paths using the Hull-White model.
Args:
NoOfPaths (int): number of paths to simulate.
NoOfSteps (int): number of time steps per path.
T (float): end of the modelling window (in years).
(Ex. a modelling window of 50 years means T=50).
P0t (function): function that calculates the price of a
zero coupon bond issued at time 0 that matures at time t, with a
notional amount 1 and discounted using the assumed term structure.
a (float): mean reversion speed parameter a of
the Hull-White model.
sigma (float): volatility parameter sigma of the Hull-White model.
epsilon (float): size of the increment used for finite
difference approximation.
Returns:
dict: A dictionary containing arrays with time steps, interest rate paths,
and bond prices.
time (array): array of time steps.
R (array): array of interest rate paths with
shape (NoOfPaths, NoOfSteps+1).
M (array): array of bond prices with
shape (NoOfPaths, NoOfSteps+1).
Implemented by Gregor Fabjan from Open-Source Modelling on 29/07/2023.
Original inspiration: https://www.youtube.com/watch?v=BIZdwUDbnDo
"""
# Initial instantaneous forward rate at time t-> 0 (also spot rate at time 0).
# r(0) = f(0,0) = - partial derivative of log(P_mkt(0, epsilon) w.r.t epsilon)
r0 = f0t(epsilon, P0t, epsilon)
# Calculation of theta = 1/a * partial derivative of f(0,t) w.r.t. t
# + f(0,t) + sigma^2/(2 a^2)* (1-exp(-2*a*t)).
theta = HW_theta(a, sigma, P0t, epsilon)
# Generate the single source of random noise.
Z = np.random.normal(0.0, 1.0, [NoOfPaths, NoOfSteps])
# Initialize arrays
# Vector of time moments.
time = np.linspace(0, T, NoOfSteps+1)
W = np.zeros([NoOfPaths, NoOfSteps+1])
# Initialize array with interest rate increments
R = np.zeros([NoOfPaths, NoOfSteps+1])
# First interest rate equals the instantaneous forward (spot)
# rate at time 0.
R[:, 0] = r0
dt = T/float(NoOfSteps) # Size of increments between two steps
for iTime in range(1, NoOfSteps+1): # For each time increment
# Making sure the samples from the normal distribution have a mean of 0
# and variance 1
if NoOfPaths > 1:
Z[:, iTime-1] = (Z[:, iTime-1]-np.mean(Z[:, iTime-1]))/np.std(Z[:, iTime-1])
# Apply the Euler-Maruyama discretisation scheme for the Hull-White model
# at each time increment.
W[:, iTime] = W[:, iTime-1] + np.power(dt, 0.5)*Z[:, iTime-1]
noise_term = sigma* (W[:, iTime]-W[:, iTime-1])
rate_term = (theta(time[iTime-1])-a*R[:, iTime-1])*dt
R[:, iTime] = R[:, iTime-1] + rate_term + noise_term
# Vectorized numeric integration using the Euler integration method .
M = np.exp(-0.5 * (R[:, :-1] + R[:, 1:]) * dt)
M = np.insert(M, 0, 1, axis=1).cumprod(axis=1)
# Output is a dataframe with time moment, the interest rate path and the price
# of a zero coupon bond issued at time 0 that matures at the selected time
# moment with a notional value of 1.
paths = {"time":time, "R":R, "M":M}
return paths
def mainCalculation(NoOfPaths, NoOfSteps, T, a, sigma, P0t, epsilon):
"""
Calculates and plots the prices of zero-coupon bonds (ZCB) calculated
using the Hull-White model`s analytical formula and the Monte Carlo simulation.
Args:
NoOfPaths (int): number of Monte Carlo simulation paths.
NoOfSteps (int): number of time steps per path.
T (float): length in years of the modelling window (Ex. 50 years means t=50).
a (float): mean reversion rate parameter a of the Hull-White model.
sigma (float): volatility parameter sigma of the Hull-White model.
P0t (function): function that calculates the price of a zero coupon bond issued.
at time 0 that matures at time t, with a notional amount 1 and discounted using
the assumed term structure.
epsilon (float): the size of the increment used for finite difference approximation.
Returns:
t XXX: time increments.
P XXX: average of the sumulated paths.
implied_term_structure XXX: term structure provided as input into the HW simulation.
"""
paths = Paths(NoOfPaths, NoOfSteps, T, P0t, a, sigma, epsilon)
M = paths["M"]
t = paths["time"]
implied_term_structure = P0t(t)
# Compare the price of an option on a ZCB from Monte Carlo and the analytical expression
P = np.zeros([NoOfSteps+1])
for i in range(0, NoOfSteps+1):
P[i] = np.mean(M[:, i])
return [t, P, implied_term_structure, M]