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")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAgVBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZmYAZpAAZrY6AAA6OgA6ZmY6kJA6kLY6kNtmAABmOgBmOpBmZjpmZmZmkJBmtv+QOgCQZgCQkDqQkGaQkNuQ2/+2ZgC2Zma225C2/7a2/9u2///bkDrbtmbb/7bb////tmb/25D//7b//9v///9hOfp/AAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3diXbbypqYUdo3zpUyyOkkcgazu3N425TN93/AECNBkFITPwsoFLX3WmfZlEVIp1j6hJmbAwAhm9zfAECpBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYLWHdANQDLpE5V8iQnlHm3gsSRvVOoFpjTDLwzg0xJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABzSv3O2L1cg8ElEhAs8qdzZPcIwElEtDiPeL/E5RBQIv3iP9PUAYBvWNZnEk3tFAIAY0vipFkYwuFEND4ohhJNrZQCAGNL4qRZGMLhRDQ+KIYSTa2UAgBjS+KkWRjC4UQ0PiiGEk2tlAIAQUIKi6gv7936ztf/5r+bAEF0iksoPvzbcZvvyY+X0CBdMoK6Ha81+3Lz2kLEFAgnaICuttsnqo/356rP//8mF5QAQXSKSmgx26+Nn87tvOl/ci0rXgBBdIpKaC7Zv2z8vbcrHqemnobAQXSKSig3Wpn+/cmpttTVG8hoEA6BQX09/fB2ua23XbfTduGF1AgnXID2pwGup92OqiAAumUFdDhJryAApkVFNB+s/0wOPxuEx7IpqSA7jb9Kuh205/H5CASkElJAa0ug6832Pu/VFd2Oo0JyKSkgFbrm73X9lKkSSugAgokVFRABwWt1jurgE68nYiAAumUFdDudkztrtDt1HuJCCiQUGkBvZOAAukIKECQgAIECShAUOEBfXv+6FJO7xwJzElAAYIeOqCXBBRIp/CATiWgQDoCChAkoABBAgoQVF5At/0B9Wk3YqoJKJBOYQHdjU5KmtpQAQXSKSqgw/uBtiadxCSgQEolBbS+gfKgmNWN6SfeEFRAgXRKCuj+IpdVUr2lB5BJSQHdXm6wH1dCvakckElBAT2ubr5cfNDbGgPZFBTQ49rm5eb63rXwQC4CChBUUEBtwgPrUlBAHUQC1qWkgF4/jelyrfQDAgqkU1JAr55IP+1SJAEF0ikpoE0xz335OWkJAgqkU1RAh7dicjMRILfCAnpwOztgNcoL6F0EFEhHQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCCo2IDuN0dPU58loEA6hQV0d6zm178Ohz8/No3Xac8XUCCdogLaZvPLz8N2swkVVECBdIoKaJfNb/84/vfrUK+QVuujtxNQIJ2SArpvsvn2fGzoS/+hSaugAgqkU1JAt81qZ1XQfr1zO+1AkoAC6RQU0D8/2vXO41/6au7aqN5IQIF0Cgro7+/t5npf0kO1DT9pJ6iAAukIKEBQQQG1CQ+sS0EB7Q8i7U8nLw1begsBBdIpKaDtaUz1RZxOYwKyKymg/fWb/Yn0eyfSA/mUFNDqMFJzKWefUpdyAvkUFdDTzUTalLqZCJBRYQE92bmdHZBZsQGNEVAgHQEFCBJQgCABBQgqPKDDO9td2lyx4DcHPDgBBQh66IBeElAgncIDOpWAAukIKECQgAIECShAUHkB3fYH1CdfCS+gQEqFBXQ3OilpakMFFEinqIC+PV+c1jntfsoCCiRUUkDr2ygPilnfFHTSe8oJKJBQSQHdX+SySqr3RAIyKSmg28sN9uNKqHflBDIpKKD9+8IPeV94IJuCAnpc27zcXN+7Fh7IRUABggoKqE14YF0KCqiDSMC6lBTQ66cxXa6VfkBAgXRKCujVE+mnXYokoEA6JQW0Kea5Lz8nLUFAgXSKCujwVkxuJgLkVlhAD25nB6xGeQG9i4AC6QgoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIEFRTQ3983337duQwBBdIpK6CbLz/vW4aAAukUFtDN5uWuZQgokE5ZAf32vzabr3/dsQwBBdIpLKC/dseV0Kf4MgQUSKe0gB7enu/ZjhdQIJ3iAno4bDfxtVABBdIpMKDNSuhm8xpYhoAC6ZQY0MNhv2lMPqtJQIF0ygxovxa6mXhUXkCBdEoN6KHdFyqgQDYFB/Tozw8BBbIpO6CTCSiQjoACBAkoQFBBAU1BQIF0BBQgSEABggQUIEhAAYIKD+jb80cn0m+uWPCbAx6cgAIEPXRALwkokE7hAZ1KQIF0BBQgSEABggQUIKi8gG77A+qBN5YTUCCdwgK6G52UNLWhAgqkU1RAT2+E1Jt0EpOAAinNG9Df36+dyx5rX/MGHsNn1UufdotQAQXSKSmg+4tcVkmd9PbwAgqkU1JAt5dPOX6BSbtBBRRIZ/aAvreCuJ8c0OPq5svFB3fTtuEFFEinoIBeXdjExQgokI6AAgQVdBqTTXhgXQoKqINIwLosG9DqTPj4e7tfP43pcq30AwIKpLNQQH9/r9LXXEn05WdwyVdPpJ+2J1VAgXSWCei+Dl0dwMgJoJ0rp5VOrLGAAuksEtBqzfNYuub9N3YTt7rPbEf9dDMRIJ9FArprStemcxu5D92J29kBK7FEQI+b7tUe0OMf9Qb39DNA0xFQIJ0lAtqeAX/8oy6ngAKPYcGA7tuNbgEFHsOCAd22t56bePFQUgIKpLPQPtCX0y7QqRcPJSWgQDqLHIXfVuncNdcRVac0TboHclICCqSz2HmgldfmJKR8W/ACCiS0zJVIzZtpVuHc3XEhUgICCqSz0LXw+01/CD7f9vtBQIGUSrqdXQICCqQjoABBBb2lRwoCCqQjoABBAgoQNHtA3yegQNkEFCBIQAGCnMYEECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKELRUQP/8aO/BlPFN4Q8CCqS0SEBHd7XL2FABBdJZIKBXbgqaLaECCqQzf0C358lst+WfUn/Z2wgokM7cAX17PtZy/MZyu+PHvvxM/YVvIaBAOjMHdP9OKauuvvd+nXMSUCAdb2sMEDRzQP/pg6NF/yygQNGcSA8QJKAAQYsFtDob9MvPw86VSMCjWCigzcn0VUDzHH3vCCiQzjIBbS9GOgZ0e+W00AUJKJDOIgGtrj769uvtuTol9FjQHOcvtQQUSGeRgO43m5fq5Pn6nPpd/SATAQXSWSSg2/pS+Dagx9XRfAeSBBRIZ4mAHpNZrXO2AT2ugubbhhdQIJ0lAtpe0NkFNM9FnA0BBdIRUB7PxQ1oM8k9Dswuwyb81j5QZpW7m73cA8HsFjqIVK1ztgE9ro9mupvyQUC5kYnCTZY6jempC2h1Tmi+M+n9XHATE4WbLHMlUn32fB3Q6mb03lSOtTNRuMmil3I2Ml6I5OeC25go3GShm4n0bwt/Zz9/f2+fXq3JBt5Xyc8FNzFRuMlit7NrE3rf6md3If0+WmM/Fyu3zNHxguR+QfhYUTdUrm7lVK117k/za+JKqAm5btk6tV65XxI+VFJAq7fyrM6Aqvao1que1VrttHVQ83HdvD4jBmTlSgrotj0BarfpzyTdT7y1k/m4bl6fEQOycgsFdH+2VRLbEXpc33xq/zwtYOJlTebjunl9RgzIyi12Q+X7A9q9yfzxz1M0J97ayXxcN6/PiAFZuUUCutskDujpWtCJdyYxH9fN6zNiQFZuoZuJpLj8vb0nSb8pXxPQh+L1GTEgK7fQ7eySXL25bcO5PduEtw/0gXh9RgzIyi0U0CTvgrRvz/p8e+5vR9Ktld7KfFw3r8+IAVm5Be8Herd+V8DpyNF24pn05uO6zXtOepFyvyR8aKGDSGnuAFqdSV9vsrfdrE6OmpZm83HdctdqhXK/JHxokYB2x8/vVhf0zMQym4/rlqNQK5f7JeFDy5xIfwxfoveC357PrqldNh/XLU+jVi33S8KHFroS6XzV8c472nWLCRzaNx/Xbfk+rV7ul4QPFXQlUgrm47rladSq5X5J+FBBVyKlYD6um9dnxICsXEFXIqVgPq6b12fEgKzcQifSr6Sf5uPKeX1GDMjKlXQl0hVvzx/tD7BLqTRenxEDsnIlXYl0hYA+Fq/PiAFZuYUOIs31VvAfB/SS+bhuXp8RA7JyS53GlOZKpLuZj+vm9RkxICu3zIn0qW5odzfzcd2Sn0VZvtwvCR9a6CDSOeeBcl2WRK1b7peEDwkoQFB5AT3dTiRwcqmAAumU9L7wh8uLQqc2VECBdIoK6OXtQKeuzAookE5JAa1v6jQoZr1rYNrRfQEF0ikpoPuLXFZJnXSGqYAC6cwb0God8TXZQaTt5ROn3qdEQIF0Cgro1UvqvS88czBRuElBAb361nR718IzAxOFmxS0D1RAWYyJwk0KCqhNeBZjonCTggLqIBKLMVG4yTKXcn75OXi4jd6Z6fppTJPu1ezngpuYKNwkS0CDpzFdPZF+2rL8XHATE4WbLB/QY/aiNxO5OJy/2Zyl+YZvzs/FZ3AxTTLJPQ7MbuaAjm/+0YjfXHk7WpKbiXBp/jLeKPdAMLuZA3plnXEzcbflmNvZASsx9yb8PnU/7yOgQDoZDiLlJKBAOgIKEFTSifQJCCiQzsw3E/mnD463/3OGt5YTUCCd2e/G9N79jifeBSQRAQXSmXkTfv/Oqe7VuxvlOBgvoEA6c+8Drd8HbrwWWp1en+e4koAC6cx/EKk58b2/+Ki+oD10EnwKAgqks8BR+CtXI8Wv5byTgALpLHIa0yih2fIpoEBKS50H2m65Z63nQUCBlJxIDxAkoABBAgoQJKAAQQIKELTM7exWcyKogALpZAporkuRBBRIJ1tA8xRUQIF0ltkH+vbc3TykiunL8APLElAgnUUCeszl6d51+zqdx5DmWAUVUCCdRQK6PYtl82iX5UCSgALpLLQPdHhH0OZe9O5ID5ROQAGClgjonx+jTfhq411AgdIttA/07CBSndOtfaBA4ZY6Ct+tb1anMVVH4XfeVA4o3TLnge7OTqF/rYvqPFCgcAvdTKR+c85GFc7jZnyWFVABBRJa7G5M7fWczZb8PtO7GgsokJDb2QEECShAUJ5N+GwEFEhnoYAO7miXNaECCqSzTEDP7gia6fhRTUCBdBYJ6J8f9dmflX3G9/M4CCiQ0iIB3Q+iWcX09fJTFiKgQDoLXQs/2PGZ6VbKDQEF0lnobkzDy47y3Eq5IaBAOtnuB5qHgALpCChAkE14gCAHkQCCnMYEEOREeoCgDJdy5rwYXkCBdNxMBCDI7ewAgtxQGSBIQAGCBBQgaN6Anh19X8NxeAEF0hFQgCABBQiyDxQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCCo0IA2t7r/9mvq8wQUSKeggB6j2RVz2705yNPEZQgokE6RAd2e3l5p4kqogALplBjQ/TGcL+1HJq6DCiiQToEBPf7Zv63nbvPl55RlCCiQToEB3bfrn5U/P6atggookE6BAT1b69xN2wsqoEA6BQZ0Owzo/rQ5fwsBBdIpMqCDZgookE2BAd1ZAwVWocCA7jeb1/6DW/tAgVwKDOjwyPvwiPwtBBRIp6yAbupT6PtqHlM6bQteQIGEygtorc5mfUXn67/7vCEBBdIpKKCHQUPrgO4m91NAgYTKCmjj7bk7HD/1ZkwCCiRUYkDvIKBAOgIKECSgAEECChBUeEDfnj86EXRzxYLfHPDgBBQg6KEDeklAgXQKD+hUAgqkI6AAQQIKECSgAEHlBXTbH1CffCW8gAIpFRbQ3eikpKkNFVAgnaIC+vZ8cVrntPspCyiQUEkBrW5APyxmfXfQSW+JJKBAQiUFdH+Ryyqpk26pLKBAOiUFdHu5wX5cCZ20G1RAgXQKCuhxdfPyHTh33tYYyKWggB7XNi831/euhQdyEVCAoIICahMeWJeCAuogErAuJQX0+mlMl2ulHxBQIJ2SAnr1RPpplyIJKJBOSQFtinnuy89JSxBQIJ2iAjq8FZObiQC5FRbQg9vZAatRXkDvIqBAOgIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECCtxgsxK5x+GcgAL/vtzd7OUeiHMCCiziEX/6BBRYxCP+9AkosIhH/OkrN6C/v29eJj/pEV9CKMMj/vQJKLCIR/zpE1BgEY/401dQQI/FvObrX1O+uQd8CaEMj/jTJ6DAIh7xp6+ggB72AgqsSUkBrddBv/06PbAPFMipqIAeDtvN5svP5q8CCmRWWEAPb8+btpsCCmRWWkAPf360m/ECCmRWXEAPh91ms3kVUCC7AgNab8Y/CSiQW4kBrTfjv/6rgEJJHvGnr8iANpvxGwGFgjziT1+hAa034wUUCvKIP32lBrQ+JVRAoRyP+NNXbkBDHvElhDI84k+fgAKLeMSfPgGFR5fizTAfSsqhTbesdompF5iSgPL55M7VCiUc22SL6paYeoEfeXv+6HZ2s44cFMKsHxHQjoDCv8esHxHQzscBvWQq8fmY9SMCGmUq8fmY9SMCGmUq8fmY9SMCGmUq8fmY9SMCGmUq8fmY9SOfOqDb/oD60/Qnm0p8Pmb9yOcN6G50UtLUhppKfD5m/chnDWh9D7tzk05iMpX4jMz6kU8a0OpG9MNiVm8Tf3qf+JuYSnw+Zv3IJw3o/iKXVVJfpyzCVOLzMetHPmlAt5cb7MeV0Em7QU0lPp/wLTceV8KxTbaobompF9g5rm5e3oJ+N20bXkD5fHLXaoUSjm2yRXVLTL3AznFt83Jzfe9aePhY7lqtUMKxTbaobompF9gRUIjIXasVSji2yRbVLTH1Ajs24SEid61WKOHYJltUt8TUC+w5iAQRuXO1OimHNt2y2iWmXmDv+mlMk97bWEAhl0f86SspoFdPpJ92KdIjvoSwhNxrjZ3c43CupIA2xTz35eekJaxt+KEQGVJ5Xe6BOFdUQIe3Ymq4mQiQT2EBPbidHbAa5QX0LgIKpCOgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKEPTpAgqQTvJGpV5gSrkHG3gsyRuVeoF8zG6JGRjUGRjUWxikhZmWMzCoMzCotzBICzMtZ2BQZ2BQb2GQFmZazsCgzsCg3sIgLcy0nIFBnYFBvYVBWphpOQODOgODeguDtDDTcgYGdQYG9RYGaWGm5QwM6gwM6i0M0sJMyxkY1BkY1FsYpIWZljMwqDMwqLcwSAszLWdgUGdgUG9hkBZmWs7AoM7AoN7CIC3MtJyBQZ2BQb2FQVqYaTkDgzoDg3oLg7Qw03IGBnUGBvUWBmlhpuUMDOoMDOotDNLCTMsZGNQZGNRbGKSFmZYzMKgzMKi3MEgLMy1nYFBnYFBvYZAAggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIE9Ko/PzZf/6r/9vv/to+//Bx+wm4z9nL5SVG/v2++/Wq+SPtd1PbVF7lit3kaPrd7yvHbeU3x3SRT5qC+PVffSPc9GNQz98zUamRfu/+HdQ3qBAJ6VT8tt81MyDYthzPunWn59jyYu6efp9E/rECJg3p8VqsdS4M6dMdMrYf29do/FEVAr+qmZTcT8k3Lwe/m69Py/Nf3drgqsGsWsxYFDuqpn/24GtSBO2ZqNVX7xysb1AkE9EPvTsvWcQKdejXPtDwt//q0HM6943cwfMJxOVe3pDIraFB3m/aft/0qlkE9X2psptafcwroSgf1BgL6oRVMy9Om0dVp2e9IOnS76wZbQ7vNGjeNyhnUwVbmbrAKalAHSw3N1G7dvv/AOgf1BgL6oezT8j8/n2bZ1Wm57X+t11Py7z+GM3Gdv9jLGdT9KQrVyv1rtyCDelpqYKY238HX/z0I6DoH9QafPKDbbiIdX8D2b8dflE/tnqV9u9PotZ1x9eOznTVXpuXZJ1WTpvpAs+z95nxXUb1huBnM5bqBT8Np+TI4vrnvNydPCz1+t93Eq9YCXv+cBTTPL/bHGdTdIATbwcqoQb1rpjZPeN0Nv0Kpq6CfPKD9a7jv58uunYaX03LbP+xdTsvzTzrOoH9sTjPo7GjEYPd+O4vbD3z7t8G0rNZ82rWgwbTsF7o7Ter6HJFRQM9n7UIeaVBP+oIZ1PsHtR7Cs4BmGdQEPnlA61/ilf48jGZCXpuW/2kzmlaHK9Ny9Enbzd/qyfg0mJXdvw0Pj76cfeA/DqZl/bz+R+dltNBuDeAw+B5Gp4osf3jz4Qa1/aYG+0oM6l2Devz6xwdnAc0yqAl88oB2vakPX9cv4PGFHKzIDfYsbZqXuz8uWxtPy/EnbbvF1v9Y/wretz8A1UbQS7uM+nO6P09bVs0WZL9ptL9Y6OXv7VFAB+tNy3m4QT30T2sY1DsHdduuk57vJFh+UBP45AHtNiyOU+J/9LPm9fq07A/CDk4ZvpiWo0/a9hNo3/+t3Yl1OkBxnFvV03bdZ1S/yQfTslrsy+FwNi27X9bnc/BwGdCLT1jCow1q8/TTCpJBvW9Q238+H8Ysg3q/zx7Q9hfjfvP1X5qzLJpfhNem5Uv/jNOP0sW0HH3Stv/I4Bfs+Bhls5DBkdHd2bSsFtb9xHTTsnv+xa/tcUDfu6huVo82qH053vlai3icQR0ceRokM8ug3u+zB7TZDqr2dv9bt0VUvbgfXN/x4bQcf1J/6Ha4i6ffndU9q14dGCyqe3a32bjrt5i6afnaP3d07HL8kfOvtZBHG9Tqk0aH5gzq2bMnDer2tGt1ENAsg3q/zx7Q5tU8vrwvzXlq7e/5dNPydPLJUH+a4Wln/ZWjFLvztYrTtBzc2mKNAX20QR3106DeNah9NwX0EdR7kuqLJOqdTO2epuTTcnBks5uWww8dFzJY7nhatptGwYBmOLj5WIM62JN35btdzIMM6imUFwEt8TD8pw9ofQnEvp0Xr93lEnNMy9H06H7Tf/2rWcgHv9fbTaNy1kAfaVAHpzd2DOodgzq+uUm349MaaKGqiVhPxmonU7ujaZYNo1HqThdu9NPyvT1L7aZROQF9oEGt1r7GBzcM6h2DKqAPZrf5+v9+1C9jc91EvVWRfFr++TE6S2Pwyc0utsFnjI5ttgv88n8upuVKj8I/0KB2x5XPGNQ7BvW9gDoKX6jjC/4/n9uTL77893aWJZ+Wg1/jzUJPn1xtIZ2dXdedp3x28vbuNNsG03Kd54E+zqCeLjwfMqgJZqrzQB9Ed3JGu6/8qftgNy27PU13Tst6R1I9Q7bNhR7drY/r+VZ9oJue9bdxMS2bHVHjabnKK5EODzOoF2tjDYOaYKa6EulRDC8u7l7Sblo2xx9fE0zLwZkgm/aM6IHu9OPGf7nYNd/963haticDnowCmu0K48cY1LPF9T/uBjXFTB0F1LXwxdp3L3f/q/Y0LeuZerYZE56Wg5NB+q/W+K/dek47L8/ucXP6PrdXpuVhfOOgUUCz7Vh6iEFtvtGLgBrUFDN1FNBCd4EK6HA/1+kq4L5D9bx8SjItu71D/bHGZp6+DC42br/a5ckh7de6nJbjLaNRQLfnG/TLeYhBHZ1T3v24G9QUM3UU0GyDeicBLdz2oy2f36Xe5zszgzqDxxxUAS3c29MnKxwAAAT7SURBVPMHBy9Lvc13bgZ1Bo85qAJaug/eEbbYX+vZGdQZPOSgCmjpjnPvvV/s5b7bdm4GdQYPOagCWry353e2ft6eizyzbhUM6gwecVAFtHy76xcR//lR5KUdK2FQZ/CAgyqgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSg5LLbnHua98v9+bH5+te8X4JPR0DJRUApnoCSi4BSPAEll2NAX3N/D3AXASUXAaV4AkouAkrxBJRcrgT07bndT/nnx+bLz+rh0+H398H+0e3m26/98fHxH4/29b7T00L2o32pZ48H+0DPntd8vN4h++1X+v9LHpqAksu1NdDjx15Of1YBPf636YtZBfQf3cP2XzZdF/vHzSLGj/uA1kU+Pa/6+L90n2qVmEkElFyuBbRZ86waV60NHhP49x9t2pqCbjd/e27XKk99bErYd7Ht4PhxF9DBx+vHf36cPq/rNNxGQMlldBpTs/28r+O4bUpWN7Kq3L7bEt/2n1h1r/6k7t92gyBWnzJ+3AV0266S7trnNQF9bT/ykmEgKJeAksvVgNbt3Le5rALafPj4t7qW2/7z9v3fjquU1b9tuxXatpTjx+0f3ZL6v1UBbXcC7GY/GZUHI6Dkcj2gxxz+vdtbeYpdlcOX0x+HfiW1sq8/uBsdBBo/bgM6iGSzo7UKaLvM49dzHIkpBJRc3jmNaXc6mDMIWrtS2q9WtntJu0976g6tnxY5fjxeMe2e1+52HX09uIWAkss7Aa3WCNt1xKZww7/2653DQ0Tt6uv27KD7xeMmoINctr0UUOIElFzeO5F++05Am0iedmBebP+f9gm8dF9g8LgPaH9FfLMWK6DECSi5vBPQ/dkm/HtroNdb16+Xvlw+tgZKegJKLtcDeuzZf/hvp2PvF/tAT5vw79xbqY7moIPd48t9oM0yBZQ4ASWX6wGtPrpvVyH7Kzv77vUBrfaUnj17sG1ex3X8+MpR+G23ZiqgBAkouVwNaN2wrmnVjs6n7sPtWfBd7Hab4aXtL8MznJpgjh9/cB6ogBIkoORyNaDb7jL3Kpz1kaIqaYMrkbrY1Vvmr83H6g/2n1P9y9Pl42tXIrVXKAkoQQJKLuM70lc53J8u2Xytg/a38zuGnJ8+f37UfXt63Jdy8PiDa+EFlCABJZcrAR3eMemYsmpFdD/s5zCgF3dfOt0VpPmc0eMP7sYkoAQJKLlcCeiu329Z/63ekq872W3rDwPaLeF0+Xqbxterjz+4H6iAEiSgrNfgPFBYIwFlvQSUlRNQ1ktAWTkBZb0ElJUTUNZLQFk5AWW9BJSVE1CAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAoP8Pe9Iv68MruK8AAAAASUVORK5CYII=)
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")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAdVBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZmYAZpAAZrY6AAA6OgA6kNtmAABmOgBmOpBmZjpmZmZmkJBmtrZmtv+QOgCQkDqQkNuQ29uQ2/+2ZgC2Zma225C2/9u2///bkDrbtmbb/7bb////tmb/25D//7b//9v///9VZdsWAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3diXbbSHigUbozndiRk1nkyZhJmknTtt7/EYdYCZKQLPwsolCle8/p46YoQcSij9i5ewEgZJf7BQCUSkABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIOjRAd3vBp8f/JsAVvbQgB52lzQUqMkDA/rjy+7aH38tG8TNAADCklfucQH99W13UcyfX0+P//x7ySByT22gLqkz97iAHm9y2ST1eckgHvCGAXxYJQV0f7vBfloJXbQbVECBdAoK6Gl18+nmi4dl2/ACCqRTUEBPa5u3m+vHZYeRBBRIR0ABggoKqE14YFsKCqiDSMC2lBTQ+dOYbtdK3yCgQDolBXT2RPpllyIJKJBOSQHtinnp0/dFQxBQIJ2iAjq9FVNn6c1EBBRIp7CAvtx5OzsBBdIpL6B3EVAgHQEFCBJQgCABBQgqPKA/vrx1Iugatz8FPi4BBQiqOqC3BBRIp/CALiWgQDoC+qDfk9o6rxtYQkAf82sEFD4AAc1pYy8HWEZAc9rYywGWEdCcNvZygGUENKeNvRxgGQFd9NPbl2pKAb8noEt+uATJphXwOwUFdOYDPRrrXYm0cgmD7hhBYBkBXfLDJbhjBIFlCgpoc+H7DAG9dMcIAsuUFNB2HXTRx8DfEFAgnaIC2hR04ecYX3EUHkinrIA2W/F//n3HzwsokE5hAX057HbPd/y4wADplBbQ00b8PaugAgqkU1pA71wFFVAgneICeh8BBdIRUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABJYVdL/frgFUJKAnsdgrKRySg3G8Mp4LysQgod5tk0wTmQxFQ7jadqqYwH4mAcjcB5aMSUO4moHxUAsrdBJSPSkC5m4DyUQkodxNQPioB5W5OY+KjKjSgP782V738+ffSn/P3/RBOpOeDKiigp2gOxdwPFw5+XjgMf9+P4VJOPqYiA7o//70uXAn1B/4g+smHVGJAj6c/1Kf+KwvXQf2FA+kUGNDTv3/81X/tsPv0fckwBBRIp8CAHvv1z8avb8tWQQUUSKfAgF6sdR6W7QUVUCCdAgO6nwb0eN6cfw8BBdIpMqCTZgookE2BAT1YAwU2ocCAHne75/GLe/tAgVwKDOj0yPv0iPx7CCiQTlkB3bWn0I/VPKV02Ra8gAIJlRfQVpvN9orO59/+3JSAAukUFNCXSUPbgB4W91NAgYTKCmjnx5fhcPzSmzEJKJBQiQG9g4AC6QgoQJCAAgQJKEBQ4QH98eWtE0F3M1Z8cUDlBBQgqOqA3hJQIJ3CA7qUgALpCChAkIACBAkoQFB5Ad2PB9QXXwkvoEBKhQX0cHVS0tKGCiiQTlEB/fHl5rTOZfdTFlAgoZIC2tyAflrM9u6giz4SSUAj5q5HuE/uMYJESgro8SaXTVIX3VLZ3+4ofRbTyz2N4DdKCuj+doP9tBK6aDeov8lR7ji+R+5pBL9RUEBPq5u3n8B58LHGMbnb+D65pxK8raCAntY2bzfXj66Fj8mdxvfJPZXgbQL6UeVu43vknkbwGwUF1CZ8Urnj+B65pxH8RkEBdRAJ2JaSAjp/GtPtWukbBBRIp6SAzp5Iv+xSJAEF0ikpoF0xL336vmgIAgqkU1RAp7di6riZCJBPYQF9cTs7YDPKC+hdBBRIR0ABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIOixAf35dfe6P/5K/avf8eIEFEhGQAGCBBQg6OEBfX7l244CChROQAGCBBQgyGlMAEECChC0bkB/fNnt/vw79W9cQECBdFYK6M+vTTebfu52n76n/pXvJ6BAOusE9Nie9fnrW74TQHsCCqSzSkCbNc/Teufpn1M7D7vdU+rf+W4CCqSzSkBPzfzc/dOkc98+yENAgXTWCOhp073ZA3r6p939mecM0I6AAumsEdD+dPrTP205BRSow4oBPXYb8gIKVGLFgO533XWdh4ynggookM5K+0CfzrtATzl1EAmowSpH4fdNOg/dRUjNKU2v3WDk8QQUSGe180Abz+12fM6LOQUUSGedK5EObT+bcB6yXogkoEBCK10Lf9yNh+Dzbb+/CCiQktvZAQQJKECQj/QACBJQgCABBQh6eEBfJ6BA2QQUIEhAAYKcxgQQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQtFZAf33r78GU8UPhXwQUSGmVgF7d1S5jQwUUSGeFgM7cFDRbQgUUSOfxAd1fJrPflv+c+te+j4AC6Tw6oD++nGp5/cFyh9PXPn1P/YvfQ0CBdB4c0OMrpWy6+trndT6SgALp+FhjgKAHB/Rf3zha9O8CChTNifQAQQIKELRaQJuzQT99fzm4EgmoxUoB7U6mbwKa5+j7QECBdNYJaH8x0img+5nTQlckoEA6qwS0ufroz79/fGlOCT0VNMf5Sz0BBdJZJaDH3e6pOXm+Paf+0D7IRECBdFYJ6L69FL4P6Gl1NN+BJAEF0lkjoKdkNuucfUBPq6D5tuEFFEhnjYD2F3QOAc1zEWdHQIF0BBQgKMMm/N4+UKAKKx1EatY5+4Ce1kcz3U35RUCBlNY6jenzENDmnNB8Z9ILKJDOOlcitWfPtwE9+FA5oBarXsrZyXghkoACCa10M5HxY+Hz9lNAgYRWu51dn9Cs+RRQICU3VAYIElCAIAEFCFopoMfdbhPH4QUUSGe1GyoLKFCbVQJ62AkoUJ+VbiaS8fL3CwIKpLPS7ezyfpjxmYAC6awU0HyfgnRJQIF0Vrwf6BYIKJDOSgeRNrILVECBhFYJaP+ZHhsgoEA665xI/+PLRjbiBRRIZ6UrkU4FdR4oUBlXIgEEuRIJIMiVSABBK51Iv5F+CiiQkCuR+FB2yeUeI3JyJRIVSJ/F9HJPIx5hpYNIbibCI+WO43vknkY8wlqnMbkSicfJ3cb3yT2VeIB1TqTfzA3tLMVVyp3G98k9lXiAlQ4iXXIeKEll6eFiuacSDyCgVCBLEBfKPY14BAEFCPK58ABBAgoQJKAAQQIKEPTYgDaHj54dRALqJKAAQQIKEGQfKECQgAIECShA0DqXcn76Pnm4z3hnJgEF0skS0PhBpPbz5Sf9/fVt2cAEFEhn/YD+/BoP6L4/jD9+Rp2AAvk8OKDXnwjfiW7C72+GIKBAPg8O6M0poK3gR8w12+/NR4MczgUVUCCfR2/CH9P187QC2u8KOGWzL6iAAvlkOIgUdYrluO9z3xdUQIF8Cgroz6+Tz/bsCyqgQD4FnUh/EdCmoJ8FFMjpwTcT+dc3jrf/+8KTmS4D2uwHfRJQIKOH343p+ZVvOy4+G3S6D7R7uHsWUCCfB2/CH3e72f2fzRlJiw/GH3YXOT4N49P/E1Agm0fvA22vvbxeC21O5AwcV7qu7nH5zUUFFEjn8QeR9pcXHzVb3tOLMZc4XP3kUUCBjFY4Cj9zNVL0Ws7D1Wn4zTqpgAKZrHIa01VC77ib3Wn99XLX6UFAgVzWOg+033K/q54JCCiQTkEn0qcgoEA6AgoQJKAAQYUH9MeXtw4izd2LdMUXB1ROQAGC1rmdXbITQa+9HdCZFyegQDKZAhq8FOluAgqkky2geQoqoEA66+wDbW6c1N08pInp0/QL6xJQIJ1VAnrK5fkCzGObzlNIc6yCCiiQzioB3V/Esnt0iB5I2t+zE0BAgXRW2gc6vSNody/65Xekbxzu3JEqoEA6RQW0vTvzpYUDEVAgnTUCevVhRvt24z0Q0PaOTpOfao/uL9sRIKBAOivtA704iNTmdL98H+jxJpftB8stenECCiSz1lH4Yc2xWWtsjsIfAh8qt79daV16MF9AgXTWOQ/08tjPc1vUxeeB3tyNvhvyohVZAQXSWelmIpPDP004m0+DW7wCOvsh8wt3pQookM5qd2Pqr+fscneMXIckoMC2FHQ7O5vwwLYUFFAHkYBtybMJHzN/GtOifakCCqSzUkAnd7SLJ3T2RPplgxNQIJ11AnpxR9D4fexmbiy6cGACCqSzSkDbVcfuCPpx8dWXF/ZX/XQzESCfVQI63Xm5+OrLa25nB2zEStfCX+65zPN5SA0BBdJZ6W5M00Pl4VspJyCgQDrZ7geah4AC6QgoQJBNeIAgB5EAgso7jekuAgqkU9iJ9PcSUCCdDJdy5juEJKBASiXdTCQBAQXSKel2dgkIKJBOSTdUTkBAgXQEFCBIQAGCHhvQmTsg590RKqBAOgIKECSgAEH2gQIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIBCzXa93K+jUgIKFdvtFPSRBBTqNYZTQR9DQKFak2xa9B9CQKFa0+Xdsv8IAgrVEtBHE1ColoA+moBCtQT00QQUqiWgjyagUC0BfTQBhWo5jenRBBTq5UT6BxNQqJhLOR9LQKFm+vlQAgrbsitA7mm0GQIKm5K7je+TeypthYDCpuRO4/vknkpbIaCwKbnT+D65p9JWCChsSu40vk/uqbQVAgqbkjuN75N7Km2FgMKm5E7j++SeSlshoFA5S/3jCCiUyopldgIKhUrfT38eSwkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBBUU0J9fZz+E4I+/lrw4AQWSEVCAoIIC+vLji4ACG1JSQNt10Ke7hiCgQDpFBbQp6Kfv9wxAQIF0ygposxX/5993/LyAAukUFtCXw273fMePCyiQTmkBPW3E37MKKqBAOqUF9M5VUAEF0ikuoPcRUCAdAQUIElCAIAEFCBJQgCABBQgqPKA/vrx1M5G5e4+s+OKAygkoQFDVAb0loEA6hQd0KQEF0hFQgCABBQgSUICg8gK6Hw+of17+wwIKpFNYQA9XJyUtbaiAAukUFdCZj+VcdBKTgAIplRTQX98ui9l+UPyy+9MLKJBOSQE93uSySeqi+9MLKJBOSQHd326wn1ZCF+0GFVAgnYICelrdfLr54mHZNryAAukUFNDT2ubt5vrRtfBALgIKEFRQQG3CA9tSUEAdRAK2paSAzp/GdLtW+gYBBdIpKaCzJ9IvuxRJQIF0SgpoV8xLn74vGoKAAukUFdDprZjcTATIrbCAvridHbAZ5QX0LgIKpCOgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgQBl2vdyvY0JAgSLsdtsrqIACJRjDuaWCCihQgEk2N/RnLKBAAaZ/u9v5OxZQoAACugHbmfDAEgK6AduZ8MASAroB25nwwBICugHbmfDAEgK6AduZ8MASTmPagA1NeWAJJ9Lnt6EpDyziUs7stjTpgUW2108BBYgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFyuBSztw2Ne2BBdxMJLstTXpgAbezy29DUx5YwA2VN2BDUx5YwEd6bMB2JjywhIBuwHYmPLCEgG7AdiY8sISAbsB2JjywhIBuwHYmPLCEgG7AdiY8sITTmDZgQ1MeWMKJ9PltaMrDB7DbvjvHL9GEOg8x9QBTElBYU+46vsOd45doQp2HmHqAKQkorCh3HN/lvhFMNanGIaYeYEoCCivK3cZ3uW8EU02qcYipB5iSgMKKcrfxXe4bwVSTahxi6gGmJKCwotxtfJf7RjDVpBqHmHqAKQkorCh3G9/lvhFMNanGIaYeYEoCCpu0rS4ueNnJh5h6gCkJKJCOgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShA0IcLKEA6yRuVeoAp5Z7YQF2SNyr1ACtW1R4FI7NVNY1MTeMyr/oRTKiqpcHIbFVNI1PTuMyrfgQTqmppMDJbVdPI1DQu86ofwYSqWhqMzFbVNDI1jcu86kcwoaqWBiOzVTWNTE3jMq/6EUyoqqXByGxVTSNT07jMq34EE6pqaTAyW1XTyNQ0LvOqH8GEqloajMxW1TQyNY3LvOpHMKGqlgYjs1U1jUxN4zKv+hFMqKqlwchsVU0jU9O4zKt+BBOqamkwMltV08jUNC7zqh/BhKpaGozMVtU0MjWNy7zqRzChqpYGI7NVNY1MTeMyr/oRTKiqpcHIbFVNI1PTuMyrfgQTqmppMDJbVdPI1DQu86ofwYSqWhqMzFbVNDI1jcu86kcwoaqWBiOzVTWNTE3jMq/6EUyoqqXByGxVTSNT07jMq34EE6pqaTAyW1XTyNQ0LvOqH0GARxFQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQiqMaC/vu3++Kv9v5//1j/+9P3223582TX+/PvVAZ2+4/P8M4fd7mny4PMrTw3PdL9q+iJO3/V8fvTza/P88GOn1zt5riblzZnL56udM+XNmOGbnofXn2nGVB3QfTdfZpeG427ULzq33rc0/PhyMYDpU/t2tnZ5vPhNzZfOg973Tw+v8mqI9Shtztw8X+ucKW3GnL/0PDvE9dQc0OPu9aVhsjDsZt9sG+9aGq7f/CZP/fzavI7zwnBeHI7T33q4efrw1pt8wQqbMzPPVzpnCpsxvWbF43kcRp4ZU2NAB28sDc0c6r/W5OuVaf+upeF6zk2eOrZPNb+g/Uozv7vhnV7RP34dvq3ZWGm+7zg+3by8p5eKFTJnZp6vfM4UMmPGF3sOaK4Z80EDephsDjQBm9+B8p6lYdwNc/tUtz002c44DO+nzQ8dhvfW/e6822d4nYfXN5JqUMacmXu+8jlTxozpdKupz+eBZJkxHzSg++lXDjc7V3rvWRr21+/F56d+fWt+yfE89NML6WZ480Ongbffd1oOzlvuw49+3BWdDc2ZuecrnzNlzJjhmT/+7ySgmWZMyQEdZ+lp2vX/186/bofOsMvmuV8a2sfDnLtYGvqZPlloTl9pvrN9on2ju14oxll+MU8vnzo912+OTPaQtzO8/aHTr2tfzWSR639tN5SSV3TqmDMzz78UPmfqmDHDDz1fHJjPM2NKDug4+Y7jqvyhn/u3S8N+fNh/380b5fzS0J+5cfV+PM7iw9ymzrgL7fp39Ath90PH7sVc7lwfloHbpawkdcyZmedfCp8z9cyYdjZcBDTPjCk5oOO622F4w+uWg7ml4Z8vD+o18/h6Ns4uDf84HA+8/O5hlp/ea6/3po9Lw+25af13n//5/HJxetvkFcwMuCB1zJnb518ZcDmqmTHdyuhFQPPMmJIDOqywNftJumnXTeCZczJ249peP6v2t/N4dmm4OUDeGWb5zLveuDR0Z2TMPDW8j3ZvrvvJQjDZTNrfLK0FqWPOzP3o3JMFqWbGdP9zeW59lhlTckCHjYHTMvC/uv/rpvPc0jAeQx3m6u3Jl/NLw80B8uGnn/p/rzf2xqXheL0a0w913FDvVgems33y/zNDLkgVc+bm+XEw5c6ZSmZM/zovh5VlxhQd0P7N7Lj74z+6MyO6AM0tDeczKMa/hP24PHRfml0ahq/sdxfvm8Msn3nTO5zfsS/faU9Da1/W+NfZvc7LgI6LwHF380ZdkCrmzM3znaLnTB0zZthav0xmlhlTdED7XSL73Z//PWyITN+vZs7JuFiVOF9cMWzW3C4Nw3cfLzdI+ll+9Yc2fer04i4XlOMrF8C9sgb6+vkgJahizrzyfNFzpo4Zs++/fhnQLDOm6IB2E/I0S566c8v6t9f3Lw2NdsPkjXMyxq8sXBq6MzJGbywMFQa0ijnzyvNlz5kaZszh4gyCkYAu1u7AaS9saPft9Dt4li0N7WzoN2PeXBqmP/iOpeFwsT2xf/Xit/mj8POvtBw1zJlXni97zlQwY86/4yag68+YsgPaXn1wbOZIu0j0lzi8Y2n48WW6u3k/7EhP+HZ6cUZGs0vptTfHyZ6b6eDKXs+pYs688nzZc6aCGXM+mNUZ/nysgS7XzP92GWj27Qxnif1+aZielfIyvJH9fofO+Uyz/q3vvDTcPDW9v1ZzaPLV3duvvGeX/WdaxZx55fnC50z5M0ZA0zns/vjPb+0UPC0T/9Wv0P9+abjau9LvejzvgTzshqVh+Lb2Hfd82fpwuHx6adzlU5MzMq5P6Lg0ey184cd6X6qYM688X/icKX/GvBZQR+GXO03o//2lP5/t0//sZ8jvl4bmYt3z7pJh2RhvizQ8fT6prfuWX9/Oy1s3xHEnzM1T++mNtt78K527G1PhZxu+1DFn5p8vfM7UMGMGzgO9V3vFxHCq2OTWgcPSMOzgud6hc9ztprN+PMV3uOzhvDQMl1uO/3ZbLv3C9GO8rOLqqfMZGeNy8poft/cDfSn7epdGBXPmlecLnzMVzJiRK5HuNr2md/LWdl5AxlvLtMZdNJebAZPlovV/vg5Lw3hh72SIu930V/Vvy1dPHae7gqZmFozbKzwKv+K6Vf6cmX+++DlT/oyZjopr4e/TTOvxUt2+P8PS0C4gu6f5czKmM2myE6V7/HNcGoZby/SDPl7/yPnWMpdPjWdkdC/i7aXh+jORXorf0fZSwZx55fni50zxM+bsIqB5ZkzpAT3vLTmO+2jGA33tnPj82kltw7vn1dkZzfdOlobJ9so4yMk73XmD5OKp7saw09/y5tJw+amc3Qt58+TuAhQ/Z155vvg5U/yMObsIaJ4ZU3pAs7u5vXYSdd/3fB3mzEZVNWME9E4/rj/gJYmib3u+EebMRlU1YwT0Xo/4PFWrOSmYMxtV04wR0Hud5lzy99NKP318ZebMRtU0YwT0btMr0FINsegzDTfDnNmoimaMgN7v9pOw7vPr5pNhiDFnNqqeGSOgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSg5HLYXfr82F/369vuj78e+yv4cASUXASU4gkouQgoxRNQcjkF9Dn3a4C7CCi5CCjFE1ByEVCKJ6DkMhPQH1/6/ZS/vu0+fW8efn75+XWyf3S/+/Pv4+nx6cmTY7vv9DyQ49W+1IvHk32gFz/Xfb3dIfvn3+nHkqoJKLnMrYGevvZ0/rcJ6Om/3VjMJqD/NTzsn9kNXRwfd4O4fjwGtC3y+eear//H8K1WiVlEQMllLqDdmmfTuGZt8JTAf/rWp60r6H73D1/6tcpzH7sSjl3sO3j9eAjo5Ovt41/fzt83dBreR0DJ5eo0pm77+djGcd+VrG1kU7njsCW+H7+x6V77TcNzh0kQm2+5fjwEdN+vkhCJ08UAAALnSURBVB76n+sC+tx/5SnDhKBcAkouswFt23nsc9kEtPvy6f/aWu7H7zuO/3dapWye2w8rtH0prx/3/wxDGv+vCWi/E+Dw8JNRqYyAkst8QE85/Kdhb+U5dk0On87/vIwrqY1j+8XD1UGg68d9QCeR7Ha0NgHth3n6fY4jsYSAkssrpzEdzgdzJkHrV0rH1cp+L+nwbZ+HQ+vnQV4/vl4xHX6u3+169fvgPQSUXF4JaLNG2K8jdoWb/u+43jk9RNSvvu4vDrrfPO4COsll30sBJU5AyeW1E+n3rwS0i+R5B+bN9v95n8DT8Asmj8eAjlfEd2uxAkqcgJLLKwE9XmzCv7YGOt+6cb306faxNVDSE1BymQ/oqWf/41/Ox95v9oGeN+FfubdSG81JB4fHt/tAu2EKKHECSi7zAW2+euxXIccrO8fujQFt9pRe/PRk27yN6/XjmaPw+2HNVEAJElBymQ1o27Chac2Ozs/Dl/uz4IfYHXbTS9ufpmc4dcG8fvzGeaACSpCAkstsQPfDZe5NONsjRU3SJlciDbFrt8yfu6+1Xxy/p3nm8+3juSuR+iuUBJQgASWX6zvSNzk8ni/ZfG6D9g+Xdwy5PH3+8qj7/vx4LOXk8RvXwgsoQQJKLjMBnd4x6ZSyZkX0OO3nNKA3d1863xWk+56rx2/cjUlACRJQcpkJ6GHcb9n+X7sl33Zy2NafBnQYwvny9T6Nz7OP37gfqIASJKBs1+Q8UNgiAWW7BJSNE1C2S0DZOAFluwSUjRNQtktA2TgBZbsElI0TUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUICg/w/gSJTV2cgt/AAAAABJRU5ErkJggg==)
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")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAhFBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZmYAZpAAZrY6AAA6OgA6kNtmAABmADpmOgBmOpBmZjpmZmZmkJBmtrZmtttmtv+QOgCQkDqQkNuQ29uQ2/+2ZgC2Zma2tma225C2/7a2/9u2///bkDrbtmbb/7bb/9vb////tmb/25D//7b//9v////zSYFjAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3dDXfaapuYUZKZM+OUdKZvW6dtmM708LYHJ/z//1f0gZBAYOv2Iz2S2HutrBgwsvXhC0lIYnMEIGST+xcAWCoBBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAoHkHdAOQTPpEJR9iQrmnNrAuyRuVeoApjfCCATwtAQUIElCAIAEFCBJQgCABBQgSUICgVQd0iqO2gOcloABBqw7oLQEF0hFQgCABBQgSUIAgAQUIEtBsHBcASyeguTiyChZPQDNpwqmgsFgCmkcrm7P5nYCBBDSP9i8ym18KGEZA8xBQWAEBzUNAYQUENA8BhRUQ0DwEFFZAQPMQUFgBAc3DYUywAgKaiQPpYfkEdNCz5y/VlALeJ6CDnj1/qaYU8D4BHfLkJUg2rYD3COiQJy9BsmkFvEdABz17/lJNKeB9Ajro2fOXakoB7xNQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUICgxQb0sDl5GfosAQXSWVhA96dqfv3zePz9Y1N5HfZ8AQXSWVRA62x++XncbTahggookM6iAnrO5h9/P/3761iukBbrox8noEA6Swroocrm27dTQ7fNXYNWQQUUSGdJAd1Vq51FQZv1zt2wN5IEFEhnQQH9/aNe7zx90VRzX0f1gwQUSGdBAf31vd5cb0p6LLbhB+0EFVAgHQEFCFpQQG3CA/OyoIA2byIdLgcvtVv6EQIKpLOkgNaHMZUncW6buxzGBGSypIA25282B9IfHEgP5LOkgBZvI1WncjYpdSonkM+iAnq5mEidUhcTATJaWEAv9i5nB2S22IDGCCiQjoACBAkoQJCAAgQtPKDtK9vd2vSY8JcDVk5AAYJWHdBbAgqks/CADiWgQDoCChAkoABBAgoQtLyA7po31AefCS+gQEoLC+j+6qCkoQ0VUCCdRQX07dvNYZ3DrqcsoEBCSwpoeRnlVjHLi4IO+kw5AQUSWlJADze5LJLqM5GATJYU0N3tBvtpJdSncgKZLCigzefCt/lceCCbBQX0tLZ5u7l+cC48kIuAAgQtKKA24YF5WVBAvYkEzMuSAtp/GNN2yCAEFEhnSQHtPZB+2KlIAgqks6SAVsXs+vJz0BAEFEhnUQFtX4rJxUSA3BYW0KPL2QGzsbyAfoqAAukIKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkICSwqaW+/eASQkoCWw2CsozElA+rwmngvJcFhbQX9+//ll+sS/Wdr78HPp8f99jaGXTBOapLCugu82mDOjhvMFY5/TD/H2PoT1VTWGeyaICuqvXOg+XXW4DV0L9eY9BQHlWSwro27fN5uVYbMfXq56/fwxdB/XnPQYB5VktKaCnFdDX4v991dHCaV10O2QQ/rzHIKA8qwUF9LS++VL/f1nt3G3++GvAMPx5j0FAeVYLCuhpy/21/v8Szf2wbXh/3mMQUJ7VMgP60tx5END8HMbEs1pQQE+b7tv6fwGdFwfS86QWFNDjrg7nrrMJbx/oDDiVk+e0pIAe6qM+375V2/LHy1rpR/kDH4l+8pSWFNDisM9yHfTyztFu4JH0/sKBdJYU0PJI+nKTve5mcUbSdtAQBBRIZ1EBrQra8fL+k9oEFEhnWQGtTodveX3/GR0CCqSztIBWZ8JXBr3/XhFQIJ3lBfRTBBRIR0ABggQUIEhAAYIWHtC3b4/Ohb8+5smpMkBKAgoQtOqA3hJQIJ2FB3QoAQXSEVCAIAEFCBJQgKDlBfRyOZGBV2IqCCiQzsICuv/c1ewEdCyOEuMpLSqgt5cDHXYQk4COxXG2PKclBbT4SI92McsL2w27pp0/8FE04VRQnsuSAnq4yWWR1EHXVPb3PYZWNk1gnsqSArq73WA/rYQO2g3q73sM7alqCvNMFhTQ3o8w9rnwMyCgPKsFBfS0tnm7uX5wLnx+AsqzElA+TUB5VgsKqE34uRJQntWCAupNpKT6LpY6N7mnEbxjSQHtP4xpO2QQ/iYbueP4EbmnEbxjSQHtPZB+2KlI/ibPcrfxY3JPJXhsSQGtitn15eegIfiTPMtQw4DcUwkeW1RA25diqriYSFSWHg6WeyrBYwsL6NHl7FLJmMUPyz2N4B3LC+in+JsMEEa4Q0B5R/p+mgmshYACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShA0LgB/fV9c9/XP1P/6A/8cgIKJCOgAEECChA0ekBf73zbQUCBhRNQgCABBQhyGBNAkIACBE0b0Ldvm80ff6X+iQMIKJDORAH99b3oZtHPzebLz9Q/8uMEFEhnmoAeyqM+f//IdwBoTUCBdCYJaLHmeVrvPP13aud+s9mm/pkfJqBAOpME9NTMl+q/7em/XXkjDwEF0pkioKdN92IP6Om/cvdnniNAKwIKpDNFQOvD6U//leUUUGAdJgzoodqQF1BgJSYM6G5Tnde5z3goqIAC6Uy0D3R72QV6yqk3kYA1mORd+F2Rzn11ElJxSNO9C4yMT0CBdCY7DrTwWm7H5zyZU0CBdKY5E2lf9rMI5z7riUgCCiQ00bnwh03zFny+7fejgAIpuZwdQJCAAgT5SA+AIAEFCBJQgKDRA3qfgALLJqAAQQIKEOQwJoAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUICgqQL6+0d9DaaMHwp/FFAgpUkCenVVu4wNFVAgnQkC2nNR0GwJFVAgnfEDuusms96Wf0n9Yz9GQIF0xg7o27dTLa8/WG5/uu/Lz9Q/+CMEFEhn5IAe7pSy6Oq9z+sck4AC6fhYY4CgkQP6rw/eLfo3AQUWzYH0AEECChA0WUCLo0G//DzunYkErMVEAa0Opi8Cmufd9zMBBdKZJqD1yUingO56DgudkIAC6UwS0OLsoz/+evtWHBJ6KmiO45dqAgqkM0lAD5vNtjh4vjymfl/eyERAgXQmCeiuPBW+DuhpdTTfG0kCCqQzRUBPydwem4CeVkHzbcMLKJDOFAGtT+g8BzTPSZwVAQXSWW5AT0PdDn6SgALpZNiE36XZByqgQGYTvYlUrHPWAT2FL8nVlAUUyGyqw5hezgEtjgmNHUnf89EghUH7AwQUSGeaM5HKo+fLgO7jn4gkoMC8THoqZ6B4bQcBBeZkoouJNB8L/6kTOYsON6uv9oECmU12Obs6oZ88gml3+YwlAQUyW9oFlYtPo9uWXwkokNnSAlpf2ekooEB2iwto9anyrwIKZDdRQA/xN85vFZvxLwIK5DbZBZUTBrQc3tf/EFAgr0kCuv/MoZu9qiFuh/9yAgokM9HFRNKc/t5SbMYLKJDVRJezG+Ea9DsBBfKaKKDb1D8lSECBdCa8HugcCCiQzkRvIqXeBRoloEA6kwS0/kyPGRBQIJ1pDqR/+zbSRvxpwA+Oieq7+t0ovwbwlCY6E6k86ijhcaCt4QookMkiz0S6eBzQnl9OQIFkFnomUpSAAuks9UykIAEF0pnoQPqZ9FNAgYSciQQQtLwzkXbNrtTAaq2AAulM9CZSqouJXL8dNbShAgqkM9VhTEnOROoeTRp5Q19AgXSmOZA+zQXtysNJW8UsPiZ+4HAFFEhnojeRkhwHerjJZZHUQeu2Agqks6SA7m6fOPQIKQEF0llQQHvfzB/4/pSAAuks6HPhey+Kd3AuPJCLgAIELSigNuGBeVlQQL2JBMzLuAEt3j56Hfkwpu2gX05AgWSWFNDeA+mHDUtAgXSWFNDbAW02X34O++UEFEhmSftAj+1LMbmYCJDbwgJ6dDk7YDaWF9BPEVAgnWlO5ezsqdwluzrocAIKpJMloD6VE1iD6QP667uAAqswckCvP4KjYhMeWIORA9pz5OZm4MlDSQkokM7Ym/CHWfVTQIGEMryJlJOAAukIKECQA+kBgka+mMi/Pni//d8yHMwkoEA6o1+N6d6nDg/8LI5EBBRIZ+RN+MOdC869fcvzZryAAumMvQ+0KOXNWuh++IU8ExFQIJ3x30TadU8+Ki8rH7oUXQoCCqQzwbvwPWcjZTuXU0CBdCY5jOkqoflOhRdQIKGpjgOtt9yz1vMooEBKDqQHCBJQgCABBQgSUIAgAQUImuZydrM5EFRAgXQyBTTXqUgCCqSTLaB5CiqgQDrT7AN9+3a+eEgR0237jmkJKJDOJAE95XLb3DiU6TyFNMcqqIAC6UwS0F0nltWtfZY3kgQUSGeifaDtK4JW16J3RXpg6QQUIGiKgP7+cbUJX2y8CyiwdBPtA+28iVTmdGcfKLBwU70Lf17fLA5jKt6F3/tQOWDppjkOdN85hP61LKrjQIGFm+hiIuWHc1aKcJ4247OsgAookNBkV2Oqz+estuQPmT7VWECBhFzODiBIQAGC8mzCZyOgQDoTBbR1RbusCRVQIJ1pAtq5Imim949KAgqkM0lAf/8oj/4sHDJ+nsdRQIGUJgnooRXNIqavt98yEQEF0pnoXPjWjs9Ml1KuCCiQzkRXY9q2bua5lHJFQIF0sl0PNA8BBdIRUIAgm/AAQd5EAghyGBNAkAPpAYIynMqZ82R4AQXScTERgCCXswMIckFlgCABBQgSUICgcQPaefd9Du/DCyiQjoACBAkoQJB9oDyV3lf0T8k9RuQkoDyT9P20RD01AWUFRuhicrmnEWMQUJYvdxs/JvdUYgQCyvLlTuPH5J5KjEBAWb7cafyY3FOJEQgoy5c7jR+TeyoxAgFl+XKn8WNyTyVGIKAsX+40fkzuqcQIBJSnooukJKA8k/T9tEQ9NQEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElBYs00t9++xUgIKK7bZKOiYBBTWqwmngo5DQGG1Wtm06I9CQGG1Tst7swVv2R+DgMJqbTpy/zZrJKCwWpv2PlDL/ggEFFarsw/Usj8CAYXVai/vlv0xCCisloCOTUBhtWzCj01AYbW8iTQ2AYXVchjT2AQUVsuB9GMTUFgtp3KOTUBhvVxMZGQCCitmD+i4BBTWTD9HJaAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKELSwgL5922w2f/zV3P79Y/P1zwHPF1Bmb7MAuafRbCwroLt69r2c7xBQ1iZrGD8s66TJ8sP7LSqgu2b2nVdCBZS1WcIymud3zJzvXksKaLH9/nr6f38pqICyNlOuR8blmTBXX8zAkgK623z5WX5xymZdUAFlbXKn8WOyTJfLl9P/+DsWFNBTLJt9n7u6oALK2uRO48fkmC79X+e1oID++l5uwFfqggooa5M7jR+TY7r0f53XUgNaFPRFQFmh3G38iCyTpf/rvBYb0GI/6FZAWZ/ccfyILJOl/+u8FhTQ9j7Q6ubmVUB5XusI40cJ6GftN+1V0OKopi//S0DhKQjoZxXHgW5btw/la6aAwhNorR/P6M94SQEtj6Bvb8YfBBSeRVPQOe1qWFRAy4JuW7eLdVIBhacwx121ywpo8c7RtnPHXkDhScyvn0sL6GfNatoDCyegAEECChAkoABBCw/o27dHbyIt7FQLYGEEFCBo1QG9JaBAOgsP6FACCqQjoABBAgoQJKAAQcsL6OXD4V/e/+ZrAgqks7CA7q8OShraUAEF0llUQIur110ZdBCTgAIpLSmgxacgtYv56/um/nj4DxNQIJ0lBfRwk8vyg+WGDEJAgXSWFNDd7Qb7aSV00G5QAQXSWVBAb65GX9gP24YXUFisGV7QYkEBPa1t3m6uH5wLD89hjpcEElBgCXwq5+fYhIfn5XPhP8ubSPC02n+78/k7XlJA+w9j2g4ZxHwmPDCEgH5W74H0w05Fms+EB4YQ0E8ri9n15eegIcxnwgNDCGgCu6t+upgIPAcBTcPl7OAJCegMzGfCA0M4jGkGZjTlgSEcSJ/fjKY8MIhTObOb06QHBplfPwUUIEpAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAgdFs5u+T45doQl2GmHqAKQkoTCh3HD/kcyOYalI1Q0w9wJQEFCaUu40f8rkRTDWpmiGmHmBKAgoTyt3GD/ncCKaaVM0QUw8wJQGFCeVu44d8bgRTTapmiKkHmJKAwoRyt/FDPjeCqSZVM8TUA0xJQGFCudv4IZ8bwVSTqhli6gGmJKAwpdxx/IBPjl+iCXUZYuoBpiSgQDoCChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBTxdQgHSSNyr1AFPKPbGBdUneqNQDXLFV7VEwMnO1ppFZ07j0W/0IJrSqpcHIzNWaRmZN49Jv9SOY0KqWBiMzV2samTWNS7/Vj2BCq1oajMxcrWlk1jQu/VY/ggmtamkwMnO1ppFZ07j0W/0IJrSqpcHIzNWaRmZN49Jv9SOY0KqWBiMzV2samTWNS7/Vj2BCq1oajMxcrWlk1jQu/VY/ggmtamkwMnO1ppFZ07j0W/0IJrSqpcHIzNWaRmZN49Jv9SOY0KqWBiMzV2samTWNS7/Vj2BCq1oajMxcrWlk1jQu/VY/ggmtamkwMnO1ppFZ07j0W/0IJrSqpcHIzNWaRmZN49Jv9SOY0KqWBiMzV2samTWNS7/Vj2BCq1oajMxcrWlk1jQu/VY/ggmtamkwMnO1ppFZ07j0W/0IJrSqpcHIzNWaRmZN49Jv9SOY0KqWBiMzV2samTWNS7/VjyDAWAQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIWEdDfPzZf/yy/+vU/69tfft5+29u3TeGPv+4O6PQdL9179qcndIe1O93zUj/02tz763vr+/blNzx4asvpia9Xd90d8L5+djUi7WF3nlI+Z7PZ1rdOU+P6J0zGnOnOme7jGeeMGdOdMedvqu5KOGOWFdBdVY3epeGwadSLzq3+paFJUaGM083S0FkY3r6VP+DBUzt33vyqdwe8K++v8tgZj+Kuy4B39cPn59W/Tw7mTHvO3Dyeb86YMd9vBlzc9dr5fRJYVEAPm/tLQ2thuH6Vu7izNLRfgMvhXC8NnXl2fvl68NTO7/5oaegM+Nf3YiwvC8NlcThsOi/m1w/vH61CjMqcaY1Tz+PZ5owZcztOxYrHazPARDNmEQE9e7A0tCbrfnN3m+TO0tAa2GnQt0tDceflFes87R88tfWt/+/H1z/33Rf4ewM+lAPen1+kL1s3p2/7p+/nV+5iY6X4vsPlh/363nlZz8CcOfY+nn3OmDHtSdHJcOfBsLUEdN+aD0Vi+ndx9C4Nf2tvVZy+42/XS0N3njU7Uh48tVYtH5v2jHsw4Gprq7WdsT8/WvzIZonanX/G6d7LnqBsG/EVc6b/8exzxow5q1ZTXy9DTDNj1hLQXfue/e1mQaV3afhv31uvvqfp+r+vloabeVZ/94OnNuqNi+6ddwZcbbkcLt9dPP56/pGnX31bD/Gy5X5eFue8ovM0c6bv8exzxoy5PPL1f3T3pF4e/IRZBbSZpaexq78q51+1Q+fQvDJVS8OhvUOlszTUM7210JzuKb6zfKCz47qYM/vOzpqXQ3dpuJpnl7ny4KkXxbN7Xvp7Bnwacr05sr1Mj/Ibyx95+u5yXFsLdD1S1SDHXdExZ96fMz2PH0efM2bM+zOm9aNbq7WJZsysAtqM4KFZ2d7Xc/92adg1N+vvu3kF7V8a6iM3Lnt/XluvVMWaf3dpuJ5nl/n/4KntX+Flf7049A+4OSLjol7EqwEcqlFtLS2XQ1WOnRfbEZgzXX1zpufx4+hzxozpujdjyp/aCWiiGTOrgDZrV/vzC161HPQtDf/pvKuk2edx87LVuzT80/kdu+qhYpq2QlS8LHWWhut59ut7a63v7lMbu9Odre+5PPF2Ybg9Nq3+WZf/Xo7H9kLQGr9f3++9CZCGOdPRO2duH7/6cgxmTMe9GVOtjHYCmmjGzCqg5+lW7kkux66aBD3HZGya9bH6dWTXnsfN4G6Xhqu3sMtp2t5F/XJsLw03mxOt160HT73YnR799Z9vX05vt1OqIzKuvm97vLyOVi+uu9ZC0NoI2938LSRlzlx/3/Z4PWduHz/2PpiUGXP9fdtjz4ypvugENNGMmVVAzyv7p2Xgb9VXh/Oq+83S0LzLeZ7++/Pra/Mq1b80dN/CLqdp8zJevl3YWhr++4/28I7Hzjx48NT3xrJnwKfnXj21/p2bV+PqR7Vne+vr7rKRnjnT0j9nbh6/+b3GYMa03Jsx9VTozoo0M2ZeAa1frA6br/9eHflQJaJvadg2z2iW1V2zPFR39S4N53t2rWnaTPDyeLXW0lBrzabrZvU/9bHeAZ9PGWlPi3LIzfJW/ahuQJtF4LC5enpi5kxnWvTNmZvHKyPPGTOmMy16Z8x5a72bzDQzZl4BrXda7DZ//N/zhkj7FaXnmIzOi/3l5IrzZs3t0nD+7nrGVdP08ibf9thdGv7467BpTej2n8qDpz7WN+Cb89cOd06vu7MGenu0SVrmTHtMHr5/e7h+p3jUOWPGtMekf8bsmmK3A5pmxswroNWollO2PHasfnn9+NJQ2J9fqPrfUjw291yWhvrG27f2YWX7+nV51z0n7Wpp6H3qY30Dro/IaDxYGLIE1Jw5G9jP0eeMGVO7O2Oabj5BQA/1TpLXat9OvYNn2NJQTqh6M+bh0lAfQ3bZY7Q77y0/Lw3nvfLNz7hdGnqf+ljfgE93blvfsrt7al3/u/D90yElc6Zyf87ceXzsOWPGVO7OmMsY3AQ0wYyZWUDL8wMOxfQtF4n6FIYPLA1v39o7hHfnGfXBl9Pqv/rkhKujgtv70XuWht6nPnZnwJcRKHZY3RtOa89N+5cZfQ3UnKlvPBxOz+MT7FwxYx7NmMv+08q2OzqfNLOAlvO/XAaKfTvn47jeXxr23R3C+/tLQ+8OnWpiHpq3MbtLQzkLmiMvrpeGvqc+1jfg9vW1ioVk2/fE4/F6gW5eQkcPqDlzfGfO3Hl8gp0rZsyjGfNcAd1vvv6fcudysVz8vZ5W7y8NV/s/6p2Dl32E9SZA6/3RXevFsP4J9Yv37dLQOgLt9sihvqe+N463A24dkdF637NH77nw478Lb84c35szdx4ffc6YMY9nzL2ArvFd+HJS/Ndv9UvTl/9Sz7v3l4bOvpFm2WguXHR++LIBcP6W85zZl5c52NYPXS0NrcG37nzw1Mf6Brxr3/Xwr7TvakzXu3fGYM68N2f6Hx99zpgx782YnqHc3IqaW0CLF5jmYK7ztG0tDecdPNc7dA6XY8TKKx9si6/2m+bEhMvScN4M6c7y8pHzK9vN0lAO4aX+vm1934OnPtYz4MsRGb9/vDNf326vB3oc/Uykoznz7py58/joc8aMee9Ppm8o6zwT6diaT+Vy0d2NUk73zWvf0nC1ot5aLkr1RbTeWif2dg8NK1DNHMgAAAYHSURBVH9a80p7szQcz1ez/n25wsuDp747itcDPrR3NLX1LBiXMb1sgo19Lvz555oz9+dM/+MTzBkz5r0/mb6hrPFc+MLlWNldU4jz0lBO99OjvcdktCfjtjWw8vavZmk4X1rmcl7b6/mL+qvepaFYsOoNmNZOontPfaxnwM0RGb9/bDr6lobrz0Q6TrEL1Jx5Z87ceXyCOWPGRAKaaMbMLqCX/RmHZh9N80Ze/cp156C286tnewruqtC0lobW9sqxNU0vb870Lg3NlvNb5+KGd5762O2A//FfzuPT/niXu0vDr86nclaj+fDg7iTMmYdz5s7jE8wZMyYS0EQzZnYBnb3d+BvLg+W+7vk8mDMzteYZI6BDvX178BqXSe5P3pkHc2am1jxjBHSwfB8ifI/VnIo5M1MrnjECOthp2s/s9XR+y2ce5sxMrXjGCOhw7XPI5qC6rA3mzGytd8YIaMDtZ1nl9Pvmk2GelzkzU6udMQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQJKLvtN18u4P+73j83XP8f9ETwdASUXAWXxBJRcBJTFE1ByOQX0NffvAJ8ioOQioCyegJKLgLJ4AkouPQF9+1bvp/z9Y/PlZ3Hz5fjre2v/6G7zx1+H0+3TgyeHct/pZSCHq32pndutfaCd51X3lztk//gr/ViyagJKLn1roKf7tpf/i4Ce/m2aYhYB/fv5Zv3I5tzF5nY1iOvbTUDLIl+eV9z/7+dvtUrMIAJKLn0BrdY8i8YVa4OnBP7zjzptVUF3m3/4Vq9VXvpYlbDpYt3B69vngLbuL2///nH5vnOn4WMElFyuDmOqtp8PZRx3VcnKRhaVO5y3xHfNNxbdK7/p/Ni+FcTiW65vnwO6q1dJ9/XzqoC+1vdsJ58MLJmAkktvQMt2HupcFgGt7j59VdZy13zfofnqtEpZPLY7r9DWpby+Xf93HlLzVRHQeifAfvSDUVkZASWX/oCecvjP572Vl9gVOdxe/js2K6mFQ3nn/upNoOvbdUBbkax2tBYBrYd5+nneR2IIASWXO4cx7S9v5rSCVq+UNquV9V7S87e9nN9avwzy+vb1iun5efVu16ufBx8hoORyJ6DFGmG9jlgVrv1ls97ZfouoXn3ddd50v7ldBbSVy7qXAkqcgJLLvQPpd3cCWkXysgPzZvv/sk9ge/4BrdtNQJsz4qu1WAElTkDJ5U5AD51N+HtroP2ta9ZLt7e3rYGSnoCSS39ATz37x3+5vPd+sw/0sgl/59pKZTRbHTzfvt0HWg1TQIkTUHLpD2hx76FehWzO7Gy61wS02FPaeXZr27yM6/Xtnnfhd+c1UwElSEDJpTegZcPOTSt2dL6c766Pgj/Hbr9pn9q+bR/hVAXz+vaD40AFlCABJZfegO7Op7kX4SzfKSqS1joT6Ry7csv8tbqvvLP5nuKRl9vbfWci1WcoCShBAkou11ekL3J4uJyy+VoG7R+6VwzpHj7ffdd9d7ndlLJ1+8G58AJKkICSS09A21dMOqWsWBE9tPvZDujN1ZcuVwWpvufq9oOrMQkoQQJKLj0B3Tf7Lcuvyi35spPnbf12QM9DuJy+Xqfxtff2g+uBCihBAsp8tY4DhTkSUOZLQJk5AWW+BJSZE1DmS0CZOQFlvgSUmRNQ5ktAmTkBBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQj6/19btBQEVfRuAAAAAElFTkSuQmCC)
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);
}