p2

#Introduction
ATM's: Affine Term Structure Models

Hull-White Analytical Model

The Hull–White (1990) general extensions of Vasicek can be fitted to the initial yield and volatility curve. However, this flexibility has its price: the model cannot be handled analytically in general. We therefore restrict to the following extension of the Vasicek model that was analyzed by Hull and White in 1994 (analytic HW model):

\[dr(t) = (b(t)+ar(t))dt + \sigma dW(t)\]

\(a\) and \(\sigma\) are chosen to fit volatility structure whereas \(b(t)\) is chosen in order to match the initial yield curve.
The model can be represented also as:
\[ r(t)=x(t)+\alpha(t)\]
where
\[dx(t)=-ax(t)dt+\sigma dW(t)\]
\[x(0)=0\]
\[\alpha(t)=fwd(0,t)+\frac{\sigma^2}{2a^2}\left(1-e^{-at}\right)^2.\]
\(fwd(0,t)\) is the market-implied istanteneous forward rate at maturity \(t>0\).

Analytical HW provide an affine term structure because the bond prices are of the form

\[ P(t,T) = A(t,T)\exp{(B(t,T)r(t))}\]

with
\[B(t,T)=\frac{1}{a}(e^{-a(T-t)}-1)\]
and
\[A(t,T)=\frac{P(0,T)}{P(0,t)}\exp{\left(-B(t,T)fwd(0,t)-\frac{1}{2}B^2(t,T)\phi(t)\right)}\]
where
\[
\phi(t)=\frac{\sigma^2}{2a}\left(1-\exp(-2at)\right)
\]
We are interested in the continuously compounded spot rate \(R(t,T)=\frac{-\log(P(t,T))}{T-t}\) and, in the case of analytic HW, we can obtain it explicitly:

\[\begin{aligned}
R(t,T) &= \
&=\frac{-\log(P(t,T))}{T-t}\
&= \frac{1}{T-t}\left(-\log(\frac{P(0,T)}{P(0,t)}) +\frac{B(t,T)^2}{2}\phi(t)+B(t,T)fwd(0,t)-B(t,T)r(t) \right)\
&= \frac{1}{T-t}\left(B(t,T)\left(-r(t)+fwd(0,t)+B(t,T)\frac{\phi(t)}{2}\right)+
R(0,T)T – R(0,t)t \right)
\end{aligned}
\]
For a more formal derivation you can see Holmgaard and Lund

We want to simulate 200 scenarios of \(R(t,T)\) every month, starting from next month to 20 years, with maturities every 6 month for 20 years from a specific yield curve.

#### Parameters' settings ####
# Frequency of simulation and interpolation
freq <- "monthly"
delta_t <- 1/12
horizon <- 10 # I take horizon = 10 years
nb.sims <- 1000 # Scenarios' Number 
mat.max <- 30 #maturities max

##Spot Term Structure Estimation
We estimate the spot term structure of interest rates using termstrc package.
As said before we do that for Moudiki/RQuantlib example commenting out the code for our data (in this part the code works!)

suppressPackageStartupMessages(library(termstrc))
##RQuantlib/Moudiki example
#maturities0 <- c(7/360,1/12,3/12,6/12,9/12,1,2,3,5,10,15)
#yields <- t(as.matrix(c(0.0382,0.0372,0.0363,0.0353,0.0348,0.0345,0.037125,0.0398,0.0443,0.05165,0.055175)))
##OUR Data
maturities0 <- c(3/12,6/12,1:3,seq(5,30,by=5))
yields <- t(as.matrix(c(-0.115064,-0.041625,0.047641,-0.118795,-0.070060,0.164114,0.844613,1.279355,1.478322,1.522483,1.527884)/100))
dates0 <- c(as.Date(Sys.Date()))
zydf <- zeroyields(maturities0,yields,dates0)
zyest <- estim_nss(zydf,method = "asv",tauconstr =  c(0.2, 7, 0.1))
times <- seq(from = 0, to = horizon, by = delta_t)
beta <- zyest$optparam
fwr <- fwr_asv(beta,times)
# that if you want to estimate using Rquantlib
#  curves <- RQuantLib::DiscountCurve(params, tsQuotes, times)
#  maturities <- curves$times
#  marketzerorates <- curves$zerorates
#  marketprices <- curves$discounts

