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

Potentially breaking change (hence version
1.0.0
): renameESGtoolkit
toesgtoolkit
(announced here). Also, more consistent with the package naming in R universe 
include
ycinter
andycextra
, respectively for yield curve interpolation and extrapolation (comes from packageycinterextra
, 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 < "semiannual"
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 = "semiannual",
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)
Comments powered by Talkyard.