Packages and sources
library(dplyr) # portfolio manipulation
# install.packages("dplyr")
library(microbenchmark) # benchmarking
# install.packages("microbenchmark")
if (!require(rTRNG)) {
devtools::install_github("miraisolutions/rTRNG",
dependencies = TRUE)
# install.packages("devtools")
library(rTRNG)
}
Rcpp::sourceCpp("code/simulationKernel.cpp") # rTRNG/RcppParallel C++ code
source("code/simulationKernel.R") # R wrapper
source("code/ES99.R") # utility for Expected Shortfall 99% calculation
Data load
# load pf, r, Z
# loaded <- load("data/inputDataSmall.RData")
loaded <- load("data/__inputDataBig.RData")
loaded
## [1] "pf" "Z" "r"
stopifnot(all(c("pf", "r", "Z") %in% loaded))
J <- nrow(pf)
M <- nrow(Z)
stopifnot(ncol(Z) == J)
stopifnot(all(dim(Z) == dim(r)))
colnames(pf)
## [1] "j" "V0" "R" "beta" "PD" "rtng"
dim(Z) %>% setNames(c("M", "J"))
## M J
## 10000 6000
# number of idiosyncratic returns to be simulated for each market scenario
K <- 100
# random seed
s <- 12358
# timer for the net simulation time in simulationKernel
timer <- function(...) {
system.time(...)["elapsed"]
}
Full simulation by rating
L_rtng <- simulationKernel(pf, Z, r, J, K, agg = pf$rtng, seed = s,
timer = timer)
## elapsed
## 298.42
ES99_rtng <- ES99(L_rtng)
sort(ES99_rtng, decreasing = TRUE)
## BBB AA AAA A BB B
## 9878920788 8620051762 7468245838 4796596354 1190520347 641062514
## D CCC CC C
## 158314661 135836081 7228975 5488510
Consistent simulation of the BBB sub-portfolio
L_BBB <- simulationKernel(pf %>% filter(rtng == "BBB"), Z, r, J, K, seed = s,
timer = timer)
## elapsed
## 144.39
all.equal(c(L_BBB), L_rtng[, "BBB"], check.attributes = FALSE)
## [1] TRUE
ES99_BBB <- ES99(L_BBB)
ES99_BBB
## PF
## 9878920788
Contribution of individual counterparties to the BBB total ES99
pfBBB <- pf %>% filter(rtng == "BBB")
L_jBBBtail <- simulationKernel(pfBBB, Z, r, J, K,
agg = pfBBB$j, mk = tail99(L_BBB), seed = s,
timer = timer)
## elapsed
## 5.63
ContrES99_jBBB <- colMeans(L_jBBBtail)
head(ContrES99_jBBB)
## 2 8 20 32 33 38
## 4277383374 999636420 212202557 21123706 135392425 50909632
all.equal(sum(ContrES99_jBBB), ES99_BBB, check.attributes = FALSE)
## [1] TRUE
Full distribution for the top 3 BBB counterparties (by highest contribution)
top3jBBB <- names(sort(ContrES99_jBBB, decreasing = TRUE))[1:3]
pftop3BBB <- pfBBB %>% filter(j %in% top3jBBB)
L_top3BBB <- simulationKernel(pftop3BBB, Z, r, J, K,
agg = pftop3BBB$j, seed = s,
timer = timer)
## elapsed
## 0.17
ES99_top3BBB <- ES99(L_top3BBB)
pftop3BBB %>%
select(j, V0, R) %>%
mutate(ES99 = ES99_top3BBB, ContrES99 = ContrES99_jBBB[top3jBBB],
Div = ContrES99/ES99, "Contr/V0" = ContrES99/V0)
## j V0 R ES99 ContrES99 Div Contr/V0
## 1 2 9444041000 278250000 4720889738 4277383374 0.9060545 0.4529188
## 2 8 3260697000 91000000 1672364832 999636420 0.5977382 0.3065714
## 3 70 298111300 13483307 963482756 436598485 0.4531461 1.4645486
What-if scenario: top 3 BBB counterparties downgraded (higher PD)
pfBBBwi <- pf %>% filter(rtng == "BBB") %>%
mutate(PD = replace(PD, j %in% top3jBBB, 0.01))
pftop3BBBwi <- pfBBBwi %>% filter(j %in% top3jBBB)
L_top3BBBwi <- simulationKernel(pftop3BBBwi, Z, r, J, K,
agg = pftop3BBBwi$j, seed = s,
timer = timer)
## elapsed
## 0.16
What-if scenario: impact on the BBB ES99
L_BBBwi <- L_BBB + (rowSums(L_top3BBBwi) - rowSums(L_top3BBB))
ES99_BBBwi <- ES99(L_BBBwi)
cbind(ES99_BBBwi, ES99_BBB)
## ES99_BBBwi ES99_BBB
## PF 11943875892 9878920788
What-if scenario: new contributions
L_jBBBtailwi <- simulationKernel(pfBBBwi, Z, r, J, K,
agg = pfBBB$j, mk = tail99(L_BBBwi), seed = s,
timer = timer)
## elapsed
## 7.41
ContrES99_jBBBwi <- colMeans(L_jBBBtailwi)
Benchmarks
# reduce to a reasonable size for benchmarking
M <- pmin(1000, M)
J <- pmin(600, J)
K <- pmin(10, K)
Z <- Z[1:M, 1:J]
r <- r[1:M, 1:J]
pf <- pf[1:J, ]
# number of repetitions
nBench <- 10
Number of threads
withThreadNr <- function(numThreads) {
RcppParallel::setThreadOptions(numThreads = numThreads)
invisible(simulationKernel(pf, Z, r, J, K, seed = s))
}
mb <- microbenchmark(withThreadNr(1),
withThreadNr(2),
withThreadNr(4),
times = nBench)
boxplot(mb, unit = "s")
Size of the sub-portfolio
set.seed(85321)
subJ <- ceiling(J/2^(0:4))
# sub-portfolios with constant PD (in order not to bias the benchmark)
subpfs <-
lapply(as.list(subJ),
function(jlen, pf) {
pf %>% sample_n(jlen)
},
pf %>% mutate(PD = mean(PD))) %>% setNames(subJ)
RcppParallel::setThreadOptions(numThreads = 1)
withSubJ <- function(subJ) {
invisible(simulationKernel(subpfs[[as.character(ceiling(subJ))]],
Z, r, J, K, seed = s))
}
mb <- microbenchmark(withSubJ(J/2^0),
withSubJ(J/2^1),
withSubJ(J/2^2),
withSubJ(J/2^3),
withSubJ(J/2^4),
times = nBench)
boxplot(mb, unit = "s")
Number of sub-simulations
set.seed(85321)
mklen <- ceiling(M*K/2^(0:4))
mks <-
lapply(mklen,function(s) {
sort(sample.int(M*K, s))
}) %>% setNames(mklen)
RcppParallel::setThreadOptions(numThreads = 1)
withSubMK <- function(mkl) {
invisible(simulationKernel(pf, Z, r, J, K, seed = s,
mk = mks[[as.character(ceiling(mkl))]]))
}
mb <- microbenchmark(withSubMK(M*K/2^0),
withSubMK(M*K/2^1),
withSubMK(M*K/2^2),
withSubMK(M*K/2^3),
withSubMK(M*K/2^4),
times = nBench)
boxplot(mb, unit = "s")
Source code
simulationKernel.R
simulationKernel <- function(pf, Z, r,
J,
K, mk = seq_len(K * nrow(Z)),
agg = factor(rep("PF", nrow(pf))),
seed,
timer = function(expr, ...) expr) {
# check arguments
stopifnot(all(c("j", "V0", "R", "PD") %in% colnames(pf)))
stopifnot(all.equal(dim(Z), dim(r)))
stopifnot(all(pf$j > 0 & pf$j <= ncol(Z)))
stopifnot(all(mk > 0 & mk <= nrow(Z) * K))
stopifnot(!is.unsorted(mk))
stopifnot(is.function(timer))
# aggregation levels
agg <- as.factor(agg)
agg_idx <- as.numeric(agg)
agg_set <- as.character(levels(agg))
# pre-allocate the aggregated output
L_agg <- matrix(0.0, length(mk), max(agg_idx),
dimnames = list(mk = as.character(mk),
agg = agg_set))
# timer wrapping the core simulation only
t <- timer(
simulationKernel_C(pf = pf, Z = Z, r = r,
J = J,
K = K, mk = mk,
agg = agg_idx,
seed = seed,
L_agg = L_agg)
)
if (!is.null(t)) {
print(t)
}
L_agg
}
ES99.R
ES99 <- function(L) {
es99 <- apply(L, 2, function(l) {
mean(l[tail99(l)])
})
es99
}
tail99 <- function(l, sorted = TRUE) {
idx <- tail(order(l), ceiling(0.01*length(l)))
if (sorted) {
idx <- sort(idx)
}
idx
}
simulationKernel.cpp
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::depends(rTRNG, RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;
#include <trng/yarn2.hpp>
#include <trng/normal_dist.hpp>
using namespace trng;
typedef yarn2 rngKind;
struct KernelWorker : public Worker {
// Input data
const int J, K;
const RVector<int> j, mk, agg;
const RVector<double> V0, R, beta, PD;
const RMatrix<double> Z, r;
const std::vector<rngKind> rngSplit; // pre-split sub-sequences
// Output object to write to
RMatrix<double> L_agg;
// Initialize from input and output objects (automatic wrap of Rcpp objects
// as RMatrix/RVector)
KernelWorker(const NumericVector V0, const NumericVector R,
const NumericVector beta, const NumericVector PD,
const IntegerVector j, const int J,
const NumericMatrix Z, const NumericMatrix r,
const int K, const IntegerVector mk,
const IntegerVector agg,
const std::vector<rngKind> rngSplit,
NumericMatrix L_agg)
: J(J), K(K), j(j), mk(mk), agg(agg), V0(V0), R(R), beta(beta), PD(PD),
Z(Z), r(r), rngSplit(rngSplit), L_agg(L_agg)
{}
// operator processing an exclusive range of indices
void operator()(std::size_t begin, std::size_t end) {
rngKind rngj;
normal_dist<> normal(0.0, 1.0); // TRNG normal distribution
unsigned int mk_begin = (unsigned int)begin;
unsigned int mk_end = (unsigned int)end;
int thismk, prevmk;
unsigned int m;
int thisj;
double sigmaj, thetaj, Yj, Lj;
// loop over the set of relevant counterparties j
for (unsigned int j_idx = 0; j_idx < j.length(); j_idx++) {
thisj = j[j_idx] - 1; // current counterparty (0-based indexing)
sigmaj = sqrt(1 - beta[j_idx]*beta[j_idx]);
thetaj = normal.icdf(PD[j_idx]);
rngj = rngSplit[j_idx]; // j-th subsequence
prevmk = -1;
// loop over the relevant combined simulation indices
for (unsigned int mk_idx = mk_begin; mk_idx < mk_end; mk_idx++) {
thismk = mk[mk_idx] - 1; // current combined simulation
rngj.jump(thismk - prevmk - 1); // jump to thismk from prevmk
m = thismk/K; // market scenario simulation index
// credit environment combined return
Yj = beta[j_idx]*Z(m, thisj) + sigmaj*normal(rngj);
// default indicator and corresponding loss
Lj = V0[j_idx] - ( Yj < thetaj ?
R[j_idx] :
V0[j_idx] * (1 + r(m, thisj)) );
L_agg(mk_idx, agg[j_idx] - 1) += Lj;
prevmk = thismk;
}
}
}
};
// [[Rcpp::export]]
void simulationKernel_C(const DataFrame pf,
const NumericMatrix Z, const NumericMatrix r,
const int J,
const int K, const IntegerVector mk,
const IntegerVector agg,
const unsigned long seed,
NumericMatrix L_agg) {
// counterparty indices in the (sub)portfolio
IntegerVector j = pf["j"];
// RNG for the full sequence
rngKind rngFull(seed);
// RNGs for the LeapFrog subsequences of individual counterparties
std::vector<rngKind> rngSplit(j.length());
for (unsigned int j_idx = 0; j_idx < j.length(); j_idx++) {
int thisj = j[j_idx] - 1; // 0-based C++ indexing
// split the full sequence based on the full set of counterparties (J)
rngSplit[j_idx] = rngFull;
rngSplit[j_idx].split(J, thisj);
}
// worker for the parallel simulation
KernelWorker w(pf["V0"], pf["R"], pf["beta"], pf["PD"],
j, J,
Z, r,
K, mk,
agg,
rngSplit,
L_agg);
// parallel execution over the set of relevant simulation indices mk
int grainSize = 100;
parallelFor(0, mk.length(), w, grainSize);
}