Dimensional Analysis

Before starting the estimation and the simulation, we want to do some consideration on measure unities of parameters and functions, in particular regarding the time.
From \(B(t,T)\) and \(R(t,T)\) definition:

  • \(a\) is \([time]^{-1}\)

  • \(B(t,T)\) is \([time]\)

  • \(fwd(0,t)\) is \([time]^{-1}\)

  • \(phi(t)\) is \([time]^2\)

  • \(\sigma^2\) is \([time]^-3\)

##Short-term rates simulation
Now is time to set the parameters and simulate \(r(t)\) using package ESGtoolkit from Moudiki

#### settings of parameters  for HW model#### 
#a <- pricing$a
#sigma <- pricing$sigma
a <- .01
sigma <-a*.1
############# Hull-White short-rates simulaton ####
# Simulation of gaussian shocks with ESGtoolkit
set.seed(4)
eps <- ESGtoolkit::simshocks(n = nb.sims, horizon = horizon, 
                             frequency = "monthly")
# Simulation of the factor x with ESGtoolkit
x <- ESGtoolkit::simdiff(n = nb.sims, horizon = horizon, 
                         frequency = "monthly",  
                         model = "OU", 
                         x0 = 0, theta1 = 0, 
                         theta2 = a, 
                         theta3 = sigma,
                         eps = eps)
# I use forward rates. With the low monthly frequency
# I consider them as being instantaneous forward rates
fwdrates <- ts(replicate(nb.sims, fwr), 
               start = start(x), 
               deltat = deltat(x))
#curves$forwards instead of fwr if you use Rquantlib estimation
# alpha
t.out <- seq(from = 0, to = horizon, by = delta_t)
param.alpha <- ts(replicate(nb.sims, 0.5*(sigma^2)*(1 - exp(-a*t.out))^2/(a^2)), 
                  start = start(x), deltat = deltat(x))
alpha <- fwdrates + param.alpha
# The short-rate
r <- x+alpha
summary(as.vector(r))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.0067190  0.0006195  0.0096760  0.0086840  0.0155000  0.0278900
summary(as.vector(x))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.147e-02 -1.248e-03  0.000e+00 -1.684e-05  1.208e-03  9.659e-03
summary(as.vector(alpha))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -0.0023710  0.0006852  0.0100600  0.0087000  0.0160400  0.0184800

0$ as shown in Brigo and Mercurio.–>

ESGtoolkit::esgplotbands(r, xlab = "Times(years)", ylab = "short-rate quantiles", 
                         main = "short-rate quantiles")

plot of chunk unnamed-chunk-5

library(ggplot2)
library(magrittr)
rg<-data.frame(r)
rg<-data.frame("r"=c(rg[,1]),"Sim"=1,"Times"=times)
for (i in 2:nb.sims){
  rg<-rbind(rg,data.frame("r"=c(r[,i]),"Sim"=i,"Times"=times))
}
rg%>%
  ggplot(aes(Times,r,group=Sim))+geom_line(aes(colour=Sim))

plot of chunk unnamed-chunk-5

Now is time to compute the \(R(t,T)\)…

