A diffusion model: G2++

[This article was first published on T. Moudiki's Webpage - R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Version 1.0.0 of esgtoolkit is out. What’s new?

  • Potentially breaking change (hence version 1.0.0): rename ESGtoolkit to esgtoolkit (announced here). Also, more consistent with the package naming in R universe

  • include ycinter and ycextra, respectively for yield curve interpolation and extrapolation (comes from package ycinterextra, which will be discontinued (my choice) and has already been unilaterally removed by the CRAN landlords)

Here is an example of use of esgtoolkit, for the simulation of a G2++ model:

rm(list=ls()) # beware, removes everything in your environment

library(esgtoolkit)

# Observed maturities
u <- 1:30

# Yield to maturities
txZC <- c(0.01422,0.01309,0.01380,0.01549,0.01747,0.01940,
          0.02104,0.02236,0.02348, 0.02446,0.02535,0.02614,
          0.02679,0.02727,0.02760,0.02779,0.02787,0.02786,
          0.02776,0.02762,0.02745,0.02727,0.02707,0.02686,
          0.02663,0.02640,0.02618,0.02597,0.02578,0.02563)

# discount factors
p <- c(0.9859794,0.9744879,0.9602458,0.9416551,0.9196671,
       0.8957363,0.8716268,0.8482628,0.8255457,0.8034710,
       0.7819525,0.7612204,0.7416912,0.7237042,0.7072136
       ,0.6922140,0.6785227,0.6660095,0.6546902,0.6441639,
       0.6343366,0.6250234,0.6162910,0.6080358,0.6003302,
       0.5929791,0.5858711,0.5789852,0.5722068,0.5653231)

# Creating a function for the simulation of G2++
# Function of the number of scenarios, the method (
# antithetic or not)
simG2plus <- function(n, methodyc = c("fmm", "hyman", "HCSPL", "SW"))
{
    # Horizon, number of simulations, frequency
    horizon <- 20
    freq <- "semi-annual" 
    delta_t <- 1/2
    # Parameters found for the G2++
    a_opt <- 0.50000000
    b_opt <- 0.35412030
    sigma_opt <- 0.09416266
    rho_opt <- -0.99855687
    eta_opt <- 0.08439934
    # Simulation of gaussian correlated shocks
    eps <- esgtoolkit::simshocks(n = n, horizon = horizon,
                     frequency = "semi-annual",
                     family = 1, par = rho_opt)
    # Simulation of the factor x
    x <- esgtoolkit::simdiff(n = n, horizon = horizon, 
                 frequency = freq,  
                 model = "OU", 
                 x0 = 0, theta1 = 0, theta2 = a_opt, theta3 = sigma_opt,
                 eps = eps[[1]])
    # Simulation of the factor y
    y <- esgtoolkit::simdiff(n = n, horizon = horizon, 
                 frequency = freq,  
                 model = "OU", 
                 x0 = 0, theta1 = 0, theta2 = b_opt, theta3 = eta_opt,
                 eps = eps[[2]])
    # Instantaneous forward rates, with spline interpolation
    methodyc <- match.arg(methodyc)
    fwdrates <- esgtoolkit::esgfwdrates(n = n, horizon = horizon, 
    out.frequency = freq, in.maturities = u, 
    in.zerorates = txZC, method = methodyc)
    fwdrates <- window(fwdrates, end = horizon)
    # phi
    t.out <- seq(from = 0, to = horizon, 
                 by = delta_t)
    param.phi <- 0.5*(sigma_opt^2)*(1 - exp(-a_opt*t.out))^2/(a_opt^2) + 
    0.5*(eta_opt^2)*(1 - exp(-b_opt*t.out))^2/(b_opt^2) +
      (rho_opt*sigma_opt*eta_opt)*(1 - exp(-a_opt*t.out))*
      (1 - exp(-b_opt*t.out))/(a_opt*b_opt)
    param.phi <- ts(replicate(n, param.phi), 
                    start = start(x), deltat = deltat(x))
    phi <- fwdrates + param.phi
    colnames(phi) <- c(paste0("Series ", 1:n))
    # The short rates
    r <- x + y + phi
    colnames(r) <- c(paste0("Series ", 1:n))
    return(r)
}

set.seed(3)
# simulations of the short rate
ptm <- proc.time()
r_SW <- simG2plus(n = 100, methodyc = "hyman")
proc.time() - ptm

par(mfrow=c(1, 2))
esgtoolkit::esgplotbands(r_SW, xlab = 'time', ylab = 'short rate')
matplot(r_SW, type = 'l', xlab = 'time', ylab = 'short rate')

image-title-here

# Stochastic deflators :
deltat_r <- deltat(r_SW)
Dt_SW <- esgtoolkit::esgdiscountfactor(r = r_SW, X = 1)
Dt_SW <- window(Dt_SW, start = deltat_r, deltat = 2*deltat_r)

# Market prices
horizon <- 20
marketprices <- p[1:horizon]
# Monte Carlo prices
montecarloprices_SW <- rowMeans(Dt_SW)
montecarlo_ci <- apply(Dt_SW, 1, function(x) t.test(x)$conf.int)

# Confidence intervals
estim_SW <- apply(Dt_SW - marketprices, 1, function(x) t.test(x)$estimate)
conf_int_SW <- apply(Dt_SW - marketprices, 1, function(x) t.test(x)$conf.int)
conf_int_SW

par(mfrow=c(1, 2))
plot(marketprices, col = "blue", type = 'l', 
     xlab = "time", ylab = "prices", main = "Prices")
points(montecarloprices_SW, col = "red")
#lines(montecarlo_ci[1, ], color = "red", lty = 2)
#lines(montecarlo_ci[2, ], color = "red", lty = 2)
matplot(t(conf_int_SW), type = "l", xlab = "time", 
        ylab = "", main = "Confidence Interval \n for the price difference")
lines(estim_SW, col = "blue")

image-title-here

difference <- (Dt_SW - marketprices)

tests on the “martingale” difference (this implementation is very slow)

lapply(1:nrow(difference), function(i) vrtest::DL.test(difference[i, ]))

boxplot(t(difference))
abline(h = 0, color = "red", lty = 2)

image-title-here

To leave a comment for the author, please follow the link and comment on their blog: T. Moudiki's Webpage - R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)