A diffusion model: G2++

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')

# 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")

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)