B <- function (t,TT,a){
  rt <- -(1.0-exp(-a*(TT-t)))/a
  return(rt)
}
phi <- function(t,a,sigma){
  phit <- sigma*sigma/(2.*a)*(1.-exp(-2.*a*t))
  return(phit)
}
system("rm HW.new.csv",ignore.stderr = T,wait = T)
TT <- times
Tmat <- seq(.5,mat.max,.5)
for (iSim in 1:nb.sims){
  for (iT in 1:length(TT)){
    rnew <- r[iT,iSim]
    #if(rnew<0) rnew=-rnew
    pp0 <- data.frame("starting.t.from.today.Years"=1,"Sim"=200,"Maturity.Years"=2,"Rate"=.0099)
    for (iT2 in 1:length(Tmat)){
      rates <- 1./(Tmat[iT2])*
        (
          B(TT[iT],TT[iT]+Tmat[iT2],a)*
            (-rnew+fwr[iT]+
             B(TT[iT],TT[iT]+Tmat[iT2],a)*phi(TT[iT],a,sigma)/2)+
          (TT[iT]+Tmat[iT2])*spr_asv(beta,(TT[iT]+Tmat[iT2]))-
          TT[iT]*spr_asv(beta,(TT[iT]+1.e-10))
      )

#Code if you are using Rquanlib estimation of forward rates      
#       rates <- 1./(Tmat[iT2])*(
#         B(TT[iT],TT[iT]+Tmat[iT2],a)*(rnew+curves$forwards[iT]+sigma^2/4*B(TT[iT],TT[iT]+Tmat[iT2],a)*B(0,(2*TT[iT]),a))-
#         (TT[iT]+Tmat[iT2])*curves$zerorates[(6*iT2)]+
#         TT[iT]*curves$zerorates[iT] 
#       )
      pp <- data.frame("starting.t.from.today.Years"=TT[iT],"Sim"=iSim,"Maturity.Years"=Tmat[iT2],"Rate"=rates)
      pp0=rbind(pp0,pp)
    }
    write.table(pp0[-1,],"HW.new.csv",row.names =F,append = T,col.names = F,dec=".",sep=",")
  }
}

Now we import saved simulations:

library(dplyr)
library(ggplot2)
HW.new <- as.tbl(data.frame(read.table("HW.new.csv",sep=",",header = F,stringsAsFactors = F,dec = ".",colClasses = "numeric")))
names(HW.new)<-c("Years.from.Today","Simulation","Maturities.Years","Rates")
HW.new
## Source: local data frame [2,892,000 x 4]
## 
##    Years.from.Today Simulation Maturities.Years        Rates
##               (dbl)      (dbl)            (dbl)        (dbl)
## 1                 0          1              0.5 0.0006907080
## 2                 0          1              1.0 0.0009920077
## 3                 0          1              1.5 0.0010296326
## 4                 0          1              2.0 0.0010603821
## 5                 0          1              2.5 0.0011802832
## 6                 0          1              3.0 0.0013978447
## 7                 0          1              3.5 0.0016919981
## 8                 0          1              4.0 0.0020368095
## 9                 0          1              4.5 0.0024097882
## 10                0          1              5.0 0.0027937470
## ..              ...        ...              ...          ...
## Error in filter(., Simulation == i): oggetto "Simulation" non trovato

Note that rates are not %

… today…

HW.new[HW.new$Years.from.Today<0.5/12 ,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-9

… 1 month from today…

HW.new[HW.new$Years.from.Today>0.5/12 & HW.new$Years.from.Today<1.5/12,] %>%  
   ggplot(aes(x=Maturities.Years,y=Rates,group=Simulation,colour=Simulation ))+ geom_line()

plot of chunk unnamed-chunk-10

… 2 month from today…

HW.new[HW.new$Years.from.Today>1.5/12 & HW.new$Years.from.Today<2.5/12 ,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation,colour=Simulation ))+geom_line()

plot of chunk unnamed-chunk-11

… 3 month from today…

HW.new[HW.new$Years.from.Today>2.5/12 & HW.new$Years.from.Today<3.5/12 ,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-12

… 4 month from today…

HW.new[HW.new$Years.from.Today>3.5/12 & HW.new$Years.from.Today<4.5/12 ,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-13

… 5 month from today…

HW.new[HW.new$Years.from.Today>4.5/12 & HW.new$Years.from.Today<5.5/12 ,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-14

… 1 year from today…

HW.new[HW.new$Years.from.Today==1,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-15

… 5 years from today…

HW.new[HW.new$Years.from.Today==5,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-16

… 10 years from today…

HW.new[HW.new$Years.from.Today==10,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-17

… 20 years from today…

HW.new[HW.new$Years.from.Today==20,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Simulation ))+geom_line(aes(colour=Simulation))

plot of chunk unnamed-chunk-18

Simulation number 10

HW.new[HW.new$Simulation==10,] %>%  
   ggplot(aes(Maturities.Years,Rates,group=Years.from.Today ))+geom_line(aes(colour=Years.from.Today))

plot of chunk unnamed-chunk-19

Standard

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s