Basic principles of pricing and hedging¶
We shall first go through basic principles of modeling in mathematical Finance before we introduce concepts of interest rate theory.
We have seen in one step bi- and tri-nomial models that pricing and hedging are in a fundamental duality relationship given a certain payoff: the largest arbitrage free price equals the smallest price of a super-hedging portfolio.
Prices are free of arbitrage if their introduction into the market does not lead to arbitrages.
A super-hedging portfolio is a self-financing portfolio (i.e. the value process is the initial value of the portfolio plus the P&L process -- all in discounted terms) dominating a certain payoff.
In case that there is only one arbitrage-free price for a payoff $X$ the super-hedging portfolio actually is a hedging portfolio leading to the remarkable relationship (in discounted terms) $$ X = \text{ price } + \text{ P & L } \, . $$ Furthermore models are free of arbitrage if and only if there exists an equivalent pricing measure, which associates in particular to P & L processes the value $0$. Whence prices of payoffs can be calculated by taking expectations (i.e. a Monte Carlo evaluation is possible) with respect to this equivalent pricing measure.
Next we shall see the concepts of pricing and hedging realized in a geometric Brownian motion market environment. We shall use code from https://github.com/yhilpisch/dx for this purpose.
import dx
import datetime as dt
import pandas as pd
import seaborn as sns; sns.set()
import numpy as np
First we shall define some market environment:
r = dx.constant_short_rate('r', 0.01)
# a constant short rate
cas_1 = dx.market_environment('cas', dt.datetime(2016, 1, 1))
cas_1.add_constant('initial_value', 100.)
# starting value of simulated processes
cas_1.add_constant('volatility', 0.2)
# volatiltiy factor
cas_1.add_constant('final_date', dt.datetime(2017, 1, 1))
# horizon for simulation
cas_1.add_constant('currency', 'EUR')
# currency of instrument
cas_1.add_constant('frequency', 'D')
# frequency for discretization
cas_1.add_constant('paths', 10000)
# number of paths
cas_1.add_curve('discount_curve', r)
# number of paths
Let us introduce a geometric Brownian motion in the above market environment.
gbm_cas_1 = dx.geometric_brownian_motion('gbm_1', cas_1)
We can obtain one trajectory generated on the predefined weekly time grid of GBM.
paths_gbm_cas_1 = pd.DataFrame(gbm_cas_1.get_instrument_values(), index=gbm_cas_1.time_grid)
%matplotlib inline
paths_gbm_cas_1.loc[:, :10].plot(legend=False, figsize=(10, 6))
Next we define a new market enviroment which inherits from the previous one but has an additional Europen call option.
cas_1_opt = dx.market_environment('cas_1_opt', cas_1.pricing_date)
cas_1_opt.add_environment(cas_1)
cas_1_opt.add_constant('maturity', dt.datetime(2017, 1, 1))
cas_1_opt.add_constant('strike', 110.)
eur_call = dx.valuation_mcs_european_single(
name='eur_call',
underlying=gbm_cas_1,
mar_env=cas_1_opt,
payoff_func='np.maximum(maturity_value - strike, 0)')
The present value of the European call is uniquely given by no arbitrage arguments.
eur_call.present_value()
The delta is the sensitivity with respect to the initial value.
eur_call.delta()
The gamma is the sensitivity of the delta with respect to the initial value.
eur_call.gamma()
The vega is the sensitivity with respect to volatility.
eur_call.vega()
The theta is sensitivity with respect to maturity.
eur_call.theta()
The rho is sensitivity with respect to interest rate.
eur_call.rho()
The previous quantities allow to understand the risks of the given derivative from the point of view of the chosen model in terms of sensitivities.
In the sequel we demonstrate the precise meaning of the option's delta by means of a running hedging portfolio.
path = gbm_cas_1.get_instrument_values()[:,0]
timegrid = gbm_cas_1.time_grid
presentvalue = eur_call.present_value()
n = len(path)
pnl = [presentvalue]
optionvalue = [np.maximum(0,path[0]-eur_call.strike)]
for i in range(n-1):
r = dx.constant_short_rate('r', 0.01)
# a constant short rate
running = dx.market_environment('running', timegrid[i])
running.add_constant('initial_value', path[i])
# starting value of simulated processes
running.add_constant('volatility', 0.2)
# volatiltiy factor
running.add_constant('final_date', dt.datetime(2017, 1, 1))
# horizon for simulation
running.add_constant('currency', 'EUR')
# currency of instrument
running.add_constant('frequency', 'W')
# frequency for discretization
running.add_constant('paths', 10000)
# number of paths
running.add_curve('discount_curve', r)
# number of paths
gbm_running = dx.geometric_brownian_motion('gbm_running', running)
opt_running = dx.market_environment('opt_running', running.pricing_date)
opt_running.add_environment(running)
opt_running.add_constant('maturity', dt.datetime(2017, 1, 1))
opt_running.add_constant('strike', 110.)
eur_call = dx.valuation_mcs_european_single(
name='eur_call',
underlying=gbm_running,
mar_env=opt_running,
payoff_func='np.maximum(maturity_value - strike, 0)')
#print(path[i])
#print(timegrid[i])
#print(eur_call.delta())
pnl = pnl + [pnl[-1]+eur_call.delta()*(path[i+1]-path[i])]
optionvalue = optionvalue + [np.maximum(0,path[i+1]-eur_call.strike)]
data = []
for j in range(n):
data = data + [[optionvalue[j],pnl[j]]]
paths_hedge = pd.DataFrame(data, index=gbm_cas_1.time_grid)
%matplotlib inline
paths_hedge.loc[:,:].plot(legend=True, figsize=(10, 6))
Interest rate terminology¶
Prices of zero-coupon bonds (ZCB) with maturity $ T $ are denoted by $P(t,T)$. Interest rates are governed by a market of (default free) zero-coupon bonds modeled by stochastic processes $ {(P(t,T))}_{0 \leq t \leq T} $ for $ T \geq 0 $. We assume the normalization $ P(T,T)=1 $.
$T$ denotes the maturity of the bond, $P(t,T)$ its price at a time $t$ before maturity $T$.
The yield $$ Y(t,T) = - \frac{1}{T-t} \log P(t,T) $$ describes the compound interest rate p.a. for maturity $T$.
The curve $f$ is called the forward rate curve of the bond market \begin{align*} P(t,T) & =\exp(-\int_{t}^{T}f(t,s)ds) \end{align*} for $0\leq t\leq T$.
Let us understand a bit the dynamics of yield curves and forward rate curves at this point:
from mpl_toolkits.mplot3d import Axes3D
import copy as copylib
from progressbar import *
%pylab
%matplotlib inline
import pandas as pandas
pylab.rcParams['figure.figsize'] = (16, 4.5)
numpy.random.seed(0)
First we load some historical data:
dataframe = pandas.DataFrame.from_csv('hjm_data.csv')
dataframe = dataframe / 100 # Convert interest rates to %
pandas.options.display.max_rows = 10
display(dataframe)
hist_timeline = list(dataframe.index)
tenors = [float(x) for x in dataframe.columns]
hist_rates = matrix(dataframe)
plot(hist_rates), xlabel(r'Time $t$'), title(r'Historical $f(t,\tau)$ by $t$'), show()
plot(tenors, hist_rates.transpose()), xlabel(r'Tenor $\tau$'), title(r'Historical $f(t,\tau)$ by $\tau$');
In a next step we look at returns for a certain maturity:
diff_rates = diff(hist_rates, axis=0)
assert(hist_rates.shape[1]==diff_rates.shape[1])
plot(diff_rates), xlabel(r'Time $t$'), title(r'$df(t,\tau)$ by $t$');
This demonstrates impressively the complexity and beauty of interest rate behaviour.
The short rate process is given through $R_{t}=f(t,t)$ for $t\geq0$ defining the ``bank acco unt process'' $$ (B(t))_{t \geq 0}:=(\exp(\int_{0}^{t}R_{s}ds))_{t \geq 0}. $$
No arbitrage is guaranteed if on the filtered probability space $ (\Omega,\mathcal{F},Q) $ with filtration $ {(\mathcal{F}_t)}_{t \geq 0} $, $$ E(\exp(-\int_t^T R_s ds)|\mathcal{F}_t) = P(t,T) $$ holds true for $ 0 \leq t \leq T $ for some equivalent (martingale) measure $ P $.
Consider a bond market $(P(t,T))_{t \leq T}$ with $P(T,T) = 1$ and $P(t,T) > 0$. Let $t \leq T \leq T^{\ast}$. We define the simple forward rate through \begin{align*} F(t;T,T^{\ast}) := \frac{1}{T^{\ast} - T} \bigg( \frac{P(t,T)}{P(t,T^{\ast})} - 1 \bigg). \end{align*} and the simple spot rate through \begin{align*} F(t,T) := F(t;t,T). \end{align*}
Apparently $P(t,T^{\ast}) F(t;T,T^{\ast})$ is the fair value at time $t$ of a contract paying $F(T,T^{\ast})$ at time $T^{\ast}$.
Indeed, note that \begin{align*} P(t,T^{\ast}) F(t;T,T^{\ast}) &= \frac{P(t,T) - P(t,T^{\ast})}{T^{\ast} - T}, \\ F(T,T^{\ast}) &= \frac{1}{T^{\ast} - T} \bigg( \frac{1}{P(T,T^{\ast})} - 1 \bigg). \end{align*} Fair value means that we can build a portfolio at time $ t $ and at price $ \frac{P(t,T) - P(t,T^{\ast})}{T^{\ast} - T} $ yielding $ F(T,T^{\ast}) $ at time $ T^{\ast} $:
-> Holding a ZCB with maturity $ T $ at time $ t $ has value $ P(t,T) $, being additionally short in a ZCB with maturity $ T^{\ast} $ amounts all together to $P(t,T) - P(t,T^{\ast})$.
-> at time $T$ we have to rebalance the portfolio by buying with the maturing ZCB another bond with maturity $ T^{\ast} $, precisely an amount $ 1/P(T,T^{\ast}) $.
-> at time $ T^{\ast} $ we receive $ 1/P(T,T^{\ast}) - 1 $.
As an example we introduce a particular short rate model, the so called Cox Ingersoll Ross model or squared diffusion model. We define a second market environment and calculate zero coupon bond prices with respect to this model:
cas_2 = dx.market_environment(name='cas_2', pricing_date=dt.datetime(2015, 1, 1))
cas_2.add_constant('initial_value', 0.01)
cas_2.add_constant('volatility', 0.1)
cas_2.add_constant('kappa', 2.0)
cas_2.add_constant('theta', 0.05)
cas_2.add_constant('paths', 10000)
cas_2.add_constant('frequency', 'W')
cas_2.add_constant('starting_date', cas_2.pricing_date)
cas_2.add_constant('final_date', dt.datetime(2015, 12, 31))
We next define our model:
ssr = dx.stochastic_short_rate('sr', cas_2)
time_list = ssr.process.time_grid
forward_rates=ssr.get_forward_rates(time_list, 10)
ZCB = ssr.get_discount_factors(time_list, 10)
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
%matplotlib inline
# short rate paths
plt.figure(figsize=(10, 6))
plt.plot(ssr.process.get_instrument_values()[:, :10]);
reversed_time_list = []
m = len(time_list)
for i in range(m):
reversed_time_list= reversed_time_list + [time_list[m-1-i]]
ZCB_df = pd.DataFrame(ZCB[1], index=reversed_time_list)
%matplotlib inline
ZCB_df.loc[:,:].plot(legend=True, figsize=(10, 6))
Caps¶
In the sequel, we fix a number of future dates \begin{align*} T_0 < T_1 < \ldots < T_n \end{align*} with $T_i - T_{i-1} \equiv \delta$.
Fix a rate $\kappa > 0$. At time $T_i$ the holder of the cap receives \begin{align*} \delta (F(T_{i-1},T_i) - \kappa)^+. \end{align*} Let $t \leq T_0$. We write \begin{align*} {\rm Cpl}(t;T_{i-1},T_i), \quad i=1,\ldots,n \end{align*} for the time $t$ price of the $i$th caplet, and \begin{align*} {\rm Cp}(t) = \sum_{i=1}^n {\rm Cpl}(t;T_{i-1},T_i) \end{align*} for the time $t$ price of the cap.
Floors¶
At time $T_i$ the holder of the floor receives \begin{align*} \delta (\kappa - F(T_{i-1},T_i))^+. \end{align*} Let $t \leq T_0$. We write \begin{align*} {\rm Fll}(t;T_{i-1},T_i), \quad i=1,\ldots,n \end{align*} for the time $t$ price of the $i$th floorlet, and \begin{align*} {\rm Fl}(t) = \sum_{i=1}^n {\rm Fll}(t;T_{i-1},T_i) \end{align*} for the time $t$ price of the floor.
Swaps¶
Fix a rate $K$ and a nominal $N$. The cash flow of a payer swap at $T_i$ is \begin{align*} (F(T_{i-1},T_i) - K) \delta N. \end{align*} The total value $\Pi_p(t)$ of the payer swap at time $t \leq T_0$ is \begin{align*} \Pi_p(t) = N \bigg( P(t,T_0) - P(t,T_n) - K \delta \sum_{i=1}^n P(t,T_i) \bigg). \end{align*} The value of a receiver swap at $t \leq T_0$ is \begin{align*} \Pi_r(t) = -\Pi_p(t). \end{align*} The swap rate $R_{\rm swap}(t)$ is the fixed rate $K$ which gives $\Pi_p(t) = \Pi_r(t) = 0$. Hence \begin{align*} R_{\rm swap}(t) = \frac{P(t,T_0) - P(t,T_n)}{\delta \sum_{i=1}^n P(t,T_i)}, \quad t \in [0,T_0]. \end{align*}
Swaptions¶
A \emph{payer (receiver) swaption} is an option to enter a payer (receiver) swap at $T_0$. The payoff of a payer swaption at $T_0$ is \begin{align*} N \delta (R_{\rm swap}(T_0) - K)^+ \sum_{i=1}^n P(T_0,T_i), \end{align*} and of a receiver swaption \begin{align*} N \delta (K - R_{\rm swap}(T_0))^+ \sum_{i=1}^n P(T_0,T_i). \end{align*}
Spot measure¶
From now on, let $P$ be a martingale measure in the bond market $(P(t,T))_{t \leq T}$, i.e. for each $T > 0$ the discounted bond price process \begin{align*} \frac{P(t,T)}{B(t)} \end{align*} is a martingale. This leads to the following fundamental formula of interest rate theory $$ P(t,T) = E(\exp(-\int_t^T R_s ds)) | \mathcal{F}_t ) $$ for $ 0 \leq t \leq T $.
Forward measures¶
For $T^{\ast} > 0$ define the \emph{$T^{\ast}$-forward measure} ${P}^{T^{\ast}} $ such that for any $T > 0$ the discounted bond price process \begin{align*} \frac{P(t,T)}{P(t,T^{\ast})}, \quad t \in [0,T] \end{align*} is a ${P}^{T^{\ast}}$-martingale.
For any $T < T^{\ast}$ the simple forward rate \begin{align*} F(t;T,T^{\ast}) = \frac{1}{T^{\ast} - T} \bigg( \frac{P(t,T)}{P(t,T^{\ast})} - 1 \bigg) \end{align*} is a $\mathbb{P}^{T^{\ast}}$-martingale.
For any time derivative $X \in \mathcal{F}_{T^*}$ paid at $ T^* $ we have that the fair value via ``martingale pricing'' is given through \begin{align*} P(t,T^*) \mathbb{E}^{T^*}[X | \mathcal{F}_t]. \end{align*} The fair price of the $i$th caplet is therefore given by \begin{align*} {\rm Cpl}(t;T_{i-1},T_i) = \delta P(t,T_i) \mathbb{E}^{T_i} [ (F(T_{i-1},T_i) - \kappa)^+ | \mathcal{F}_t]. \end{align*} By the martingale property we obtain therefore \begin{align*} \mathbb{E}^{T_i}[F(T_{i-1},T_i) | \mathcal{F}_t] = F(t;T_{i-1},T_i), \end{align*} what was proved by trading arguments before.
Black's formula¶
Let $X \sim N(\mu,\sigma^2)$ and $K \in \mathbb{R}$. Then we have \begin{align*} \mathbb{E}[ (e^X - K)^+ ] &= e^{\mu + \frac{\sigma^2}{2}} \Phi \bigg( -\frac{\log K - (\mu + \sigma^2)}{\sigma} \bigg)
- K \Phi \bigg( -\frac{\log K - \mu}{\sigma} \bigg), \ \mathbb{E}[ (K - e^X)^+ ] &= K \Phi \bigg( \frac{\log K - \mu}{\sigma} \bigg) - e^{\mu + \frac{\sigma^2}{2}} \Phi \bigg( \frac{\log K - (\mu + \sigma^2)}{\sigma} \bigg). \end{align*}
Black's formula for caps and floors¶
Let $t \leq T_0$. From our previous results we know that \begin{align*} {\rm Cpl}(t;T_{i-1},T_i) &= \delta P(t,T_i) \mathbb{E}_t^{{T_i}} [ (F(T_{i-1},T_i) - \kappa)^+ ], \\ {\rm Fll}(t;T_{i-1},T_i) &= \delta P(t,T_i) \mathbb{E}_t^{{T_i}} [ (\kappa - F(T_{i-1},T_i))^+ ], \end{align*} and that $F(t;T_{i-1},T_i)$ is an $P^{T_i}$-martingale.
We assume that under ${P}^{T_i}$ the forward rate $F(t;T_{i-1},T_i)$ is an exponential Brownian motion \begin{align*} F(t;T_{i-1},T_i) & = F(s;T_{i-1},T_i) \\ & \exp \bigg( -\frac{1}{2} \int_s^t \lambda(u,T_{i-1})^2 du + \int_s^t \lambda(u,T_{i-1}) dW_u^{T_i} \bigg) \end{align*} for $s \leq t \leq T_{i-1}$, with a function $\lambda(u,T_{i-1})$.
We define the variance $\sigma^2(t)$ as \begin{align*} \sigma^2(t) := \frac{1}{T_{i-1} - t} \int_t^{T_{i-1}} \lambda(s,T_{i-1})^2 ds. \end{align*} The $P^{T_i}$-distribution of $\log F(T_{i-1},T_i)$ conditional on $\mathcal{F}_t$ is $N(\mu,\sigma^2)$ with \begin{align*} \mu &= \log F(t;T_{i-1},T_i) - \frac{\sigma^2(t)}{2} (T_{i-1} - t), \\ \sigma^2 &= \sigma^2(t) (T_{i-1} - t). \end{align*} In particular \begin{align*} \mu + \frac{\sigma^2}{2} &= \log F(t;T_{i-1},T_i), \\ \mu + \sigma^2 &= \log F(t;T_{i-1},T_i) + \frac{\sigma^2(t)}{2} (T_{i-1} - t). \end{align*}
We have \begin{align*} {\rm Cpl}(t;T_{i-1},T_i) &= \delta P(t,T_i) (F(t;T_{i-1},T_i) \Phi(d_1(i;t)) - \kappa \Phi(d_2(i;t))), \\ {\rm Fll}(t;T_{i-1},T_i) &= \delta P(t,T_i) ( \kappa \Phi(-d_2(i;t)) - F(t;T_{i-1},T_i) \Phi(-d_1(i;t)) ), \end{align*} where \begin{align*} d_{1,2}(i;t) = \frac{\log \big( \frac{F(t;T_{i-1},T_i)}{\kappa} \big) \pm \frac{1}{2} \sigma(t)^2 (T_{i-1} - t)}{\sigma(t) \sqrt{T_{i-1} - t}}. \end{align*}
We shall now use Quantlib library to valuate a cap in a constant and non-constant volatility market environment. We use here code kindly provided at http://gouthamanbalaraman.com/blog/interest-rate-cap-floor-valuation-quantlib-python.html. First we have to specify the market environment, then we specify the cap as a sum of caplets.
import QuantLib as ql
calc_date = ql.Date(14, 6, 2016)
ql.Settings.instance().evaluationDate = calc_date
dates = [ql.Date(14,6,2016), ql.Date(14,9,2016),
ql.Date(14,12,2016), ql.Date(14,6,2017),
ql.Date(14,6,2019), ql.Date(14,6,2021),
ql.Date(15,6,2026), ql.Date(16,6,2031),
ql.Date(16,6,2036), ql.Date(14,6,2046)
]
yields = [0.000000, 0.006616, 0.007049, 0.007795,
0.009599, 0.011203, 0.015068, 0.017583,
0.018998, 0.020080]
day_count = ql.ActualActual()
calendar = ql.UnitedStates()
interpolation = ql.Linear()
compounding = ql.Compounded
compounding_frequency = ql.Annual
term_structure = ql.ZeroCurve(dates, yields, day_count, calendar,
interpolation, compounding, compounding_frequency)
ts_handle = ql.YieldTermStructureHandle(term_structure)
start_date = ql.Date(14, 6, 2016)
end_date = ql.Date(14, 6 , 2026)
period = ql.Period(3, ql.Months)
calendar = ql.UnitedStates()
buss_convention = ql.ModifiedFollowing
rule = ql.DateGeneration.Forward
end_of_month = False
schedule = ql.Schedule(start_date, end_date, period,
calendar, buss_convention, buss_convention,
rule, end_of_month)
ibor_index = ql.USDLibor(ql.Period(3, ql.Months), ts_handle)
ibor_index.addFixing(ql.Date(10,6,2016), 0.0065560)
ibor_leg = ql.IborLeg([1000000], schedule, ibor_index)
Now all required pieces are set and we can evaluate the cap:
strike = 0.02
cap = ql.Cap(ibor_leg, [strike])
vols = ql.QuoteHandle(ql.SimpleQuote(0.547295))
engine = ql.BlackCapFloorEngine(ts_handle, vols)
cap.setPricingEngine(engine)
print(cap.NPV())
Reality, however, is more complicated, and whence volatility not constant.
strikes = [0.01,0.015, 0.02]
expiries = [ql.Period(i, ql.Years) for i in range(1,11)] + [ql.Period(12, ql.Years)]
vols = ql.Matrix(len(expiries), len(strikes))
data = [[47.27, 55.47, 64.07, 70.14, 72.13, 69.41, 72.15, 67.28, 66.08, 68.64, 65.83],
[46.65,54.15,61.47,65.53,66.28,62.83,64.42,60.05,58.71,60.35,55.91],
[46.6,52.65,59.32,62.05,62.0,58.09,59.03,55.0,53.59,54.74,49.54]
]
for i in range(vols.rows()):
for j in range(vols.columns()):
vols[i][j] = data[j][i]/100.0
Cap-Floor volatilities are provided by the market, Quantlib strips volatilities out of it.
calendar = ql.UnitedStates()
bdc = ql.ModifiedFollowing
daycount = ql.Actual365Fixed()
settlement_days = 2
capfloor_vol = ql.CapFloorTermVolSurface(settlement_days, calendar, bdc, expiries, strikes, vols, daycount)
optionlet_surf = ql.OptionletStripper1(capfloor_vol, ibor_index)
ovs_handle = ql.OptionletVolatilityStructureHandle(
ql.StrippedOptionletAdapter(optionlet_surf)
)
The procedure is visualized below:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
tenors = np.arange(0,10,0.25)
strike = 0.015
capfloor_vols = [capfloor_vol.volatility(t, strike) for t in tenors]
opionlet_vols = [ovs_handle.volatility(t, strike) for t in tenors]
plt.plot(tenors, capfloor_vols, "--", label="CapFloor Vols")
plt.plot(tenors, opionlet_vols,"-", label="Optionlet Vols")
plt.legend(bbox_to_anchor=(0.5, 0.25))
The resulting cap price is slightly different from before
engine2 = ql.BlackCapFloorEngine(ts_handle, ovs_handle)
cap.setPricingEngine(engine2)
print(cap.NPV())
The HJM approach and another way of cap pricing¶
We consider the yield curve time series from above and build a model from it directly. We first calculate the empirical covariance matrix, annualize it and look at its principal components, i.e. the eigenvectors corresponding to the largest eigenvalues (here three).
sigma = cov(diff_rates.transpose())
print("Sigma shape : " + str(sigma.shape))
sigma *= 252
eigval, eigvec = linalg.eig(sigma)
eigvec=matrix(eigvec)
assert type(eigval) == ndarray
assert type(eigvec) == matrix
print(eigval)
factors=3
index_eigvec = list(reversed(eigval.argsort()))[0:factors] # highest principal component first in the array
princ_eigval = array([eigval[i] for i in index_eigvec])
princ_comp = hstack([eigvec[:,i] for i in index_eigvec])
print("Principal eigenvalues")
print(princ_eigval)
print()
print("Principal eigenvectors")
print(princ_comp)
plot(princ_comp, marker='.'), title('Principal components'), xlabel(r'Tenor $\tau$');
Next we calculate the volatilities as a function of time: volatilities are just principal component multiplied by the overall volatility - there are different choice what to do here.
sqrt_eigval = matrix(princ_eigval ** .5)
tmp_m = vstack([sqrt_eigval for i in range(princ_comp.shape[0])]) # resize matrix (1,factors) to (n, factors)
vols = multiply(tmp_m, princ_comp) # multiply matrice element-wise
print('vols shape: ' + str(vols.shape))
plot(vols, marker='.'), xlabel(r'Tenor $\tau$'), ylabel(r'Volatility $\sigma$'), title('Discretized volatilities');
The discretized volatalities have to be made into curves for building the model, so we interpolate (by cubic splines), the result is drawn below.
def get_matrix_column(mat, i):
return array(mat[:,i].flatten())[0]
class PolynomialInterpolator:
def __init__(self, params):
assert type(params) == ndarray
self.params = params
def calc(self, x):
n = len(self.params)
C = self.params
X = array([x**i for i in reversed(range(n))])
return sum(multiply(X, C))
fitted_vols = []
def fit_volatility(i, degree, title):
vol = get_matrix_column(vols, i)
fitted_vol = PolynomialInterpolator(polyfit(tenors, vol, degree))
plot(tenors, vol, marker='.', label='Discretized volatility')
plot(tenors, [fitted_vol.calc(x) for x in tenors], label='Fitted volatility')
plt.title(title), xlabel(r'Time $t$'), legend();
fitted_vols.append(fitted_vol)
subplot(1, 3, 1), fit_volatility(0, 0, '1st component');
subplot(1, 3, 2), fit_volatility(1, 3, '2nd component');
subplot(1, 3, 3), fit_volatility(2, 3, '3rd component');
We prepare now for solving the HJM equation by simulation: we choose as a dynamics time dependent volatilities for each principal components. First we set up an integration routine (trapezoidal rule) and we discretize the interpolated volatilities.
def integrate(f, x0, x1, dx):
n = (x1-x0)/dx+1
out = 0
for i, x in enumerate(linspace(x0, x1, n)):
if i==0 or i==n-1:
out += 0.5 * f(x)
else:
out += f(x) # not adjusted by *0.5 because of repeating terms x1...xn-1 - see trapezoidal rule
out *= dx
return out
mc_tenors = linspace(0,25,51)
# Discretize fitted volfuncs for the purpose of monte carlo simulation
mc_vols = matrix([[fitted_vol.calc(tenor) for tenor in mc_tenors] for fitted_vol in fitted_vols]).transpose()
plot(mc_tenors, mc_vols, marker='.'), xlabel(r'Time $t$'), title('Volatilities');
def m(tau, fitted_vols):
#This funciton carries out integration for all principal factors.
#It uses the fact that volatility is function of time in HJM model
out = 0.
for fitted_vol in fitted_vols:
assert isinstance(fitted_vol, PolynomialInterpolator)
out += integrate(fitted_vol.calc, 0, tau, 0.01) * fitted_vol.calc(tau)
return out
mc_drift = array([m(tau, fitted_vols) for tau in mc_tenors])
plot(mc_drift, marker='.'), xlabel(r'Time $t$'), title('Risk-neutral drift');
curve_spot = array(hist_rates[-1,:].flatten())[0]
plot(mc_tenors, curve_spot.transpose(), marker='.'), ylabel('$f(t_0,T)$'), xlabel("$T$");
def simulation(f, tenors, drift, vols, timeline):
assert type(tenors)==ndarray
assert type(f)==ndarray
assert type(drift)==ndarray
assert type(timeline)==ndarray
assert len(f)==len(tenors)
vols = array(vols.transpose()) # 3 rows, T columns
len_tenors = len(tenors)
len_vols = len(vols)
yield timeline[0], copylib.copy(f)
for it in range(1, len(timeline)):
t = timeline[it]
dt = t - timeline[it-1]
sqrt_dt = sqrt(dt)
fprev = f
f = copylib.copy(f)
random_numbers = [normal() for i in range(len_vols)]
for iT in range(len_tenors):
val = fprev[iT] + drift[iT] * dt
#
sum = 0
for iVol, vol in enumerate(vols):
sum += vol[iT] * random_numbers[iVol]
val += sum * sqrt_dt
#
iT1 = iT+1 if iT<len_tenors-1 else iT-1 # if we can't take right difference, take left difference
dfdT = (fprev[iT1] - fprev[iT]) / (iT1 - iT)
val += dfdT * dt
#
f[iT] = val
yield t,f
proj_rates = []
proj_timeline = linspace(0,5,500)
progressbar = ProgressBar("One simulation path", len(proj_timeline))
for i, (t, f) in enumerate(simulation(curve_spot, mc_tenors, mc_drift, mc_vols, proj_timeline)):
progressbar.update(i)
proj_rates.append(f)
proj_rates = matrix(proj_rates)
plot(proj_timeline.transpose(), proj_rates), xlabel(r'Time $t$'), ylabel(r'Rate $f(t,\tau)$');
title(r'Simulated $f(t,\tau)$ by $t$'), show()
plot(mc_tenors, proj_rates.transpose()), xlabel(r'Tenor $\tau$'), ylabel(r'Rate $f(t,\tau)$');
title(r'Simulated $f(t,\tau)$ by $\tau$'), show()
We can now price a simple caplet with this model: strike $K=3\%$, expiring in $t_s=1.0$, maturing in $t_e=2.0$, notional $N=1,000,000$.
class Integrator:
def __init__(self, x0, x1):
assert x0 < x1
self.sum, self.n, self.x0, self.x1= 0, 0, x0, x1
def add(self, value):
self.sum += value
self.n += 1
def get_integral(self):
return (self.x1 - self.x0) * self.sum / self.n
t_exp = 1.
t_mat = 2.
K = .03
notional = 1e6
n_simulations = 500
n_timesteps = 50
proj_timeline = linspace(0,t_mat, n_timesteps)
progressbar = ProgressBar('Performing simulation', n_simulations)
simulated_forecast_rates = []
simulated_df = []
simulated_pvs = []
pv_convergence_process = []
for i in range(0, n_simulations):
progressbar.update(i)
rate_forecast = None # Forecast rate between t_exp and t_mat for this path
rate_discount = Integrator(0, t_exp) # cont.compounded discount rate for this path
for t, curve_fwd in simulation(curve_spot, mc_tenors, mc_drift, mc_vols, proj_timeline):
f_t_0 = interp(0., mc_tenors, curve_fwd) # rate $f_t^0$
rate_discount.add(f_t_0)
if t>=t_exp and rate_forecast is None: # t is expiration time
Tau = t_mat-t_exp
rate_forecast = Integrator(0, Tau) # integrate all inst.fwd.rates from 0 till 1Y tenor to get 1Y spot rate
for s in linspace(0, Tau, 15): # $\int_0^T f(t,s)ds$
f_texp_s = interp(s, mc_tenors, curve_fwd)
rate_forecast.add(f_texp_s)
rate_forecast = rate_forecast.get_integral()
plot(mc_tenors, curve_fwd), xlabel(r'Tenor $\tau$'), ylabel(r'Rate $f(t_{exp},\tau)$'); # Plot forward curve at expiration
df = exp(-rate_discount.get_integral()) # Discount factor
simulated_forecast_rates.append(rate_forecast)
simulated_df.append(df)
pv = max(0, rate_forecast - K) * (t_mat-t_exp) * notional * df
simulated_pvs.append(pv)
pv_convergence_process.append(average(simulated_pvs))
show()
#
scatter(simulated_df, simulated_forecast_rates), xlabel('Discount Factor'), ylabel('Forecast Rate')
show()
#
plot(pv_convergence_process[10:]), title('Convergence of PV'), xlabel("Simulations"), ylabel("V");
print("Final value: %f" % pv_convergence_process[-1])
