Efficient, Consistent and Flexible Credit Default Simulation: TRNG & RcppParallel

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

}

Efficient, Consistent and Flexible Credit Default Simulation: TRNG & RcppParallel

Riccardo Porreca

Roland Schmid

Packages and sources

Data load

Full simulation by rating

Consistent simulation of the BBB sub-portfolio

Contribution of individual counterparties to the BBB total ES99

Full distribution for the top 3 BBB counterparties (by highest contribution)

What-if scenario: top 3 BBB counterparties downgraded (higher PD)

What-if scenario: impact on the BBB ES99

What-if scenario: new contributions

Benchmarks

Number of threads

Size of the sub-portfolio

Number of sub-simulations

Source code

simulationKernel.R

ES99.R

simulationKernel.cpp