# The Dynamic Splitting Method with an application to portfolio credit risk

Supervisor: Dr. Zdravko Botev

Kevin Lam

Acknowledgements I would like to sincerely thank my supervisor Dr. Zdravko
Botev for his guidance, tutelage and intellectual discussions he has provided
throughout my Honours year. I am grateful for him taking time to share his
interesting research ideas that has perked my interests. Also, I must express
my gratitude to my colleagues for accompanying me on the same journey in 2015.
I would also like to thank the School of Mathematics and Statistics at UNSW
for the enriching experience over the last 5 years. Finally, I would like to
thank my family for their encouragement and support.

Abstract We consider the problem of accurately measuring the credit risk of a
portfolio consisting of loss exposures such as loans, bonds and other
financial assets. We are particularly interested in the probability of large
portfolio losses. We describe the popular models in the credit risk framework
including _factor models_ and _copula models_. To this end, we revisit the
most efficient probability estimation algorithms within current copula credit
risk literature, namely importance sampling. We illustrate the workings and
developments of these algorithms for large portfolio loss probability
estimation and quantile estimation. We then propose a modification to the
dynamic splitting method which allows application to the credit risk models
described. Our proposed algorithm for the unbiased estimation of rare-event
probabilities, exploits the quasi-monotonic property of functions to embed a
static simulation problem within a time-dependent Markov process. A study of
our proposed algorithm is then conducted through numerical experiments with
its performance benchmarked against current popular importance sampling
algorithms.

###### Keywords: Rare-event probability estimation; Monte Carlo methods;
Importance sampling; Splitting Method; Markov Processes

##  Chapter 1 Introduction

Simulation from intractable multidimensional distributions and estimates of
corresponding real-valued quantities are hallmark problems in Monte Carlo
methods; (see [26]). In this chapter we examine these two related problems
more closely as they frequently arise in Monte Carlo applications. The first
problem is to simulate from the conditional density in $\mathbb{R}^{d}$

$f^{*}(\mathbf{x})=\frac{1}{\ell}f(\mathbf{x})\mathbb{I}\\{S(\mathbf{x})\geq\gamma\\},\quad\mathbf{x}=(x_{1},\ldots,x_{d})^{\top},$
(1.1)

where we assume that: $S:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a real-valued
function, which we refer to as an _importance function_ ; $f$ is a density
function on $\mathbb{R}^{d}$ such that $X_{1},\ldots,X_{d}$ are independent;
$\gamma$ is a real parameter, and

$\ell=\mathbb{P}_{f}(S(\mathbf{X})\geq\gamma),\quad\mathbf{X}\sim f$ (1.2)

is a normalizing constant. $S$ can be interpreted as a measure of performance
of $\mathbf{X}$ with $\gamma$ being a loss threshold which, if exceeded,
triggers an event whose associated probability we wish to estimate. To
estimate the probability of this event, we must also solve a second problem,
and that is to estimate the normalizing constant $\ell$ accurately and
efficiently. $\ell$ can be interpreted the probability of default for a
financial institution or a probability of ruin for an insurance company if
$S(\mathbf{X})$ is the aggregate losses or claims incurred.

Note that, despite the specification of (1.1), the definition of the
conditional density function $f^{*}$ comes as a special case of any complex
high-dimensional density function
$\tilde{f}(\mathbf{x})=\breve{f}(\mathbf{x})/\mathcal{Z}$ on $\mathbb{R}^{d}$
with a known or unknown normalizing constant $\mathcal{Z}$. In other words,
without loss of generality, the conditional density (1.1) includes many models
arising in Bayesian inference and statistics, econometrics and finance.

In the unlikely case that $\ell$ is not a _rare-event probability_ (that is,
it is not too small; say, larger than $10^{-4}$), one can simulate from (1.1)
exactly with the acceptance-rejection algorithm. That is, we simulate
$\mathbf{X}\sim f$ until $S(\mathbf{X})\geq\gamma$ is satisfied. For the
simulated $\mathbf{X}$ for which this condition is satisfied, we accept this
as a sample from the conditional distribution (1.1). Unfortunately, more often
than not, $\ell$ is a rare-event probability and in such cases the only
practicable approach to simulate from (1.1) is approximate Markov Chain Monte
Carlo (MCMC) sampling as described in [9].

In this thesis, we utilize Monte Carlo methods to estimate these probabilities
of the form 1.2. In particular, we focus on a classical Monte Carlo technique,
Dynamic Splitting (DS), (see [23, 26]). In the original formulation of dynamic
splitting, the state space of a Markov process is decomposed into nested
subsets so that the rare event is represented as the intersection of
sequentially decreasing event subsets. Within each subset the sample paths of
the Markov process are split into multiple copies; the rationale behind this
is to promote and capture more occurrences of the rare event. As a result, the
probability of the rare event is calculated as the product of conditional
probabilities. Note that by splitting the Markov process into more copies, it
allows more accurate estimation of each conditional probability.

The remaining chapters of the thesis are organized as follows. In Chapter 2,
we cover essential background knowledge of Monte Carlo methods and describe
algorithms for the efficient estimation of $\ell$, namely _Importance
Sampling_ and _Cross-Entropy_. In Chapter 3, we introduce and survey the
framework of copula credit risk models in current literature with worked
examples of the algorithms presented in Chapter 2. Further, in Chapter 4, we
illustrate the workings of our proposed algorithm through numerical
experiments benchmark its performance against existing robust Monte Carlo
estimators. Moreover, critical analysis of advantages and limitations of the
proposed DS algorithm is given. Finally, in Chapter 5, we provide concluding
remarks and directions for future research. To preserve the flow of the
thesis, we delegate verifications and proofs of key results to the appendix.

##  Chapter 2 Background on Monte Carlo methods

The distribution of losses $F_{L}$ is the structure of interest in credit risk
modelling. There is often no closed form for $F_{L}$ making direct
calculations of $\ell(l)$ for a given $l$, or $l$ for a given $\ell$. To
estimate these measures we require a method to draw independent and
identically distributed (iid) samples of $L$ from $F_{L}$ and a method to
estimate probabilities and quantiles given an iid sample $L_{1},\dots,L_{N}$,
particularly those in the upper tail of $F_{L}$ such as the Value-at-Risk
(VaR) and Conditional Value-at-Risk (CVaR). The purpose of this chapter is to
describe the estimators of probabilities and quantiles using Monte Carlo
methods, namely _Importance Sampling_.

###  2.1 Efficiency

In this section, we describe the idea of efficiency as we require a criteria
to benchmark performance of estimators illustrated in this paper. For a rare-
event estimator, this is often measured by its _relative error_ (RE)[26]. This
is the normalized standard deviation of the estimator. Suppose we have an
unbiased estimator $\hat{\ell}$ of the rare-event probability defined as

$\begin{split}\ell(\gamma)&=\mathbb{P}(S(\mathbf{X})\geq\gamma)\\\
&=\mathbb{E}[\mathbb{I}\\{S(\mathbf{X})\geq\gamma\\}]\end{split}$

where $f$ is a probability density function (pdf),$S$ is a real-valued
function, $\mathbf{X}$ is a random vector and $\gamma$ is a threshold
parameter. The event of interest is $\\{S(\mathbf{X})\geq\gamma\\}$ occurring
under $f$. RE is defined as

$RE(\hat{\ell})=\frac{\sqrt{\mathrm{Var}(\hat{\ell})}}{\ell\sqrt{N}}$

where $N$ is the number of iid estimates of $\ell$.

###  2.2 Crude Monte Carlo

_Crude Monte Carlo_ (CMC) [32] approximates the cumulative distribution
function (cdf) $F_{L}$ with the empirical distribution function (edf)
$\hat{F}_{L}$. Let $L_{1},\dots,L_{N}$, be an iid sample of size $N$, then
$\hat{F}_{L}$ defined as follows

$\hat{\ell}(\gamma)=\frac{1}{N}\sum_{k=1}^{N}\mathbb{I}(L_{k}\geq\gamma).$

We define $\hat{v}_{\alpha}$ as the CMC estimate of the $\alpha$-VaR or
equivalently the $\alpha$-quantile. By using the above definition, CMC sets
$\hat{v}_{\alpha}$ as the solution to the problem

$\hat{v}_{\alpha}=\inf\\{\gamma:\hat{\ell}(\gamma)\leq 1-\alpha\\}.$

As a result, we can also obtain $\hat{v}_{\alpha}$ by sorting the
$\\{L_{k}\\}$ in ascending order and taking the $\lceil\alpha N\rceil$-th
largest value. Once $\hat{v}_{\alpha}$ is known, we can calculate the CMC
estimator for the CVaR as follows

$\hat{c}_{\alpha}=\frac{1}{N(1-\alpha)}\sum_{k=1}^{N}L_{k}\mathbb{I}(L_{k}\geq\hat{v}_{\alpha}).$

It is well-known that the VaR of a portfolio is not additive, that is, it is
not the sum of the VaR of sub-portfolios or individual losses as it is not the
sum of independent random variables. Nevertheless,[21] provides some insight
into the asymptotic distribution of these estimators through the following
central limit theorems.

######  Theorem 2.2.1 (Central Limit Theorems for CMC estimators).

If $L$ is positive and continuously differentiable density $f_{L}$ around
$v_{\alpha}$ and $\mathbb{E}[L^{2}]<\infty$, then as $N\rightarrow\infty$

$\sqrt{N}(\hat{v}_{\alpha}-v_{\alpha})\xrightarrow{d}\mathsf{N}\left(0,\frac{\alpha(1-\alpha)}{f_{L}(v_{\alpha})^{2}}\right),$
$\sqrt{N}(\hat{c}_{\alpha}-c_{\alpha})\xrightarrow{d}\mathsf{N}\left(0,\frac{\mathrm{Var}(L\mathbb{I}(L>v_{\alpha}))}{(1-\alpha)^{2}}\right).$

###  2.3 Importance Sampling

Our aim is to efficiently estimate $\ell(\gamma)$ for a given $\gamma$, or
$\gamma$ for a given $\ell$. CMC estimators heavily rely on the computational
power of a large sample size to achieve accuracy. However, if $\ell$ is a
rare-event probability say, smaller than $10^{-4}$, then generating large
samples is costly to simulate large values of $\gamma$. Much of the work on
simulation methods has been done on variance reduction algorithms to improve
the efficiency of the CMC estimator. We now present an algorithm that is
particularly well-suited to rare-event problems, _Importance Sampling_ (IS).
IS has been shown to accurately and efficiently sample from the upper tail of
a general loss distribution, and hence reduce the variance of $\hat{\ell}$
[18].

Suppose that $L$ can be simulated under another density $g_{L}$, which we
denote as the IS density. Define,

$W(\gamma)=\frac{f_{L}(\gamma)}{g_{L}(\gamma)}$

as the likelihood ratio obtained by the change in probability measures. Note
that

$\displaystyle\mathbb{P}(L>\gamma)$
$\displaystyle=\mathbb{E}_{f}[\mathbb{I}\\{L>\gamma\\}]$
$\displaystyle=\mathbb{E}_{g}[W(L)\mathbb{I}\\{L>\gamma\\}]$

where $\mathbb{E}_{f}$ and $\mathbb{E}_{g}$ denote expectations with respect
to pdf $f_{L}$ and $g_{L}$ respectively. As we are interested in the upper
tail of loss distribution, we can thus estimate tail probabilities as follows

$\hat{\ell}=\frac{1}{N}\sum_{k=1}^{N}W(L_{k})\mathbb{I}\\{L_{k}>\gamma\\}.$

[19] suggests the IS approach to quantile estimation by defining the IS cdf as

$\begin{split}\hat{F}_{L}^{IS}(\gamma)&=1-\hat{\ell}\\\
&=1-\frac{1}{N}\sum_{k=1}^{N}W(L_{k})\mathbb{I}\\{L_{k}>\gamma\\}.\end{split}$

With this definition, the IS estimators for the VaR and CVaR are given by

$\hat{v}_{\alpha}^{IS}=\inf\\{\gamma\hat{F}_{L}^{IS}\geq\alpha\\},$
$\hat{c}_{\alpha}^{IS}=\frac{1}{N(1-\alpha)}\sum_{k=1}^{N}W(L_{k})L_{k}\mathbb{I}\\{L_{k}\geq\hat{v}_{\alpha}^{IS}\\}.$

Similar to the CMC estimators, [22] show that the IS estimators asymptotically
follow Normal distributions and give the following central limit theorems.

######  Theorem 2.3.1 (Central Limit Theorems for IS estimators).

If $L$ is positive and continuously differentiable density $f_{L}$ around
$v_{\alpha}$ and there exists $\epsilon>0$ and $p>2$ such that $W(\gamma)$ is
bounded for all $\gamma\in(v_{\alpha}-\epsilon,v_{\alpha}+\epsilon)$ and
$\mathbb{E}_{g}[\mathbb{I}(L\geq v_{\alpha}-\epsilon)(W(L))^{p}]<\infty$, then
as $N\rightarrow\infty$

$\sqrt{N}(\hat{v}_{\alpha}^{IS}-v_{\alpha})\xrightarrow{d}\mathsf{N}\left(0,\frac{\mathrm{Var}_{g}(W(L))\mathbb{I}(L\geq
v_{\alpha}))}{f_{L}(v_{\alpha})^{2}}\right),$
$\sqrt{N}(\hat{c}_{\alpha}^{IS}-c_{\alpha})\xrightarrow{d}\mathsf{N}\left(0,\frac{\mathrm{Var}_{g}(W(L))L\mathbb{I}(L>v_{\alpha}))}{(1-\alpha)^{2}}\right).$

where $\mathrm{Var}_{g}$ denotes the variance under the pdf $g_{L}$. Note that
$g_{L}$ is not known so the above theorems only describe attractive properties
of a good choice of $g_{L}$. A good choice for the IS pdf depends on the
distribution of $\boldsymbol{X}$, properties of the set
$\\{S(\mathbf{X})\geq\gamma\\}$ and more importantly, the tail behaviour of
$S(\mathbf{x})$. A light-tailed random variable $X$ is defined as one which
has a finite moment generating function (mgf), that is $\mathbb{E}[e^{\theta
X}]<\infty$ for $\theta>0$. [8] suggests that a good IS pdf $g_{L}$ in a
light-tailed setting is an exponentially twisted pdf derived from $f_{L}$,
defined as

$g_{L}(\gamma)=\frac{\exp(\theta\gamma)f_{L}(\gamma)}{\mathbb{E}[e^{\theta
L}]}.$

Here $\mathbb{E}[e^{\theta L}]$ acts as the normalizing constant to ensure
$g_{L}$ is a density. The likelihood ratio of an exponentially twisted pdf is
thus given by

$W(\gamma)=E[e^{\theta L}]\exp(-\theta\gamma).$ (2.1)

Attractive properties of likelihood ratios of this form are described in [10,
12].

####  2.3.1 Adaptive Importance Sampling

Adaptive importance sampling methods aim to avoid theoretical complications
and computational issues in rare-event probability estimation by deriving
parameters for an optimal IS density by using sub-samples of sample data. As
described previously, there is often no closed form for $f_{L}$ so we
represent the portfolio loss in the form $L=S(\mathbf{X})$ and seek to
parameterize a prespecified IS density. As this is discussed within a credit
risk framework, our interest is to simulate from a density conditional on the
event $\\{S(\mathbf{X})\geq\gamma\\}$ where $\gamma$ is often chosen to be
large loss threshold. If an initial sample
$\mathbf{X}_{1},\ldots,\mathbf{X}_{M}$ can be generated directly from the
optimal IS density, that is, the zero-variance density
$g^{*}(\mathbf{x})=f(\mathbf{x}|S(\mathbf{x})\geq\gamma)$ then the parameters
can be computed to approximate $g^{*}$. A popular and versatile adaptive
importance sampling method is _Cross-Entropy_ (CE)[28, 29]. The goal of CE is
to specify a density $g$ ‘close’ to $g^{*}$ so that both would behave
similarly and give reasonably accurate IS estimators. We consider the family
of distributions $\mathcal{G}=\\{g(\boldsymbol{x};\boldsymbol{v})\\}$ where
$\boldsymbol{v}$ denotes a vector of parameters. A convenient measure of the
difference between two densities $g_{1}$ and $g_{2}$ is the _Kullback-Leibler
divergence_ also known as the _cross-entropy distance_ :

$\mathcal{D}(g_{1},g_{2})=\int
g_{1}(\mathbf{x})\log\frac{g_{1}(\mathbf{x})}{g_{2}(\mathbf{x})}\mathrm{d}\mathbf{x}.$

Every density in $\mathcal{G}$ can be represented as $g(\cdot;\boldsymbol{v})$
for some $\boldsymbol{v}$, so we obtain the optimal IS density by solving the
following problem:

$\boldsymbol{v}^{*}_{CE}=\operatorname*{\arg\\!\min}_{\boldsymbol{v}}\mathcal{D}(g^{*},g(\cdot;\boldsymbol{v})).$

This is equivalent [28] to solving

$\boldsymbol{v}^{*}_{CE}=\operatorname*{\arg\\!\max}_{\boldsymbol{v}}\mathbb{E}[f(\mathbf{X})\mathbb{I}\\{S(\mathbf{X})\geq\gamma\\}\log
g(\mathbf{X};\boldsymbol{v})].$

This problem does not typically have an explicit solution. Instead we can
estimate $\boldsymbol{v}^{*}_{CE}$ by solving the following problem:

$\hat{\boldsymbol{v}}^{*}_{CE}=\operatorname*{\arg\\!\max}_{\boldsymbol{v}}\frac{1}{N}\sum_{i=1}^{N}\mathbb{I}\\{S(\mathbf{X}_{i})\geq\gamma\\}\log
g(\mathbf{X}_{i};\boldsymbol{v})]$

where $\mathbf{X}_{1},\ldots.\mathbf{X}_{N}$ are simulated from $f$. This is
further simplified if we are able to draw approximately from $g^{*}$, reducing
the problem to

$\hat{\boldsymbol{v}}^{*}_{CE}=\operatorname*{\arg\\!\max}_{\boldsymbol{v}}\frac{1}{N}\sum_{i=1}^{N}\log
g(\mathbf{X}_{i};\boldsymbol{v})]$

where $\mathbf{X}_{1},\ldots.\mathbf{X}_{N}$ are drawn approximately from
$g^{*}$. Suppose we want to sample approximately from the zero-variance
density $g^{*}(\mathbf{x})=f(\mathbf{x}|S(\mathbf{x})>v_{\alpha})$. Given a
generated sample of portfolio losses $L_{1},\ldots,L_{N}$ with corresponding
vectors $\mathbf{X}_{L_{1}},\ldots,\mathbf{X}_{L_{N}}$ from $f$, we can order
the portfolio losses in ascending order as $L_{(1)}\leq\cdots\leq L_{(N)}$ and
choose $L_{(\lceil\alpha N\rceil)},\ldots,L_{(N)}$ with corresponding vectors
$\mathbf{X}_{L_{(\lceil\alpha N\rceil)}},\ldots,\mathbf{X}_{L_{(N)}}$ as an
approximate sample from $g^{*}$. With this approximate sample, we can use
standard maximum likelihood estimation to parameterize a prespecified density
$g$ which approximates $g^{*}$. The approximate sampling algorithm can be
summarized as follows.

Algorithm 1 : Sampling approximately from $g^{*}$

0: distribution of $\mathbf{X}$ $f$; importance function $S$; sample size $N$;
loss threshold $x$

$t\leftarrow 0$

repeat

$t\leftarrow t+1$

Simulate $\mathbf{X}_{t}$ from $f$

$B_{k}\leftarrow\mathbb{I}\\{X_{k}>x_{k}\\}$ for $k=1,\ldots,d$

Generate the portfolio loss $L_{t}\leftarrow S(\mathbf{X}_{t})$

until $t=N$

Sort the losses $\boldsymbol{L}\leftarrow(L_{[1]},\ldots,L_{[N]})^{\top}$

return $L_{(\lceil\alpha N\rceil)},\ldots,L_{(N)}$ and
$\mathbf{X}_{L_{(\lceil\alpha N\rceil)}},\ldots,\mathbf{X}_{L_{(N)}}$

##  Chapter 3 Introduction to Copula Credit Risk Models

###  3.1 Summary of the Model

The model can be summarized as follows. Let

$L=c_{1}B_{1}+\cdots+c_{d}B_{d}=\boldsymbol{c}^{\top}\boldsymbol{B}$

be the total loss incurred by a portfolio of $d$ obligors, where $c_{k}$ is
the loss incurred from the $k$-th obligor and $B_{i}\sim\mathsf{B}(P_{k})$ is
a Bernoulli random variable indicating whether the $k$-th obligor has
defaulted. The distribution of the column vector $\boldsymbol{B}$ is
implicitly defined under factor models and copula models which are described
later in this chapter. For now, we note that dependence of default events is
captured in the model through the default probabilities $\\{P_{k}\\}$ which
contain a set of common factors, say $\Psi$, that affect all obligors.
Conditional on $\Psi$, the problem is simplified to modelling the sum of
independent Bernoulli variables $\\{B_{k}\\}$ scaled by the losses
$\\{c_{k}\\}$. These models have popular applications in finance [5, 13, 14],
particularly in the valuation of credit risk such as the estimation of Value-
at-Risk for a given confidence level $\alpha$. Popular values for $\alpha$ are
$0.95$, $0.99$ and $0.995$ since the loss values of interest lie in the upper
tail of the loss distribution. The purpose of this chapter is to describe the
most popular credit risk models, _factor models_ and _copula models_ , with
applications of the algorithms presented in Chapter 2.

###  3.2 Factor Models

In factor models, $\boldsymbol{B}$ are defined as

$B_{k}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbb{I}\\{X_{k}>x_{k}\\},\quad
k=1,\ldots,d$

where $\\{x_{k}\\}$ are given fixed thresholds and $\boldsymbol{X}$ has a
continuous joint density $f(\boldsymbol{x})$. For example,
$\boldsymbol{X}\sim\mathsf{N}(\boldsymbol{\mu},\mathrm{\Sigma})$ for some mean
vector $\boldsymbol{\mu}$ and covariance matrix $\mathrm{\Sigma}$. Note that
$\mathrm{\Sigma}$ can be singular.

####  3.2.1 Gaussian factor model

A popular Gaussian factor model for $\boldsymbol{X}$ is

$X_{k}=a_{k1}Z_{1}+\cdots+a_{km}Z_{m}+b_{k}\epsilon_{k},\quad k=1,\ldots,d$
(3.1)

where
$Z_{1},\cdots,Z_{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathsf{N}(0,1)$
are the so called systematic risk factors which affect all obligors;
$a_{k1},\cdots,a_{kd}$ are default factor loadings for the $k$-th obligor with

$a_{k1}^{2}+\cdots+a_{km}^{2}\leq 1,$
$b_{k}=\sqrt{1-(a_{k1}^{2}+\cdots+a_{km}^{2})}.$

and $\epsilon_{k}\sim\mathsf{N}(0,1)$ is risk specific to the $k$-th obligor.
We thus have the marginal distribution $X_{k}\sim\mathsf{N}(0,1)$. We can
write

$\boldsymbol{X}=\mathrm{A}\boldsymbol{Z}+\mathrm{diag}(\boldsymbol{b})\boldsymbol{\epsilon},$

where $\mathrm{A}$ is an $d\times m$ matrix.

Here, probabilities are conditionally independent on $\boldsymbol{Z}$, that is
$\Psi=\boldsymbol{Z}$, the default probability of the $k$-th obligor
$P_{k}(\boldsymbol{Z})$ is

$\begin{split}P_{k}(\boldsymbol{Z})&=\mathbb{P}(X_{k}>x_{k}|\boldsymbol{Z}=\boldsymbol{z})\\\
&=\mathbb{P}\left(\epsilon_{k}>\frac{x_{k}-(a_{k1}z_{1}+\cdots+a_{km}z_{m})}{b_{k}}\right)\\\
&=\Phi\left(\frac{(a_{k1}z_{1}+\cdots+a_{km}z_{m})-x_{k}}{b_{k}}\right).\end{split}$
(3.2)

#####  3.2.1.1 Crude Monte Carlo

Suppose that the factor loadings matrix $\mathrm{A}$ and default thresholds
$\boldsymbol{x}$ are known. For a chosen sample size $N$ and confidence level
$\alpha$, the CMC algorithm is as follows.

Algorithm 2 : Generating $L$ under a Gaussian factor model using CMC

0: factor loadings $\mathrm{A}$; cost vector $\boldsymbol{c}$; default
thresholds $\boldsymbol{x}$; sample size $N$; confidence level $\alpha$; loss
threshold $\gamma$

$t\leftarrow 0$

repeat

$t\leftarrow t+1$

Simulate $\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{0},\mathrm{I}_{m})$

Simulate $\boldsymbol{\epsilon}\sim\mathsf{N}(\boldsymbol{0},\mathrm{I}_{d})$

$\boldsymbol{X}\leftarrow\mathrm{A}\boldsymbol{Z}+\mathrm{diag}(\boldsymbol{b})\boldsymbol{\epsilon}$

$B_{k}\leftarrow\mathbb{I}\\{X_{k}>x_{k}\\}$ for $k=1,\ldots,d$

$L_{t}\leftarrow\boldsymbol{c}^{\top}\boldsymbol{B}$

until $t=N$

Sort the losses $\boldsymbol{L}\leftarrow(L_{[1]},\ldots,L_{[N]})^{\top}$

return
$\hat{\ell}(\gamma)\leftarrow\frac{1}{N}\sum_{k=1}^{N}\mathbb{I}(L_{k}>\gamma)$
or $\hat{v}_{\alpha}\leftarrow L_{\lceil\alpha\times N\rceil}$

#####  3.2.1.2 Importance Sampling

Suppose that the conditional default probabilities
$\boldsymbol{P}(\boldsymbol{Z})$ are known, [18] describes an IS algorithm
which changes the default indicators to
$B_{k}\sim\mathsf{B}(P_{k}(\boldsymbol{Z}))$ before applying exponential
twisting to the conditional probabilities $P_{k}(\boldsymbol{Z})$ as follows

$P_{k,\theta}(\boldsymbol{Z})=\frac{P_{k}(\boldsymbol{Z})e^{\theta
c_{k}}}{1+P_{k}(\boldsymbol{Z})(e^{\theta c_{k}}-1)}.$ (3.3)

By applying this change in probability measure, conditional default
probabilities are increased if $\theta>0$. Exponential twisting are well known
for the significant variance reduction in range of contexts [1]. It is
difficult to apply exponential twisting directly on $f_{L}$ so it is applied
on the default indicators $B_{k}$ instead and this can be shown to give
equivalent results [18]. Conditional on $\Phi=\boldsymbol{Z}$, $L$ becomes a
sum of independent scaled Bernoulli random variables $c_{k}B_{k}$. Now each
$c_{k}B_{k}$ has mgf

$\mathbb{E}[e^{\theta
c_{k}B_{k}}]=(1-P_{k}(\boldsymbol{Z}))+P_{k}(\boldsymbol{Z})e^{\theta
c_{k}}<\infty$

so $\\{c_{k}B_{k}\\}$ are light-tailed with a good IS pdf being an exponential
twisted density obtained by applying (3.3). It can easily be verified that
this leads to the likelihood ratio for $L$ given by (2.1). The algorithm,
which [18] calls the one-step algorithm, can be summarized as follows. Given
$\boldsymbol{Z}$ and $\\{x_{k}\\}$, we calculate the conditional probabilities
$P_{k}(\boldsymbol{Z})$ and apply exponential twisting to each probability.
This generates an exponential twisted density for each random variable
$c_{k}B_{k}$; the product of these densities forms the exponentially twisted
density $g_{L}$. The likelihood ratio is then used to estimate upper tail
probabilities from which quantiles in the tail can then be estimated. The one-
step IS algorithm is as follows.

Algorithm 3 : Generating $L$ using one-step IS algorithm

0: factor loadings $\mathrm{A}$; cost vector $\boldsymbol{c}$; default
thresholds $\boldsymbol{x}$; sample size $N$; confidence level $\alpha$; loss
threshold $\gamma$

$t\leftarrow 0$

repeat

$t\leftarrow t+1$

Simulate $\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{0},\mathrm{I}_{m})$

Calculate vector $\boldsymbol{P}(\boldsymbol{Z})$ with
$P_{k}(\boldsymbol{Z})\leftarrow\Phi\left(\frac{\sum_{i}\mathrm{A}_{ki}Z_{i}-x_{k}}{b_{k}}\right)$

if $\boldsymbol{c}^{\top}\boldsymbol{P}(\boldsymbol{Z})\geq\gamma$ then

$\theta\leftarrow 0$

else

$\theta$ is assigned to be the solution to
$\gamma=\sum_{k=1}^{m}\frac{P_{k}(\boldsymbol{Z})c_{k}\exp(\theta
c_{k})}{1+P_{k}(\boldsymbol{Z})(\exp(\theta c_{k})-1)}$

for $k=1,\ldots,d$ do

$P_{k,\theta}(\boldsymbol{Z})\leftarrow\frac{P_{k}(\boldsymbol{Z})\exp(\theta
c_{k})}{1+P_{k}(\boldsymbol{Z})(\exp(\theta c_{k})-1)}$

Simulate $B_{k}\sim\mathsf{B}(P_{k,\theta}(\boldsymbol{Z}))$, independently

Set $L_{t}\leftarrow\boldsymbol{c}^{\top}\boldsymbol{B}$ and
$W(L_{t})\leftarrow\exp\left(-\theta
L_{t}+\sum_{k=1}^{m}\log\left(1+P_{k}(\boldsymbol{Z})(\exp(\theta
c_{k})-1)\right)\right)$

until $t=N$

return
$\hat{\ell}(\gamma)\leftarrow\frac{1}{N}\sum_{t=1}^{N}\mathbb{I}\\{L_{t}>\gamma\\}W(L_{t})$
or
$\hat{v}_{\alpha}^{IS}\leftarrow\inf\\{\gamma:\hat{F}_{L}^{IS}(\gamma)\geq\alpha\\}$

The efficiency of the algorithm depends on the level of dependence between
obligors. [18] states that when dependence is weak, such as the case of the
Gaussian factor model, increasing conditional default probabilities by
exponential twisting already reduces variance in the estimators of tail
probabilities and quantiles. The algorithm aims to simulate large values for
$L$ centred around a carefully chosen loss threshold $\gamma$ that lies in the
tail of $f_{L}$. This is done by solving for the unique value of $\theta$ as
described in the algorithm and applying the exponential twisting with this
value of $\theta$. An attractive property of exponential twisted pdf $g_{L}$,
that can be easily verified, is that for each simulated $L$ we have

$\mathbb{E}_{g}[L|\boldsymbol{Z}]\geq\gamma,$

which comes a result of solving for $\theta$ such that

$\gamma=\sum_{k=1}^{m}\frac{P_{k}(\boldsymbol{Z})c_{k}\exp(\theta
c_{k})}{1+P_{k}(\boldsymbol{Z})(\exp(\theta c_{k})-1)}.$

Hence, the algorithm can sample from the tail of $f_{L}$ if we chose our loss
threshold $\gamma$ to be a quantile in the tail, such as $\hat{v}_{0.95}$,
which can be initially estimated with CMC. With this sampling, small portfolio
losses are now rare events while large portfolio losses are now frequent.

The vector of likelihood ratios, $\boldsymbol{W}$, can be used to estimate
upper tail probabilities. After sorting the simulated losses in ascending
order, the IS estimator $\hat{v}_{\alpha}^{IS}$ can be computed as
$\mathrm{VaR}=L_{(j)}$ where $j$ is the solution to

$\min\limits_{j}\left(\frac{1}{N}\sum_{k=j}^{N}W(L_{(k)})\leq 1-\alpha\right)$

[18] states a further extension to form a two-step algorithm. This is
motivated by the variance decomposition of the estimator $\hat{\ell}$

$\mathrm{Var}(\hat{\ell})=\mathbb{E}\left(\mathrm{Var}(\hat{\ell|\boldsymbol{P}})\right)+\mathrm{Var}\left(\mathbb{E}(\hat{\ell}|\boldsymbol{P})\right).$

The one-step IS algorithm minimizes the variability of $\hat{\ell}$, that is,
it minimizes $\mathrm{Var}(\hat{\ell|\boldsymbol{P}})$. The two-step IS
algorithm aims to minimize
$\mathrm{Var}\left(\mathbb{E}(\hat{\ell}|\boldsymbol{P})\right)$. This is
equivalent to minimizing the variance of the CMC estimator $\hat{q}$ [8] of

$q=\mathbb{P}(L>\gamma|\boldsymbol{P}(\boldsymbol{Z})).$

The corresponding zero-variance density $g^{*}$ [26] is given by

$g^{*}_{\boldsymbol{Z}}(\boldsymbol{z})\propto\mathbb{P}(L>\gamma|\boldsymbol{P}(\boldsymbol{Z}))f_{\boldsymbol{Z}}(\boldsymbol{z}).$

However, we must note that the normalizing constant $\ell$ is the same
constant we wish to estimate so this is not a practical IS density.
Nevertheless, this provides a direction in the searching for a good IS
density. A common approach to by applying the one-step algorithm after a
change of measure for $\boldsymbol{Z}$. In particular, we change the mean of
$\boldsymbol{Z}$. The underlying rationale of this approach is to generate
more defaults by shifting the mean of the factors $\boldsymbol{Z}$ by
increasing each of its components. This leads to high values for the default
factor loadings $\\{X_{k}\\}$ which are more likely to exceed the thresholds
$\\{x_{k}\\}$. The challenge in the two-step algorithm is to describe a
suitable IS distribution $g_{Z}$. [18] and [15] propose using a Normal
distribution $\mathsf{N}(\boldsymbol{\mu},\mathrm{I}_{m})$ with the same mode
as optimal pdf $g_{Z}^{*}$. The mode $\boldsymbol{\mu}^{*}$ is also the mean
of the Normal distribution and provided that we have loss threshold $x$, is
given as the solution to the following problem

$\boldsymbol{\mu}^{*}=\operatorname*{\arg\\!\max}\limits_{\boldsymbol{z}}\mathbb{P}\left(L>\gamma|\boldsymbol{Z}=\boldsymbol{z}\right)\exp{(-\frac{1}{2}\boldsymbol{z}^{\top}\boldsymbol{z})}.$
(3.4)

Hence the two-step algorithm applies a change in distribution to the factors
$\boldsymbol{Z}$ and simulates
$\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{\mu}^{*},\mathrm{I}_{m})$. The
difficulty now lies in solving (3.4). [18] states several approximations to
simplify this problem. We focus on the _constant approximation_ and the _tail
bound approximation_ as it provides the most convenient way to combine IS
applied on the probabilities $P_{k}(\boldsymbol{Z})$ with IS applied on the
factors $\boldsymbol{Z}$. For details on other approximations used in solving
(3.4), interested readers may refer to [18].

The _constant approximation_ involves replacing $L$ with
$\mathbb{E}[L|\boldsymbol{Z}=\boldsymbol{z}]$ and
$\mathbb{P}(L>\gamma|\boldsymbol{Z}=\boldsymbol{z})$ with
$\mathbb{I}(\mathbb{E}[L|\boldsymbol{Z}=\boldsymbol{z}]>\gamma)$. This
approximation replaces
$\mathbb{P}\left(L>\gamma|\boldsymbol{Z}=\boldsymbol{z}\right)$ with a
constant and so (3.4) becomes

$\operatorname*{\arg\\!\min}\limits_{\boldsymbol{z}}\\{\boldsymbol{z}^{\top}\boldsymbol{z}:\mathbb{E}[L|\boldsymbol{Z}=\boldsymbol{z}]>\gamma\\}.$
(3.5)

The _tail bound approximation_ is an approach which aims to approximate
$\mathbb{P}\left(L>\gamma|\boldsymbol{Z}=\boldsymbol{z}\right)$ by its upper
bound. It then proceeds by maximising this upper bound which in turn, could
maximise the probability
$\mathbb{P}\left(L>\gamma|\boldsymbol{Z}=\boldsymbol{z}\right)$. Using this
approximation, (3.4) becomes

$\operatorname*{\arg\\!\max}\limits_{\boldsymbol{z}}\left\\{\sum_{k=1}^{m}\log[1+P_{k}(\boldsymbol{Z})(\exp(\theta
c_{k})-1)]-\theta\gamma-\frac{1}{2}\boldsymbol{z}^{\top}\boldsymbol{z}\right\\}.$
(3.6)

The two-step algorithm for dependent obligors is as follows.

Algorithm 4 : Generating $L$ using Glasserman and Li’s two-step algorithm

0: factor loadings $\mathrm{A}$; cost vector $\boldsymbol{c}$; default
thresholds $\boldsymbol{x}$; sample size $N$; confidence level $\alpha$; loss
threshold $\gamma$; shifted mean vector $\boldsymbol{\mu}^{*}$

$t\leftarrow 0$

repeat

$t\leftarrow t+1$

Simulate $\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{\mu},\mathrm{I}_{m})$

Calculate vector $\boldsymbol{P}(\boldsymbol{Z})$ with
$P_{k}(\boldsymbol{Z})\leftarrow\Phi\left(\frac{\sum_{j}\mathrm{A}_{k,j}Z_{j}-x_{k}}{b_{k}}\right)$

if $\boldsymbol{c}^{\top}\boldsymbol{P}(\boldsymbol{Z})\geq\gamma$ then

$\theta\leftarrow 0$

else

$\theta$ is assigned to be the solution to
$\gamma=\sum_{k=1}^{m}\frac{P_{k}(\boldsymbol{Z})c_{k}\exp(\theta
c_{k})}{1+P_{k}(\boldsymbol{Z})(\exp(\theta c_{k})-1)}$

for $k=1,\ldots,d$ do

$P_{k,\theta}(\boldsymbol{Z})\leftarrow\frac{P_{k}(\boldsymbol{Z})\exp(\theta
c_{k})}{1+P_{k}(\boldsymbol{Z})(\exp(\theta c_{k})-1)}$

$B_{k}\leftarrow\mathsf{Ber}(P_{k,\theta})$, independently

Set $L_{t}\leftarrow\boldsymbol{c}^{\top}\boldsymbol{B}$ and
$W(L_{t})\leftarrow\exp\left(-\theta
L+\sum_{k=1}^{m}\log\left(1+P_{k}(\boldsymbol{Z})(\exp(\theta
c_{k})-1)\right)\right)\exp\left(-\boldsymbol{\mu}^{*\top}\boldsymbol{Z}+\frac{\boldsymbol{\mu}^{*\top}\boldsymbol{\mu}^{*}}{2}\right)$

until $t=N$

return
$\hat{\ell}(\gamma)\leftarrow\frac{1}{N}\sum_{t=1}^{N}\mathbb{I}(L_{t}>\gamma)W(L_{t})$
or
$\hat{v}_{\alpha}^{IS}\leftarrow\inf\\{\gamma:\hat{F}_{L}^{IS}(\gamma)\geq\alpha\\}$

####  3.2.2 $t$ Factor Model

The $t$ factor model [11] differs from the Gaussian factor model as the
factors now have a multivariate $t$ distribution rather than a multivariate
Normal distribution. Following the desired properties and notation from (3.1),
$\boldsymbol{X}$ generated from a $t$ factor model usually has the
representation

$X_{k}=\sqrt{\frac{r}{V}}\left(a_{k1}Z_{1}+\cdots+a_{km}Z_{m}+b_{k}\epsilon_{k}\right),\quad
k=1,\ldots,d$

where
$Z_{1},\cdots,Z_{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathsf{N}(0,1)$
and $V\sim\mathsf{\chi}^{2}(r)$. Note $X_{k}$ is in the form

$X_{k}=\frac{Z}{\sqrt{V/r}},$

where $Z\sim\mathsf{N}(0,1)$ and $V\sim\mathsf{\chi}^{2}(r)$. This implies
that we have the marginal distribution $X_{k}\sim\mathsf{T}(r)$, which we
denote as a $t$ distribution with $r$ degrees of freedom.

Here, probabilities are conditionally independent on
$\boldsymbol{Z}=\boldsymbol{z}$ and $V=v$, that is
$\Psi=\\{\boldsymbol{Z},V\\}$, the default probability of the $k$-th obligor
$P_{k}(\boldsymbol{Z},V)$ is

$\begin{split}P_{k}(\boldsymbol{Z},V)&=\mathbb{P}(X_{k}>x_{k}|\boldsymbol{Z}=\boldsymbol{z},V=v)\\\
&=\mathbb{P}\left(\epsilon_{k}>\frac{\sqrt{\frac{v}{r}}x_{k}-(a_{k1}z_{1}+\cdots+a_{km}z_{m})}{b_{k}}\right)\\\
&=\Phi\left(\frac{(a_{k1}z_{1}+\cdots+a_{km}z_{m})-\sqrt{\frac{v}{r}}x_{k}}{b_{k}}\right).\end{split}$
(3.7)

#####  3.2.2.1 Crude Monte Carlo

Suppose that the factor loadings matrix $\mathrm{A}$ and default thresholds
$\boldsymbol{x}$ are known. For a chosen sample size $N$, confidence level
$\alpha$ and degrees of freedom $r$, the CMC algorithm is as follows

Algorithm 5 : Generating $L$ under a $t$ factor model using CMC

0: factor loadings $\mathrm{A}$; cost vector $\boldsymbol{c}$; default
thresholds $\boldsymbol{x}$; sample size $N$; confidence level $\alpha$; loss
threshold $\gamma$; degrees of freedom $r$

$k\leftarrow 0$

repeat

$k\leftarrow k+1$

Simulate $V\sim\mathsf{\chi}^{2}(r)$

Simulate $\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{0},\mathrm{I}_{m})$

Simulate $\boldsymbol{\epsilon}\sim\mathsf{N}(\boldsymbol{0},\mathrm{I}_{d})$

$\boldsymbol{X}\leftarrow\sqrt{\frac{r}{V}}\left(\mathrm{A}\boldsymbol{Z}+\mathrm{diag}(\boldsymbol{b})\boldsymbol{\epsilon}\right)$

$B_{i}\leftarrow\mathbb{I}\\{X_{i}>x_{i}\\}$ for $i=1,\ldots,d$

$L_{k}\leftarrow\boldsymbol{c}^{\top}\boldsymbol{B}$

until $k=N$

Sort the losses $\boldsymbol{L}\leftarrow(L_{[1]},\ldots,L_{[N]})^{\top}$

return $\hat{\ell}(\gamma)=\frac{1}{N}\sum_{k=1}^{N}\mathbb{I}(L_{k}>\gamma)$
or $\hat{v}_{\alpha}\leftarrow L_{\lceil\alpha\times N\rceil}$

#####  3.2.2.2 Importance Sampling

From the well-founded IS framework under the Gaussian factor model, much of
the literature consider exponential twisting to the $t$ factors by first
conditioning on $V$. Under $t$ factor models we note that the factors
$\boldsymbol{Z}$ contribute little to the occurrence of defaults as opposed to
the case of Gaussian factor models. Rather, it is the value of $V$ that plays
a much bigger role as a common shock factor. [24] applies exponential twisting
to $V$, with the twisting parameter $\theta_{V}$ found as the solution to a
linearly constrained optimization problem. The approach makes use of the
observation that, conditional on $V$, one can apply the same IS algorithm as
in [18] with modified thresholds for each obligor. By asymptotically optimal
results, it forms an approximate zero-variance IS pdf for $V$ and implicitly
applies the _constant approximation_ through the constrained sets specified in
the optimization problem when solving for the twisting parameter $\theta_{V}$.
The algorithm is however, computationally expensive as it requires solving the
optimization problem multiple times per sample. [24] suggests combining
stratified sampling with IS to reduce this cost.

For a general single factor model where $V$ need not follow a Chi-squared or
Gamma distribution, [4] presents two IS algorithms which apply exponential
twisting on $V$ and $W=\frac{1}{V}$ with the twisting parameter $\theta$ found
by solving optimization problems on the uniform upper bound of the likelihood
ratio estimator in a similar manner to the _tail bound approximation_. Another
recent advancement applies a shift in factors as described in the two-step IS
algorithm and adjusting the degrees of freedom of $V$ to maintain independence
[30, 31].

It can be shown that for $\lambda^{2}=\frac{V}{r}$ where
$V\sim\mathsf{\chi}^{2}(r)$ then
$\lambda^{2}\sim\mathsf{G}(\frac{r}{2},\frac{r}{2})$. [11] utilizes ordered
values of $\\{\frac{X_{k}}{x_{k}}\\}$ and applies the cross-entropy method to
efficiently sample large loss probabilities from a general $t$ copula $m$
factor model. Conditional on $\boldsymbol{Z}$ and $\boldsymbol{\epsilon}$, we
can arrange the order statistics of $\\{\frac{X_{k}}{x_{k}}\\}$ with
corresponding costs $\\{c_{i}\\}$. The event $\\{L>x\\}$ occurs when
$\lambda<\frac{X_{(i)}}{x_{(i)}}$ where
$i=\min\\{j:\sum_{k=j+1}^{d}c_{(k)}\leq\gamma\\}$ so we can write

$\begin{split}\mathbb{P}(L>x|Z,\boldsymbol{\epsilon})&=\mathbb{P}\left(\lambda<\frac{X_{(i)}}{x_{(i)}}\Bigg{|}Z,\boldsymbol{\epsilon}\right)\\\
&=F_{\mathsf{G}}\left(\frac{X_{(i)}^{2}}{x_{(i)}^{2}}\Bigg{|}Z,\boldsymbol{\epsilon}\right).\end{split}$

With this formulation, the cross-entropy method is applied to choose an
optimal IS density
$g^{*}(\boldsymbol{Z},\boldsymbol{\epsilon};\boldsymbol{v}^{*})$ which belongs
to the parametric density family $\mathcal{F}$ defined as

$\begin{split}\mathcal{F}=\left\\{f(\boldsymbol{Z},\boldsymbol{\epsilon};\boldsymbol{v})=\prod_{j=1}^{m}f(Z_{j};\mu_{Z},\sigma^{2}_{Z})\prod_{k=1}^{d}f(\epsilon_{k};\mu_{\boldsymbol{\epsilon}},\sigma^{2}_{\boldsymbol{\epsilon}})\right\\}\end{split}$

where
$\boldsymbol{v}^{*}=(\mu^{*}_{Z},\mathrm{\sigma}^{2*}_{Z},\mu^{*}_{\boldsymbol{\epsilon}},\sigma^{2*}_{\boldsymbol{\epsilon}})$
and
$\boldsymbol{v}=(\mu_{Z},\mathrm{\sigma}^{2}_{Z},\mu_{\boldsymbol{\epsilon}},\sigma^{2}_{\boldsymbol{\epsilon}})$.
This formulation gives the following IS estimator

$\hat{\ell}=\frac{1}{N}\sum_{k=1}^{N}F_{\mathsf{G}}\left(\frac{X_{(i)}^{2}}{x_{(i)}^{2}}\Bigg{|}\boldsymbol{Z}^{*}_{k},\boldsymbol{\epsilon}^{*}_{k}\right)\frac{f(\boldsymbol{Z}^{*}_{k},\boldsymbol{\epsilon}^{*}_{k};\boldsymbol{v})}{g^{*}(\boldsymbol{Z}^{*}_{k},\boldsymbol{\epsilon}^{*}_{k};\boldsymbol{v}^{*})}.$

####  3.2.3 Numerical Example: Gaussian and $t$ Factor Models

We illustrate CMC and IS with an example from [26] and [18] where we apply the
two-step IS algorithm for the Gaussian and CE for the $t$ factor model. Here
we consider a portfolio of size $d=1000$ under a $m=21$ factor model with
costs, marginal probabilities and thresholds as follows.

$\begin{split}c_{k}&=\left(\lceil\frac{5k}{d}\rceil\right)^{2},\\\
P_{k}&=0.01\times\left(1+\sin\left(\frac{16\pi k}{d}\right)\right),\\\
x_{\mathsf{G},k}&=\Phi^{-1}(1-P_{k}),\\\
x_{\mathsf{T}(r),k}&=F^{-1}_{\mathsf{T}(r)}(1-P_{k}),\quad
k=1,\ldots,d\end{split}$

where $x_{\mathsf{G},k}$ and $x_{\mathsf{T}(r),k}$ are the default thresholds
for the $k$-th obligor under the Gaussian and $t$ factor model with $r=3$
degrees of freedom respectively. The factor loadings matrix $\mathrm{A}$ has
the block structure

$\mathrm{A}=\begin{pmatrix}\boldsymbol{r}\begin{bmatrix}\boldsymbol{f}&&\\\
&\ddots&\\\ &&\boldsymbol{f}\end{bmatrix}\begin{array}[]{c}\mathrm{G}\\\
\vdots\\\ \mathrm{G}\end{array}\end{pmatrix},\quad\textrm{with
$\mathrm{G}=\begin{pmatrix}\boldsymbol{g}&&\\\ &\ddots&\\\
&&\boldsymbol{g}\end{pmatrix}$},$

where $\boldsymbol{r}$ is a column vector of $1000$ entries, all equal to
$0.8$; $\boldsymbol{f}$ is a column vector of $100$ entries, all equal to
$0.4$; $\mathrm{G}$ is a $100\times 10$ matrix with $\boldsymbol{g}$ a column
vector of $10$ entries, all equal to $0.4$. The conditional probabilities used
in the two-step IS algorithm are calculated as in (3.2).

Table 3.1: Estimation of risk measures for a Gaussian factor model $\alpha$ | $N$ | $\hat{\ell}^{IS}$ | $RE(\%)$ | $\hat{v}_{\alpha}$ | $\hat{v}_{\alpha}^{IS}$ | $\hat{c}_{\alpha}$ | $\hat{c}_{\alpha}^{IS}$
---|---|---|---|---|---|---|---
$0.95$ | $10^{4}$ | $0.0482$ | $1.08$ | $530$ | $548$ | $1607$ | $1650$
$0.99$ | $10^{5}$ | $0.01$ | $0.59$ | $2310$ | $2361$ | $3720$ | $3862$
$0.995$ | $10^{5}$ | $0.005$ | $0.63$ | $3376$ | $3039$ | $4863$ | $5585$
Table 3.2: Estimation of risk measures for a $t$ factor model $\alpha$ | $N$ | $\hat{\ell}^{CE}$ | $RE(\%)$ | $\hat{v}_{\alpha}$ | $\hat{v}_{\alpha}^{CE}$ | $\hat{c}_{\alpha}$ | $\hat{c}_{\alpha}^{CE}$
---|---|---|---|---|---|---|---
$0.95$ | $10^{4}$ | $0.0480$ | $1.26$ | $388$ | $352$ | $2144$ | $1934$
$0.99$ | $10^{5}$ | $0.0106$ | $0.48$ | $2934$ | $3072$ | $4733$ | $5171$
$0.995$ | $10^{5}$ | $0.0060$ | $0.77$ | $4272$ | $4684$ | $5931$ | $6539$

We first use CMC to generate a sample of size $N$ for our initial estimates
$\hat{v}_{\alpha}$ and $\hat{c}_{\alpha}$. We proceed to apply the two-step IS
algorithm where we set $\hat{v}_{\alpha}$ as the loss threshold $\gamma$. In
effect, we aim to simulate around the $\alpha$-VaR for the significance level
specified in the above table. To apply the two-step IS algorithm, we have used
the tail bound approximation to find the shifted mean vector for the risk
factors $\boldsymbol{\mu}^{*}$. Note that the sample size has been increased
for higher levels of $\alpha$ to allow for generation of more loss values in
the upper tail of the loss distribution. After obtaining our parameter
estimates for the two-step IS algorithm and CE, we have run $10$ iterations of
both methods; with each iteration generating an elite sample of size $10^{4}$
to calculate $\hat{\ell}^{IS}$ and $\hat{\ell}^{CE}$. The mean and relative
error of the $10$ values of $\hat{\ell}^{IS}$ and $\hat{\ell}^{CE}$ are
calculated with the results shown in Table 3.1. From Table 3.1, we can see
that the mean of the values $\\{\hat{\ell}^{IS}\\}$ and
$\\{\hat{\ell}^{CE}\\}$ are close to the desired values of $1-\alpha$ with RE
of $0.59\%$ to $1.08\%$ and $0.48\%$ to $1.26\%$ respectively.

###  3.3 Copula Models

A copula is defined as a multivariate distribution in the form

$C\left(u_{1},\dots,u_{d}\right)=\mathbb{P}\left(U_{1}\leq
u_{1},\dots,U_{d}\leq u_{d}\right),$

where $U_{1},\cdots,U_{d}$ are marginal uniformly distributed variables. Now
$U_{1},\cdots,U_{d}$ can be written with respect to random variables
$X_{1},\cdots,X_{n}$ with marginal distributions $F_{1},\cdots,F_{d}$. Thus we
have

$\left(U_{1},\dots,U_{d}\right)=\left(F_{1}(X_{1}),\dots,F_{d}(X_{d})\right).$

Hence the dependency of $\\{X_{k}\\}$ can be described individually through
their marginal distributions by setting

$\left(X_{1},\dots,X_{d}\right)=\left(F_{1}^{-1}(U_{1}),\dots,F_{n}^{-1}(U_{d})\right).$

We focus our discussion on a popular class of copulas known as Archimedean
Copulas. Archimedean copulas have the following form

$C\left(u_{1},\dots,u_{d}\right)=\psi^{-1}\left(\psi(u_{1})+\cdots+\psi(u_{d})\right),$

where the function $\psi:[0,1]\rightarrow[0,\infty]$ is strictly decreasing
with $\psi(0)=\infty$, $\psi(1)=0$ and its inverse $\psi^{-1}$ monotonic. This
class of copulas includes the Gumbel copula, where $\psi_{\eta}(u)=\left(-\log
u\right)^{\eta}$, and the Clayton copula, where
$\psi_{\eta}(u)=\frac{1}{\eta}(u^{-\eta}-1)$. Note that the Gumbel copula has
dependence in the upper tail while the Clayton copula has dependence in the
lower tail.

In the case of the Archimedean copulas, we can simulate the requisite vector
$\boldsymbol{U}=(U_{1},\ldots,U_{d})$ as follows. First, simulate $\Lambda\geq
0$ from the distribution $F_{\Lambda}(\lambda)$, where the pdf $f_{\Lambda}$
has the Laplace transform

$\int_{0}^{\infty}f_{\Lambda}(\lambda)\exp(-u\lambda)\mathrm{d}\lambda=\psi^{-1}(u),\quad
u\geq 0.$

Note that $\psi^{-1}(u)$, with $\psi^{-1}(0)=1$ and $\psi^{-1}(\infty)=0$, is
then a decreasing completely monotone function, as required for the copula
definition. Then, given $\Lambda$, simulate
$E_{1},\ldots,E_{n}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathsf{Exp}(1)$
and output

$(U_{1},\ldots,U_{d})=\left(\psi^{-1}\left(\frac{E_{1}}{\Lambda}\right),\ldots,\psi^{-1}\left(\frac{E_{d}}{\Lambda}\right)\right),$
$(X_{1},\ldots,X_{d})=\left(F_{1}^{-1}\left(\psi^{-1}\left(\frac{E_{1}}{\Lambda}\right)\right),\ldots,F_{d}^{-1}\left(\psi^{-1}\left(\frac{E_{d}}{\Lambda}\right)\right)\right),$

Here, probabilities are conditionally independent on $\Lambda=\lambda$, that
is $\Psi=\Lambda$, the default probability of the $k$-th obligor
$P_{k}(\Lambda)$ is

$\begin{split}P_{k}(\Lambda)&=\mathbb{P}(X_{k}>x_{k}|\Lambda=\lambda)\\\
&=\mathbb{P}(E_{k}<\lambda\psi(F_{k}(x_{k})))\quad\textrm{ since $\psi$ is
invertible and decreasing}\\\ &=1-\exp(\lambda\psi(F_{k}(x_{k}))).\end{split}$

where $\\{x_{k}\\}$ are fixed thresholds similar to those present in factor
models. Note that the nature of thresholds in copula models are often default
times [27] and although this is also possible to model under factor models,
the thresholds under factor models are more easily interpreted as economic and
market risk thresholds.

####  3.3.1 Crude Monte Carlo

Suppose that the distribution $F_{\Lambda}$ is known. For a chosen sample size
$N$ and confidence level $\alpha$, the CMC algorithm is as follows

Algorithm 6 : Generating $L$ under a Archimedean copula model using CMC

0: cost vector $\boldsymbol{c}$; default thresholds $\boldsymbol{x}$, sample
size $N$; confidence level $\alpha$; loss threshold $x$; distribution
$F_{\Lambda}$

$t\leftarrow 0$

repeat

$t\leftarrow t+1$

Simulate $\Lambda\sim F_{\Lambda}$

Simulate $E_{i}\sim\mathsf{Exp}(1)$ for $i=1,\ldots,d$

$\boldsymbol{U}\leftarrow\frac{\boldsymbol{E}}{\Lambda}$

$X_{k}\leftarrow F_{k}^{-1}(U_{k})$ for $k=1,\ldots,d$

$B_{k}\leftarrow\mathbb{I}\\{X_{k}>x_{k}\\}$ for $k=1,\ldots,d$

$L_{t}\leftarrow\boldsymbol{c}^{\top}\boldsymbol{B}$

until $t=N$

Sort the losses $\boldsymbol{L}\leftarrow(L_{[1]},\ldots,L_{[N]})^{\top}$

return $\hat{\ell}(x)=\frac{1}{N}\sum_{k=1}^{N}\mathbb{I}(L_{k}>x)$ or
$\hat{v}_{\alpha}\leftarrow L_{\lceil\alpha\times N\rceil}$

####  3.3.2 Importance Sampling

Under an Archimedean copula model, $\\{X_{k}\\}$ are independent conditional
on $\Lambda=\lambda$ with probabilities

$P_{k}(X_{k}>x_{k})=1-\exp(\lambda\psi(F_{k}(x_{k}))).$

Thus, we can apply the one-step IS algorithm similar to the case of the
Gaussian factor model by now generating the default indicators as
$B_{k}\sim\mathsf{B}(P_{k}(\Lambda))$. We then proceed by applying exponential
twisting to the probabilities $P_{k}(\Lambda)$ with twisting parameter
$\theta$, computing $L$ and the likelihood ratio

$W(L)=\exp\left(-\theta
L+\sum_{k=1}^{m}\log[1+P_{k}(\boldsymbol{Z})(\exp(\theta c_{k})-1)]\right)$

for a sample size $N$. Our IS estimator for $\ell$ remains unchanged

$\hat{\ell}(\gamma)=\frac{1}{N}\sum_{t=1}^{N}\mathbb{I}(L_{t}>\gamma)W(L_{t}).$

####  3.3.3 Numerical Example: a Clayton Copula Model

We now illustrate CMC and IS for a Clayton copula model with marginal
$\mathsf{Exp}(1)$ factors. The copula parameter $\eta$, default indicators and
costs are given as follows

$\begin{split}\eta&=5.5,\\\ B_{k}&=\mathbb{I}(X_{k}>3),\\\ c_{k}&=1\quad
k=1,\ldots,d.\end{split}$

Table 3.3: Estimation of risk measures for a Clayton copula model $\alpha$ | $N$ | $\hat{\ell}^{IS}$ | $RE(\%)$ | $\hat{v}_{\alpha}$ | $\hat{v}_{\alpha}^{IS}$ | $\hat{c}_{\alpha}$ | $\hat{c}_{\alpha}^{IS}$
---|---|---|---|---|---|---|---
$0.95$ | $10^{4}$ | $0.0476$ | $5.30$ | $55$ | $55$ | $94$ | $92$
$0.99$ | $10^{5}$ | $0.0093$ | $2.92$ | $121$ | $119$ | $162$ | $158$
$0.995$ | $10^{5}$ | $0.0041$ | $5.16$ | $154$ | $146$ | $197$ | $187$

We first use CMC to first generate a sample of size $N$ and used to give the
estimates $\hat{v}_{\alpha}$ and $\hat{c}_{\alpha}$. We proceed to apply the
one-step IS algorithm where we use $\hat{v}_{\alpha}$ as the loss threshold
$\gamma$. In effect, we aim to simulate around the $\alpha$-VaR for the
significance level specified in the above table. Note that the sample size has
once again, been increased for higher levels of $\alpha$ to allow for
generation of more loss values in the upper tail of the loss distribution.
$10$ iterations of the one-step algorithm are then run; with each iteration
generating an elite sample of size $10^{4}$ to calculate $\hat{\ell}^{IS}$.
The mean and relative error of the $10$ values of $\hat{\ell}^{IS}$ are
calculated with the results shown in Table 3.2. From Table 3.2, we can see
that although the mean of the values $\\{\hat{\ell}^{IS}\\}$ are close to the
desired values of $1-\alpha$, the relative error of the one-step IS algorithm
is $2.92\%$ to $5.30\%$ as opposed to $0.59\%$ to $1.08\%$ from the two-step
IS algorithm. From this, we can see that the shift in factors for the Gaussian
factor does indeed provide greater variance reduction for
$\\{\hat{\ell}^{IS}\\}$. We now wish to compare this performance with our
proposed algorithm, which we refer to as the _Dynamic Splitting Method_.

##  Chapter 4 Splitting Simulation

We now assume that $S$ is a quasi-monotone function, that is, $x_{i}\leq
y_{i},\;i=1,\ldots,d$ implies $S(\mathbf{x})\leq S(\mathbf{y})$. In other
words, a quasi-monotone function is one in which an increase in one of its
components cannot reduce its value. It turns out that a large number of
applied models either possess this property, or can be transformed into models
possessing it. In the absence of any special properties of $f$ or $S$ in
(1.1), there is little hope that one can do better than MCMC.

The purpose of this chapter is to show how we can simulate efficiently from
(1.1) whenever the importance function $S$ is a _quasi-monotone_ function. $S$
is a quasi-monotone function if $x_{i}\leq y_{i},\;i=1,\ldots,d$, implies
$S(\mathbf{x})\leq S(\mathbf{y})$. In other words, a quasi-monotone function
is one in which an increase in one of its components cannot reduce its value.
It turns out that a large number of applied models either possess this
property, or can be transformed into models possessing it.

The proposed algorithm is an ingenious application of the DS algorithm for
simulation of Markov processes conditional on a rare event [23]. In its
original form, DS cannot be applied to simulate from (1.1), because there is
no underlying Markov process that we can split. Our algorithm can be viewed as
a way of transforming the problem of simulation from (1.1) to one which
involves the simulation of a Markov process. The underlying idea is to embed
the static density (1.1) within a continuous time Markov process in such a way
that, at a particular instant of time, the Markov process has the exact same
density as (1.1). Given this embedding, we can then apply the original
splitting method of [23].

We emphasize that, unlike _generalized splitting_ [7, 6] and _subset
simulation_ [2], the classical splitting algorithm does not employ MCMC
sampling.

###  4.1 Classical Splitting Method

In classical splitting, we consider a Markov process $\\{\mathbf{X}_{t},t\geq
0\\}$ with the importance function $S$ over a state space $\mathscr{X}$. It is
assumed that $S(\mathbf{X}_{0})=0$ and that for any threshold $\gamma>0$,
there are unique entry times to the sets $\\{S(\mathbf{X}_{t})\geq\gamma\\}$
and $\\{S(\mathbf{X}_{t})\leq 0\\}$ respectively given by

$\tau_{\gamma}=\min\\{t:S(\mathbf{X}(t))\geq\gamma\\},$
$\tau_{0}=\min\\{t:S(\mathbf{X}(t))\leq 0\\}.$

Note that $\tau_{0}$ only exists in the absence of the quasi-monotonic
property in $S$, which is the case in classical splitting. Here the
probability of interest is $\ell=\mathbb{P}(\tau_{\gamma}<\tau_{0})$; the
probability that the process reaches the threshold $\gamma$ before reaching
$0$. Hence, $\ell$ depends on the distribution of $\mathbf{X}_{0}$. Suppose
there exists thresholds $\gamma_{1}$ and $\gamma_{2}$ such that
$\gamma_{2}>\gamma_{1}$. The cornerstone of the splitting method is the
observation that the process must first reach $\gamma_{1}$ to reach
$\gamma_{2}$, and thus, giving us a sequence of nested event subsets

$E_{\gamma_{2}}=\\{\tau_{\gamma_{2}}<\tau_{0}\\}\subset
E_{\gamma_{1}}=\\{\tau_{\gamma_{1}}<\tau_{0}\\}.$

Hence
$\ell=\mathbb{P}(E_{\gamma_{2}}|E_{\gamma_{1}})\mathbb{P}(E_{\gamma_{1}})$, a
product of conditional probabilities. This can similarly be extended to more
threshold levels, say $0=\gamma_{0}<\gamma_{1}\cdots<\gamma_{L}=\gamma$ giving
us

$E_{\gamma_{0}}\supseteq E_{\gamma_{1}}\supseteq\cdots\supseteq
E_{\gamma_{L}}.$

Let $c_{i}=\mathbb{P}(E_{\gamma_{i}}|E_{\gamma_{i-1}})$ then we have
$\ell=\prod_{i=1}^{L}c_{i}$. Each $c_{i}$ is then estimated as follows.

Define
$\mathscr{X}_{i}=\\{\mathbf{X}_{t_{i}}:S(\mathbf{X}_{t_{i}})\geq\gamma_{i}\\}$,
that is, the set of states in the Markov process that reach threshold
$\gamma_{i}$ at time $t_{i}$ for $t_{1},\cdots,t_{L}$. At each
$t_{i}\in\\{t_{1},\cdots,t_{L}\\}$, we run $s_{i}$ copies of the Markov
process $\\{\mathbf{X}_{t_{i-1}}\\}$ for each
$\mathbf{X}_{t_{i-1}}\in\mathscr{X}_{i-1}$; giving $s_{i}|\mathscr{X}_{i-1}|$
copies of $\\{\mathbf{X}_{t_{i}}\\}$ and the corresponding evolved process
$\\{S(\mathbf{X}_{t_{i}})\\}$. Each copy of $\\{\mathbf{X}_{t_{i-1}}\\}$ is
run until $\\{S(\mathbf{X}_{t_{i}})\\}$ either reaches the set
$\\{S(\mathbf{X}_{t_{i}})\geq\gamma_{i}\\}$ or $\\{S(\mathbf{X}_{t_{i}})\leq
0\\}$ with ending state $\mathbf{X}_{t_{i}}$. Each state $\mathbf{X}_{t_{i}}$
that reaches the set $\\{S(\mathbf{X}_{t_{i}})\geq\gamma_{i}\\}$ before the
set $\\{S(\mathbf{X}_{t_{i}})\leq 0\\}$, referred to as an entrance state, is
then stored as the set $\mathscr{X}_{i}$. It is from these states that the
next iteration of the splitting method will begin from; giving rise to
$s_{i+1}|\mathscr{X}_{i}|$ sample paths. An unbiased estimate of $c_{i}$ is
given by $\hat{c}_{i}=\frac{|\mathscr{X}_{i}|}{s_{i}|\mathscr{X}_{i-1}|}$. It
is clear that $|\mathscr{X}_{i}|$ is dependent on the entrance states
$\mathscr{X}_{i-1}$. Despite this dependence, [1, 17, 26] notes that the
following estimate remains unbiased

$\begin{split}\hat{\ell}&=\prod_{i=1}^{L}\frac{|\mathscr{X}_{i}|}{s_{i}|\mathscr{X}_{i-1}|}\\\
&=\frac{|\mathscr{X}_{L}|}{|\mathscr{X}_{0}|}\prod_{i=1}^{L}\frac{1}{s_{i}}.\end{split}$

Note we begin splitting at $t_{1}$ rather than $t_{0}$. In the above
framework, the integer-valued splitting factors $\\{s_{i}\\}$ are distinct and
thus may vary for each $t_{i}\in\\{t_{1},\cdots,t_{L}\\}$. This need not be
the case as a predetermined splitting factor $s_{i}=s$ for all $i$, may be
used; this version of the algorithm is referred to as _Fixed Factor
Splitting_. We will use this in combination with _Fixed Effort Splitting_
described later for our results in Chapter 4.3. The splitting process is
repeated for all $t_{i}\in\\{t_{1},\cdots,t_{L}\\}$ with
$s_{i}|\mathscr{X}_{i-1}|$ being the _simulation effort_ at $t_{i}$. Potential
problems when using the splitting method are large growths in the simulation
effort and inefficiency. These problems come from inappropriate choices in the
number of levels $L$, intermediate threshold levels
$\\{\gamma_{1},\ldots,\gamma_{L-1}\\}$ and splitting factors
$\\{s_{1},\ldots,s_{L}\\}$. [26] notes that ideally, levels should be chosen
in such a way that the conditional probabilities $\\{c_{i}\\}$ can be easily
estimated with CMC. Under the assumption of independence between the
computational cost and time from running the Markov process, the total
simulation effort is a random variable with expected value

$\begin{split}\sum_{i=1}^{L}s_{i}\mathbb{E}[|\mathscr{X}_{i-1}|]&=\sum_{i=1}^{L}s_{i}N_{0}\prod_{j=1}^{i-1}c_{j}s_{j}\\\
&=N_{0}\sum_{i=1}^{L}\frac{1}{c_{i}}\prod_{j=1}^{i}c_{j}s_{j}\end{split}$

If $c_{j}s_{j}>1$ for all $j$, the simulation effort would become large as it
increases with the number of levels $L$. This phenomenon is referred to as an
explosion in [16, 26]. If $c_{j}s_{j}<1$ for all $j$, most sample paths will
not reach the threshold levels $\gamma_{j}$ and as a consequence, the
algorithm will be inefficient. Thus, an ideal choice is
$s_{j}=\frac{1}{c_{j}}$ for all $j$. An alternative approach to avoiding
explosions is _Fixed Effort Splitting_ , where the simulation effort is fixed
for all $t_{i}\in\\{t_{1},\cdots,t_{L}\\}$ to say N, giving us a corresponding
estimator

$\hat{\ell}_{\mathrm{FE}}=\prod_{i=1}^{L}\frac{|\mathscr{X}_{i}|}{N}.$

Our results in Section 4.3 will be based on this approach. Now that a basic
description of the classical splitting method is given, we present the
workings of our proposed algorithm including the embedding of the static
density $\eqref{cond}$ in the Markov process $\mathbf{X}_{t}$.

###  4.2 The Dynamic Splitting Method for Static Problems

We first require a way to apply the classical splitting method of [23] to the
problem of sampling from (1.1). To achieve this, we induce an artificial
Markov process whose paths we can then repeatedly split to encourage entry
into a desired rare-event set, say $\\{\mathbf{X}:S(\mathbf{X})\geq\gamma\\}$
.

A suitable Markov process for our algorithm is the multivariate Lévy
subordinator [3]. A $d$-dimensional Lévy subordinator
$\\{\boldsymbol{\Lambda}(t),t\in\mathbb{R}_{+}\\}$ with
$\boldsymbol{\Lambda}(0)=0$ is an almost surely increasing stochastic process
on a probability space $(\Omega,\mathbb{P},\mathcal{F})$ with a continuous
index set on $\mathbb{R}_{+}$ and with a continuous state space
$\mathbb{R}^{d}$ defined by the following properties: (1) the increments of
$\\{\boldsymbol{\Lambda}(t)\\}$ are stationary and non-negative, that is,
$(\boldsymbol{\Lambda}(t+s)-\boldsymbol{\Lambda}(t))\geq\mathbf{0}$ has the
same distribution as $\boldsymbol{\Lambda}(s)\geq\mathbf{0}$ for all $t,s\geq
0$; we denote the density of this stationary distribution as
$\nu_{s}(\boldsymbol{\lambda})$; (2) the increments of
$\\{\boldsymbol{\Lambda}(t)\\}$ are independent, that is,
$\boldsymbol{\Lambda}(t_{i})-\boldsymbol{\Lambda}(t_{i-1}),\;i=1,2,\ldots$ are
independent for any $0\leq t_{0}<t_{1}<t_{2}<\cdots$; and (3) for any
$\epsilon>0$, we have
$\mathbb{P}(\|\boldsymbol{\Lambda}(t+s)-\boldsymbol{\Lambda}(t)\|\geq\epsilon)=0$
as $s\downarrow 0$. The distributional properties of the Lévy subordinator are
characterized by the _characteristic exponent_ which is expressed as the
logarithm of the characteristic function of the random column vector
$\boldsymbol{\Lambda}(1)$:

$\log\mathbb{E}[\exp(\mathrm{i}\boldsymbol{s}^{\top}\boldsymbol{\Lambda}(1))]=\mathrm{i}\boldsymbol{s}^{\top}\boldsymbol{\mu}+\int_{\mathbb{R}^{d}}\left(\exp(\mathrm{i}\boldsymbol{s}^{\top}\mathbf{x})-1-\mathrm{i}\boldsymbol{s}^{\top}\mathbf{x}\,\mathbb{I}_{\\{\|\mathbf{x}\|\leq
1\\}}\right)\nu(\operatorname{d\\!}\mathbf{x})\;,\quad\boldsymbol{s}\in\mathbb{R}^{d},$

for some $\boldsymbol{\mu}\in\mathbb{R}^{d}$, and measure $\nu$ such that
$\nu(\\{\mathbf{0}\\})=0$ and
$\int_{\mathbb{R}^{d}}\min\\{1,\|\mathbf{x}\|^{2}\\}\nu(\operatorname{d\\!}\mathbf{x})<\infty$.

One of the simplest multivariate subordinators we can have is the _gamma
process with independent components_ [3, Section 3.2] in which the components
$\Lambda_{1}(t),\ldots,\Lambda_{d}(t)$ of vector
$\boldsymbol{\Lambda}(t)=(\Lambda_{1}(t),\ldots,\Lambda_{d}(t))^{\top}$ are
independent one-dimensional subordinators with characteristic exponent
$\log((1-is)^{-1})$. In other words, each $\Lambda(t)$ has gamma distribution
with shape parameter $t$ and scale parameter $1$, which denote by
$\mathsf{G}(t,1)$. Note that other multivariate subordinators are possible as
long as they possess the quasi-monotonicity property, but the Gamma process
with independent components will suffice for our illustrations.

######  Example 4.1 (Simulating Gamma process).

Consider simulating a Gamma process at distinct times.

Algorithm 7 : Simulating Gamma Process at distinct times

0: Intermediate times $0=t_{0}<t_{1}<t_{2}<\cdots<t_{L}=1$

$\boldsymbol{\Lambda}(t_{0})\leftarrow 0$

for $i=1,\ldots,L$ do

Simulate $\boldsymbol{\Lambda}^{*}(t_{i}-t_{i-1})$ independently from
subordinator distribution.

$\boldsymbol{\Lambda}(t_{i})\leftarrow\boldsymbol{\Lambda}(t_{i-1})+\boldsymbol{\Lambda}^{*}(t_{i}-t_{i-1})$

return $\boldsymbol{\Lambda}(t_{1}),\ldots,\boldsymbol{\Lambda}(t_{L})$

Figure 4.1: Simulation of 1D Gamma processes Figure 4.2: Simulation of a 3D
Gamma process

Now that we have defined one of the simplest continuous-time processes (the
gamma subordinator above), we can proceed to embed the distribution of
$\mathbf{X}=(X_{1},\ldots,X_{d})^{\top}$ in (1.1) within a continuous time
process as follows. Let $F_{k}(x)$ be the cdf of $X_{k}$ for $k=1,\ldots,d$.
Define the vector $\mathbf{X}(t)=(X_{1}(t),\ldots,X_{d}(t))^{\top}$ through
the random variables

$X_{k}(t)=F^{-1}_{k}(\exp(-\Lambda_{k}(t)))$ (4.1)

With this formulation, the non-decreasing property of the gamma subordinator
implies that we will have an almost surely partial (and in fact total)
ordering for the vectors $\mathbf{X}(t)\succ\mathbf{X}(t+s)$, meaning that
$X_{i}(t)\geq X_{i}(t+s)$ for all $i$ and $s,t>0$. In addition, it can be
easily verified that at time $t=1$ each $X_{k}(1)$ has the desired
distribution $\mathbb{P}(X_{k}(1)\leq x)=F_{k}(x)$. This observation is what
connects our induced Markov process to the static distribution of
$\mathbf{X}$. In fact, we can view the realization of the original vector
$\mathbf{X}$ as a snapshot of the state of a multivariate continuous-time
process $\\{\mathbf{X}(t),t\geq 0\\}$ at the instant $t=1$. Notice that
$\Lambda_{k}(0)=0$ so we have $X_{k}(0)=\infty$. Consequently, we begin
simulations in the set $\\{S(\mathbf{X}\geq\gamma)\\}$ and evolve the gamma
subordinator towards the set $\\{S(\mathbf{X}<\gamma)\\}$ as a result of the
ordering for the vectors $\mathbf{X}(t)\succ\mathbf{X}(t+s)$ for all $i$ and
$s,t>0$. As this property suggests, we will refer to (4.1) as a _monotonically
decreasing embedding transformation_. Our numerical experiment results will be
based on this transformation.

The above formulation now has a number of implications. First, the quasi-
monotonicity of the importance function $S(\mathbf{X})$ implies that the exit
time exists and is unique

$\tau_{\gamma}=\sup\\{t:S(\mathbf{X}(t))\geq\gamma\\}$ (4.2)

and in fact
$\mathbb{P}(S(\mathbf{X}(t))\geq\gamma)=\mathbb{P}(\tau_{\gamma}>t).$ Second,
if $0=t_{0}<t_{1}<\cdots<t_{L}=1$ is any sequence of increasing times, the
state space can be decomposed into the decreasing sequence of events:

$\mathbb{R}^{d}\equiv\\{S(\mathbf{X}(t_{0}))\geq\gamma\\}\supseteq\\{S(\mathbf{X}(t_{1}))\geq\gamma\\}\cdots\supseteq\\{S(\mathbf{X}(t_{L}))\geq\gamma\\}$

and therefore we have the following decomposition of (1.2)

$\displaystyle\ell=\mathbb{P}(S(\mathbf{X}(1))\geq\gamma)$
$\displaystyle=\prod_{i=1}^{L}\mathbb{P}(S(\mathbf{X}(t_{i}))\geq\gamma|S(\mathbf{X}(t_{i-1}))\geq\gamma)$
$\displaystyle=\prod_{i=1}^{L}\mathbb{P}(\tau_{\gamma}>t_{i}|\tau_{\gamma}>t_{i-1})\,.$

At each intermediate point in time $t_{i}\in\\{t_{1},\cdots,t_{L-1}\\}$ for
$i=1,\dots,L-1$, we consider $s$ splits in $s^{i-1}$ sample paths of the gamma
subordinator for each component of $\boldsymbol{\Lambda}(t_{i})$. From this,
it is clear the estimation of each
$\mathbb{P}(\tau_{\gamma}>t_{i}|\tau_{\gamma}>t_{i-1})$ can be done
empirically as

$\mathbb{P}(\tau_{\gamma}>t_{i}|\tau_{\gamma}>t_{i-1})=\frac{|\mathscr{X}_{i}|}{s|\mathscr{X}_{i-1}|}.$

Note that since we begin the splitting of the sample paths at $t=t_{1}$ so we
have

$\mathbb{P}(\tau_{\gamma}>t_{1}|\tau_{\gamma}>t_{0})=\frac{|\mathscr{X}_{1}|}{|\mathscr{X}_{0}|}.$

$\ell$ will thus be calculated as

$\begin{split}\ell&=\frac{|\mathscr{X}_{1}|}{|\mathscr{X}_{0}|}\prod_{i=1}^{L-1}\frac{|\mathscr{X}_{i+1}|}{s|\mathscr{X}_{i}|}\\\
&=\frac{|\mathscr{X}_{L}|}{s^{L-1}|\mathscr{X}_{0}|}\\\
&=\frac{|\mathscr{X}_{L}|}{s^{L-1}}\quad\textrm{if
$|\mathscr{X}_{0}|=1$}.\end{split}$

For completeness, we will also define the _monotonically increasing embedding
transformation_. We consider the vector
$\mathbf{X}(t)=(X_{1}(t),\ldots,X_{d}(t))^{\top}$ with random variables is
defined as

$X_{k}(t)=F^{-1}_{k}(1-\exp(-\Lambda_{k}(t))).$ (4.3)

This transformation utilizes the observation that for any Uniformly
distributed random variable $U$, $1-U$ is also uniform. This comes from the
observation that $\Lambda_{k}(1)\sim\mathsf{Exp}(1)$ so its corresponding cdf
value $1-\exp(-\Lambda_{k}(1))$ is Uniformly distributed and thus,
$\exp(-\Lambda_{k}(1))$ is also Uniformly distributed. It can be easily
verified that (4.3) will also provide us with the desired distribution
$\mathbb{P}(X_{k}(1)\leq x)=F_{k}(x)$. However, the consequences of this
transformation differs from those of (4.1).

The non-decreasing property of the Gamma subordinator will now lead to
increases in $X_{k}(t)$ as $t$ increases. As a result, we will now have the
partial or total ordering for the vectors $\mathbf{X}(t)\prec\mathbf{X}(t+s)$,
that is, $X_{i}(t)\leq X_{i}(t+s)$ for all $i$ and $s,t>0$. As
$\Lambda_{k}(0)=0$ we will have $X_{k}(0)=\min\\{-\infty,0\\}$ depending on
the marginal distribution $F_{k}$. We will assume $X_{k}(0)=0$ for our
purposes as we aim our discussion to the models described in Chapter 3. We
note that all random variables that require this transformation in these
models are indeed non-negative. Under this transformation we have
$S(\mathbf{X}(0))=0$ and begin the Markov process in the set
$\\{S(\mathbf{X}(t))<\gamma\\}$. We wish to estimate $\ell$ by modelling the
entry time to the set $\\{S(\mathbf{X}(t))\geq\gamma\\}$

$\tau_{\gamma}=\min\\{t:S(\mathbf{X}(t))\geq\gamma\\},$

Hence, under the transformation (4.3) we arrive at the classical splitting
method as described in section 4.1 with the modification that $\tau_{0}$ does
not exist as a consequence of the quasi-monotonicity of the importance
function and that we now seek

$\displaystyle\ell=\mathbb{P}(S(\mathbf{X}(1))\geq\gamma)$
$\displaystyle=\prod_{i=1}^{L}\mathbb{P}(S(\mathbf{X}(t_{i}))\geq\gamma|S(\mathbf{X}(t_{i-1}))\geq\gamma)$
$\displaystyle=\prod_{i=1}^{L}\mathbb{P}(\tau_{\gamma}<t_{i}|\tau_{\gamma}<t_{i-1}).$

We can summarize the fixed factor dynamic splitting algorithm as follows.

Figure 4.3: Splitting of paths Algorithm 8 : Fixed Factor Dynamic Splitting,
returning $W$, an unbiased estimate of $\ell$

0: Splitting factor $s$ and intermediate times $0<t_{1}<t_{2}\cdots<t_{L}=1$

Generate
$\boldsymbol{\Lambda}(t_{1})=(\Lambda_{1}(t_{1}),\ldots,\Lambda_{d}(t_{1}))^{\top}$
from the subordinator distribution.

if $S(\boldsymbol{\Lambda}(t_{1}))>\gamma$ then

$\mathscr{X}_{1}\leftarrow\\{\boldsymbol{\Lambda}(t_{1})\\}$

else

return $W\leftarrow 0$

for $i=2,\ldots,L$ do

if $\mathscr{X}_{i-1}=\emptyset$ then

return $W\leftarrow 0$

else

$\mathscr{X}_{i}\leftarrow\emptyset$

for all $\boldsymbol{\Lambda}(t_{i-1})\in\mathscr{X}_{i-1}$ do

for $j=1,\ldots,s$ do

For $k=1,\ldots,d$ sample independently
$\Lambda_{k}^{*}(t_{i}-t_{i-1})\sim\mathsf{G}(t_{i}-t_{i-1},1).$

$\boldsymbol{\Lambda}^{*}(t_{i}-t_{i-1})\leftarrow(\Lambda_{1}^{*}(t_{i}-t_{i-1}),\ldots,\Lambda_{d}^{*}(t_{i}-t_{i-1}))^{\top}$

$\boldsymbol{\Lambda}(t_{i})\leftarrow\boldsymbol{\Lambda}(t_{i-1})+\boldsymbol{\Lambda}^{*}(t_{i}-t_{i-1})$

if $S(\boldsymbol{\Lambda}(t_{i}))>\gamma$ then

add $\boldsymbol{\Lambda}(t_{i})$ to $\mathscr{X}_{i}$

return $W\leftarrow|\mathscr{X}_{L}|/s^{L-1}$ as an unbiased estimate.

In the above formulation, the splitting factor $s$ is chosen arbitrarily but
under two idealizing assumptions, a near optimal value can be chosen for $s$
(see Appendix A.2).

As noted previously, 8 has the risk of explosions so we implement a _fixed
effort_ variant of the above algorithm to avoid this. At each $t_{i}$, we fix
the simulation effort to $s$ splits for a randomly chosen
$\boldsymbol{\Lambda}(t_{i-1})\in\mathscr{X}_{i-1}$. The fixed effort dynamic
splitting algorithm can be summarized as follows.

Algorithm 9 : Fixed Effort Dynamic Splitting, returning $W$, an unbiased
estimate of $\ell$

0: Total sample at each level $s$; intermediate times
$0<t_{1}<t_{2}\cdots<t_{L}=1$

$\mathscr{X}_{1}\leftarrow\emptyset$

for $k=1,\ldots,s$ do

Generate
$\boldsymbol{\Lambda}(t_{1})=(\Lambda_{1}(t_{1}),\ldots,\Lambda_{d}(t_{1}))^{\top}$
from the subordinator distribution.

if $S(\boldsymbol{\Lambda}(t_{1}))>\gamma$ then

Add $\boldsymbol{\Lambda}(t_{1})$ to $\mathscr{X}_{1}$

for $i=2,\ldots,L$ do

if $|\mathscr{X}_{i-1}|=0$ then

return $W\leftarrow 0$

else

$\mathscr{X}_{i}\leftarrow\emptyset$

for $j=1,\ldots,s$ do

Let $\boldsymbol{\Lambda}(t_{i-1})$ be a randomly chosen member of
$\mathscr{X}_{i-1}$.

For $k=1,\ldots,d$ sample independently
$\Lambda_{k}^{*}(t_{i}-t_{i-1})\sim\mathsf{G}(t_{i}-t_{i-1},1).$

$\boldsymbol{\Lambda}^{*}(t_{i}-t_{i-1})\leftarrow(\Lambda_{1}^{*}(t_{i}-t_{i-1}),\ldots,\Lambda_{d}^{*}(t_{i}-t_{i-1}))^{\top}$

$\boldsymbol{\Lambda}(t_{i})\leftarrow\boldsymbol{\Lambda}(t_{i-1})+\boldsymbol{\Lambda}^{*}(t_{i}-t_{i-1})$

if $S(\boldsymbol{\Lambda}(t_{i}))>\gamma$ then

add $\boldsymbol{\Lambda}(t_{i})$ to $\mathscr{X}_{i}$

return $W\leftarrow\prod_{i=1}^{L}|\mathscr{X}_{i}|/s^{L}$ as an estimate of
$\ell$.

###  4.3 Numerical Experiments

We now illustrate the performance of dynamic splitting (DS) in the credit risk
models described in Chapter 3 through numerical experiments. For these
results, we first verify that the importance function
$S(\mathbf{X})=L(\mathbf{X})$ is indeed quasi-monotonic and specify suitable
transformation (4.1) and (4.3) to embed the static random variables
$\\{X_{k}\\}$ into time-dependent random variables $\\{X_{k}(t)\\}$. We first
note that all of these models have the form

$\begin{split}S(\mathbf{X})&=L(\mathbf{X})\\\
&=\boldsymbol{c}^{\top}\boldsymbol{B}(\mathbf{X})\end{split}$

where the $k$-th component of $\boldsymbol{B}$ is generated as
$\mathbb{I}(X_{k}>x_{k})$. An increase in one of the components of
$\mathbf{X}$ will lead to the value of $L(\mathbf{X})$ either staying the same
or increased due to the activation of the corresponding indicator variable. It
is clear that the function value of $L(\mathbf{X})$ will not be reduced so the
importance function $S(\mathbf{X})=L(\mathbf{X})$ does indeed possess the
quasi-monotonicity property required for our algorithm. We now present the
results of the application of our algorithm to the _Gaussian factor model_ ,
_$t$ factor model_ and _Clayton copula model_. In the following results, we
have applied two-step IS for the Gaussian factor model, the CE-based estimator
from [11] for the $t$ factor model and one-step IS for the Clayton copula
module. We note that the threshold $\gamma$ is chosen with an initial run of
IS before our splitting algorithm is applied.

In the following results, $\gamma$ is the $\alpha$-quantile of the loss
distribution estimated with IS from Chapter 3 while $\hat{\ell}^{IS}$ and
$\hat{\ell}^{DS}$ are the probabilities of the a loss exceeding $\gamma$.
Hence, ideally we expect $\alpha+\hat{\ell}^{IS}$ and $\alpha+\hat{\ell}^{DS}$
to both equal to $1$. Disregarding the initial simulation effort for $\gamma$,
all algorithms implemented below have been constrained to a simulation effort
of $10^{5}$ and the same number of runs $R=10$. For IS and CE estimators, this
is simply $NR$ where $N=10^{4}$ is the sample size in each run of the IS
algorithm. For the DS estimator, it is $sTR$ where $s=1000$ is the splitting
factor and $T=10$ is the number of intermediate time levels. Note, we have
made appropriate adjustments to $s$ and $T$ to match the simulation effort of
both algorithms in each comparison below.

####  4.3.1 Factor Models

We illustrate the performance of the splitting algorithm with factor models.
We assume that $\mathrm{A}$ has the matrix structure given in our worked
example

$\mathrm{A}=\begin{pmatrix}\boldsymbol{r}\begin{bmatrix}\boldsymbol{f}&&\\\
&\ddots&\\\ &&\boldsymbol{f}\end{bmatrix}\begin{array}[]{c}\mathrm{G}\\\
\vdots\\\ \mathrm{G}\end{array}\end{pmatrix},\quad\textrm{with
$\mathrm{G}=\begin{pmatrix}\boldsymbol{g}&&\\\ &\ddots&\\\
&&\boldsymbol{g}\end{pmatrix}$},$

where $\boldsymbol{r}$ is a column vector of $1000$ entries, all equal to
$0.8$; $\boldsymbol{f}$ is a column vector of $100$ entries, all equal to
$0.4$; $\mathrm{G}$ is a $100\times 10$ matrix with $\boldsymbol{g}$ a column
vector of $10$ entries, all equal to $0.4$. $\boldsymbol{b}$ is calculated as

$b_{k}=\sqrt{1-(a_{k1}^{2}+\cdots+a_{km}^{2})}$

#####  4.3.1.1 Gaussian Factor Model

We have

$X_{k}=a_{k1}Z_{1}+\cdots+a_{km}Z_{m}+b_{k}\epsilon_{k}$

where
$Z_{1},\cdots,Z_{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathsf{N}(0,1)$
and $\epsilon_{k}\sim\mathsf{N}(0,1)$. By using (4.1) we have

$X_{k}(t)=a_{k1}\Phi^{-1}(\Lambda_{1}(t))+\cdots+a_{km}\Phi^{-1}(\Lambda_{m}(t))+b_{k}\Phi^{-1}(\Lambda_{m+k}(t)).$

Here we require gamma processes for $m$ standard Normal random variables which
are the systematic risk factors $\boldsymbol{Z}$ and $d$ standard Normal
random variables for the obligor-specific risk factors
$\boldsymbol{\epsilon}$.

Table 4.1: Estimation of $\ell$ for a Gaussian factor model $\alpha$ | $\gamma$ | $\hat{\ell}^{IS}$ | $\hat{\ell}^{DS}$ | $RE^{IS}(\%)$ | $RE^{DS}(\%)$
---|---|---|---|---|---
$0.95$ | $548$ | $0.0493$ | $0.0497$ | $0.73$ | $4.35$
$0.99$ | $2361$ | $0.0098$ | $0.0103$ | $0.60$ | $6.82$
$0.995$ | $3039$ | $0.0062$ | $0.0056$ | $0.58$ | $7.28$

We can see that both $\alpha+\hat{\ell}^{IS}$ and $\alpha+\hat{\ell}^{DS}$ do
indeed equal to $1$ approximately. However, the relative error of the DS
estimator has a much larger relative error $RE^{DS}$ than the relative error
for the two-step IS estimator $RE^{IS}$. Thus, in this study the DS estimator
does not perform as efficiently as the two-step IS estimator.

#####  4.3.1.2 $t$ Factor Model

We have

$\begin{split}X_{k}&=\sqrt{\frac{r}{V}}\left(a_{k1}Z_{1}+\cdots+a_{km}Z_{m}+b_{k}\epsilon_{k}\right)\\\
&=\sqrt{\frac{1}{G}}\left(a_{k1}Z_{1}+\cdots+a_{km}Z_{m}+b_{k}\epsilon_{k}\right)\end{split}$

where
$Z_{1},\cdots,Z_{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathsf{N}(0,1)$,
$\epsilon_{k}\sim\mathsf{N}(0,1)$ and
$G\sim\mathsf{G}(\frac{v}{2},\frac{v}{2})$. By using (4.1) and (4.3) we have

$X_{k}(t)=\frac{\left(a_{k1}\Phi^{-1}(\exp(-\Lambda_{1}(t)))+\cdots+a_{km}\Phi^{-1}(\exp(-\Lambda_{m}(t)))+b_{k}\Phi^{-1}(\exp(-\Lambda_{m+k}(t)))\right)}{\sqrt{F_{\mathsf{G}}^{-1}(1-\exp(-\Lambda_{d+m}(t)))}}.$

Here we require gamma processes for $m$ standard Normal random variables which
are the systematic risk factors $\boldsymbol{Z}$, $d$ standard Normal random
variables for the obligor-specific risk factors $\boldsymbol{\epsilon}$ and
another for the common random variable $G$.

Table 4.2: Estimation of $\ell$ for a $t$ factor model $\alpha$ | $\gamma$ | $\hat{\ell}^{CE}$ | $\hat{\ell}^{DS}$ | $RE^{CE}(\%)$ | $RE^{DS}(\%)$
---|---|---|---|---|---
$0.95$ | $352$ | $0.0500$ | $0.0478$ | $0.36$ | $3.43$
$0.99$ | $3072$ | $0.0100$ | $0.0096$ | $1.00$ | $6.55$
$0.995$ | $4684$ | $0.0050$ | $0.0047$ | $0.72$ | $5.76$

The performance of the CE estimator for the $t$ factor model has similar
performance efficiency as that of the two-step IS estimator for the Gaussian
factor model. This supports our theoretical motivations of these estimators
since both estimators aim to significantly reduce the variance of $\hat{\ell}$
(see Chapter 2) where CE algorithm aims to sample approximately from the zero-
variance density by learning near optimal parameters for the likelihood
function while two-step IS reduces the variance via exponential twisting and
change in parameter measures. Hence, we expect the relative error of both
estimators to be quite small as supported by the above empirical results. We
observe that $\alpha+\hat{\ell}^{IS}$ and $\alpha+\hat{\ell}^{DS}$ equal to
$1$ approximately. However, the relative error of the DS estimator remains
larger than the relative error for the CE estimator $RE^{CE}$. Thus, in this
study the DS estimator does not perform as efficiently as the one-step IS
estimator.

####  4.3.2 Clayton Copula Model

In this example, we consider the variables $\\{X_{k}\\}$ with marginal
exponential distributions $\mathsf{Exp}(1)$. Note that $\\{X_{k}\\}$ are
static variables computed as

$\begin{split}X_{k}&=F_{k}^{-1}(U_{k})\\\
&=-\log\left(1-\psi^{-1}\left(\frac{E_{k}}{G}\right)\right)\end{split}$

where $F_{k}$ is the marginal distribution of $X_{k}$, $U_{k}$ is a uniform
variable generated from the Clayton copula,
$G\sim\mathsf{G}(\frac{1}{\eta},1)$ and $\psi^{-1}(t)=(1+\eta
t)^{\frac{1}{\eta}}$. Let $F_{\eta}$ be the cdf of the
$\mathsf{G}(\frac{1}{\eta},1)$ then an appropriate embedding transformation is

$X_{k}(t)=-\log\left(1-\psi^{-1}\left(\frac{\Lambda_{k}(t)}{F_{\eta}^{-1}(\exp(-\Lambda_{d+1}(t)))}\right)\right).\
$

where $\Lambda_{k}(t),\Lambda_{d+1}(t)\sim\mathsf{G}(t,1)$. Here we require
gamma processes for $d$ independent $\mathsf{Exp}(1)$ random variables and
another for the common random variable $G$.

Table 4.3: Estimation of $\ell$ for a Clayton copula model $\alpha$ | $\gamma$ | $\hat{\ell}^{IS}$ | $\hat{\ell}^{DS}$ | $RE^{IS}(\%)$ | $RE^{DS}(\%)$
---|---|---|---|---|---
$0.95$ | $55$ | $0.0493$ | $0.0513$ | $1.83$ | $4.87$
$0.99$ | $119$ | $0.0094$ | $0.0099$ | $4.14$ | $6.67$
$0.995$ | $146$ | $0.0049$ | $0.0047$ | $5.49$ | $8.14$

Once again, we can see that $\alpha+\hat{\ell}^{IS}$ and
$\alpha+\hat{\ell}^{DS}$ equal to $1$ approximately. However, the relative
error of the DS estimator is still larger than the relative error for the one-
step IS estimator $RE^{IS}$. Thus, in this study the DS estimator does not
perform as efficiently as the one-step IS estimator.

###  4.4 Critical Analysis of the Dynamic Splitting Method

In this section, we give a critical analysis of the proposed DS method when
applied to rare-event probability estimation. The proposed DS method has
advantages in terms of simplicity, versatility and simulation effort.
Conversely, the DS approach can be inapplicable if the problem does not
possess a quasi-monotonic importance function $S$.
Advantages:

  1. 1.

Interpretability $-$ In the DS method, we estimate $\ell$ by computing the
product of sequential conditional probabilities. Each conditional probability
can be interpreted as a survival probability from the current time level to
the next time level. This gives us a simple interpretation to our estimate of
$\ell$ as a survival probability from the starting time to the terminal time.

  2. 2.

Versatility $-$ One of the biggest strengths of the proposed DS method is that
it is versatile as it can generate any set of continuous random variables
$\mathbf{X}$ by using the inverse-transform method. In particular, we can
generate sequential values of each variable such that
$\\{S(\mathbf{X})\geq\gamma\\}$. This is possible as the monotonicity of the
Gamma process is preserved from under the embedding transformation, thus
allowing us generate the set $\\{\mathbf{X}:S(\mathbf{X})\geq\gamma\\}$. We
also note that unlike IS algorithms, DS does not require the computation of
likelihood ratios and can be applied on any continuous pdf including heavy-
tailed distributions. This is a major advantage for DS when compared to the
one-step and two-step IS algorithms which are exclusively applicable to light-
tailed distributions such as the Normal distribution. The same argument cannot
hold against CE as the only main restriction in CE is that it only considers
densities within the family of distributions which is not a strong assumption
in itself.

  3. 3.

Simulation effort $-$ A natural question that would arise is the simulation
effort and computational cost of DS since it simulates split paths of the
Markov chain over many time levels. The danger of explosions can be avoided by
a fixed effort implementation and paths $\mathbf{X}$ for which
$\\{\mathbf{X}:S(\mathbf{X})<\gamma\\}$ are not simulated further. We also
note that unlike generalized splitting or particle methods, our proposed DS
method does not require approximate MCMC sampling at each level of splitting.
These observations ensure that the implementation of our proposed DS method
would not be computationally expensive.

Limitations:

  1. 1.

Quasi-monotonicity $-$ Our proposed DS algorithm hinges on the existence of a
quasi-monotonic importance function $S$. If the problem of interest has an
importance function $S$ that does not possess this property then we cannot
apply DS as the absence of quasi-monotonicity in $S$ implies we can no longer
apply the embedding transformations to the problem and the decomposition of
1.2 into nested subsets no longer holds. Consequently, there is no connection
between the static distributions of interest and a time-dependent Markov
process, making DS inapplicable.

  2. 2.

Efficiency $-$ As shown by the results of numerical experiments earlier in
this chapter, we can see that the RE of the DS method is much larger than
current hallmark rare-event probability estimation algorithms such two-step IS
and CE based algorithms [11]. This would be attributed to the decomposition of
time levels to estimate each conditional probability to be inadequate form of
variance reduction in comparison to two-step IS and CE, both of which seek to
obtain a density with reduced variance for efficient estimation of $\ell$.

##  Chapter 5 Conclusions and Suggestions for Future Research

We have developed and illustrated a splitting method which is simple and
effective to estimate rare event probabilities in the framework of popular
copula credit risk models. This method is designed to estimate tail
probabilities by decomposing the state space of the risk variables into nested
subsets so that the rare event can be expressed as an intersection of these
subsets. This decomposition allows us to achieve greater accuracy in
estimating each conditional probability than estimating the rare probability
itself. We have also shown that despite the inapplicability of the classical
splitting method on static problems, one can always find an appropriate
embedding transformation such that the static density being modelled can be
taken as a snapshot of a continuous time Markov process at a particular
instance in time. The illustrations of importance sampling estimators in this
thesis, namely _exponential twisting_ , have relied on conditional
independence of the model to simplify the problem to a sum of scaled Bernoulli
random variables and attractive forms of the likelihood ratio. Our proposed
modification to the dynamic splitting algorithm does not rely on conditional
independence, rather it relies on the quasi-monotonicity property of the
importance function. A large number of models either possess this property or
can achieve it through a transformation. Hence, although the algorithm does
not perform as efficiently as current hallmark importance sampling estimators,
it is a versatile method that is applicable to any static problem for which
random variables required in the model have prespecified distributions.

A possible direction for future research include developing more efficient
adaptive dynamic splitting algorithms and variance reduction techniques to
dynamic splitting, perhaps by applying exponential twisting to each candidate
sample path that will be split. The latter is based on the observation that
each sample path that is split in fact Binomial distributed with $s$ trials,
that is $\mathsf{Bin}(s,c)$. Possible practical applications for future work
may include financial and insurance risk management. This is motivated by the
observation that the splitting framework can be interpreted as a model for
survival probabilities which is analogous to survival probabilities of
insurance policyholders and insurers; and also default and bankruptcy rates of
financial institutions. Specifically, it may be useful in the pricing of
insurance products as well as risk and solvency capital allocation where the
VaR and CVaR are possible risk measures.

##  Appendix A Monte Carlo methods

###  A.1 One-step IS algorithm: Shift in Probabilities

We simulate $L$ under the following model. Let

$L=c_{1}B_{1}+\cdots+c_{m}B_{m}=\boldsymbol{c}^{\top}\boldsymbol{B}$

be the total loss incurred by a portfolio of $m$ obligors, where $c_{i}$ is
the loss incurred from the $i$-th obligor and $B_{i}\sim\mathsf{B}(P_{i})$ is
a Bernoulli random variable indicating whether the $i$-th obligor has
defaulted. The distribution of the vector $\boldsymbol{B}$ is implicitly
defined via

$B_{k}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbb{I}\\{X_{k}>x_{k}\\},\quad
k=1,\ldots,m$

where $\\{x_{k}\\}$ are given fixed thresholds and $\boldsymbol{X}$ has a
continuous joint density $f(\boldsymbol{x})$. Thus, the likehihood of the
total loss $L$ is given by

$f_{L}(l)=\prod_{k=1}^{m}P_{k}^{b_{k}}(1-P_{k})^{1-b_{k}}$

To simulate from the upper tail of the distribution, an exponential twist is
applied to the probabilities $\\{P_{i}\\}$ for $\theta\geq 0$ as the first
step in importance sampling. It is defined as follows

$P_{k,\theta}=\frac{P_{k}e^{\theta c_{k}}}{1+P_{k}(e^{\theta c_{k}}-1)}.$

This change in probability measure results in a density that simulates higher
values of $L$ with default indicators now generated by
$B_{i}\sim\mathsf{B}(\tilde{P}_{i})$. Note that the notation $L$ and $B_{k}$
have been preserved as the desired likelihood ratio should be described under
the notation of the original distribution. As the probability of a particular
default outcome is being compared across two densities, the likelihood ratio
aims to describe the ratio of probabilities for the same outcome. The
likelihood of $L$ under the exponentially twisted density $g_{L}$ is given by

$g_{L}(l)=\prod_{k=1}^{m}\tilde{P}_{k}^{b_{k}}(1-\tilde{P}_{k})^{1-b_{k}}.$

Hence the likelihood ratio $W$ is given by

$\begin{split}W(l)&=\frac{f_{L}(l)}{g_{L}(l)}\\\
&=\prod_{k=1}^{m}\left(\frac{P_{k}}{\tilde{P}_{k}}\right)^{b_{k}}\left(\frac{1-P_{k}}{1-\tilde{P}_{k}}\right)^{1-b_{k}}\\\
&=\prod_{k=1}^{m}\left(\frac{1+P_{k}(e^{\theta c_{k}}-1)}{e^{\theta
c_{k}}}\right)^{b_{k}}(1+P_{k}(e^{\theta c_{k}}-1))^{1-b_{k}}\\\
&=\prod_{k=1}^{m}\left(1+P_{k}(e^{\theta c_{k}}-1)\right)e^{-\theta
b_{k}c_{k}}\\\ &=\exp\left(-\theta
l+\sum_{k=1}^{m}\log\left(1+P_{k}(\exp(\theta
c_{k})-1)\right)\right)\end{split}$

###  A.2 Two-step IS algorithm: Shift in Factors

As the second step in importance given in [18], the mean of the factors
$\boldsymbol{Z}$ are shifted from $\boldsymbol{0}$ to $\boldsymbol{\mu}$.
Assuming at Normal copula, this means we now have
$\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{\mu},\mathrm{I}_{d})$ instead of
$\boldsymbol{Z}\sim\mathsf{N}(\boldsymbol{0},\mathrm{I}_{d})$. Hence, the
change in mean must also form part of the estimator. We note that the
normalization constants of Normal densities are unchanged by a shift in the
mean so the likelihood ratio is given by

$\frac{\exp{(-\frac{1}{2}\boldsymbol{Z}^{\top}\boldsymbol{Z})}}{\exp{(-\frac{1}{2}\boldsymbol{(\mu-Z)}^{\top}\boldsymbol{(\mu-Z)})}}=\exp{(\frac{1}{2}\boldsymbol{\mu}^{\top}\boldsymbol{\mu}-\boldsymbol{\mu}^{\top}\boldsymbol{Z})}$

###  A.3 Archimedean copula sampling algorithm

To verify that the algorithm in [25] draws
$\boldsymbol{U}=(U_{1},\cdots,U_{d})$ from an Archimedean copula, we must show
that

$P(\boldsymbol{U}\leq\boldsymbol{u})=\psi^{-1}\left(\sum_{i=1}^{d}\psi(u_{i})\right),$

where

$(U_{1},\ldots,U_{d})=\left(\psi^{-1}\left(\frac{E_{1}}{\Lambda}\right),\ldots,\psi^{-1}\left(\frac{E_{d}}{\Lambda}\right)\right).$

The probability of $\boldsymbol{U}$ can be written as

$\begin{split}\mathbb{P}(\boldsymbol{U}\leq\boldsymbol{u})&=\mathbb{P}(U_{1}\leq
u_{1},\ldots,U_{d}\leq u_{d})\\\
&=\mathbb{P}(E_{1}\geq\Lambda\psi(u_{1}),\ldots,E_{d}\geq\Lambda\psi(u_{d}))\quad\textrm{
since $\psi$ is invertible and decreasing}\\\
&=\mathbb{E}_{\Lambda}\mathbb{P}(E_{1}\geq\Lambda\psi(u_{1}),\ldots,E_{d}\geq\Lambda\psi(u_{d})|\Lambda=\lambda)\\\
&=\mathbb{E}_{\Lambda}\prod_{i}\mathbb{P}(E_{i}\geq\lambda\psi(u_{i}))\\\
&=\mathbb{E}_{\Lambda}\prod_{i}\exp(-\lambda\psi(u_{i}))\\\
&=\mathbb{E}_{\Lambda}\exp\left(-\lambda\sum_{i=1}^{d}\psi(u_{i})\right)\\\
&=\int_{0}^{\infty}\exp\left(-\lambda\sum_{i=1}^{d}\psi(u_{i})\right)f_{\Lambda}(\lambda)\mathrm{d}\lambda\\\
&=\psi^{-1}\left(\sum_{i=1}^{d}\psi(u_{i})\right),\qquad\textrm{as
required}.\end{split}$

##  Appendix A Dynamic Splitting

###  A.1 Embedding transformations

We now verify that the embedding transformations (4.1) and (4.3) do indeed
have the desired distribution at $t=1$. For (4.1) we have

$\begin{split}\mathbb{P}(X_{k}(1)\leq x)&=\mathbb{P}(\exp(-\Lambda_{k}(1))\leq
F_{k}(x))\\\ &=\mathbb{P}(\Lambda_{k}(1)\geq-\log(F_{k}(x)))\\\
&=\exp(\log(F_{k}(x)))\\\ &=F_{k}(x).\end{split}$

For (4.3) we have

$\begin{split}\mathbb{P}(X_{k}(1)\leq
x)&=\mathbb{P}(1-\exp(-\Lambda_{k}(1))\leq F_{k}(x))\\\
&=\mathbb{P}(\exp(-\Lambda_{k}(1))\geq 1-F_{k}(x))\\\
&=\mathbb{P}(\Lambda_{k}(1)\leq-\log(1-F_{k}(x)))\\\
&=1-\exp(\log(1-F_{k}(x)))\\\ &=F_{k}(x).\end{split}$

Hence both transformations yield the desired distribution $F_{k}(x)$ at $t=1$
as required.

###  A.2 Ideal Case Analysis of Fixed Factor Splitting Algorithm

We now present an analysis of the performance of the Fixed Factor Splitting
algorithm under an ideal assumption. The assumption is that the time levels
$t_{i}\in\\{t_{1},\cdots,t_{L}\\}$ are selected such that the conditional
probabilities $c=\mathbb{P}(S(\mathbf{X}(t_{i}))|S(\mathbf{X}(t_{i-1}))$ are
exactly, rather than approximately, equal to $s$ for all $i$.

Let $N_{i}=|\mathscr{X}_{i}|$ be the random number of states in the set
$\mathscr{X}_{i}=\\{X(t_{i}):S(X(t_{i}))\geq\gamma\\}$. At time $t_{0}=0$ we
have $N_{1}=1$ as the algorithm begins with a single path for the Markov
process $\boldsymbol{\Lambda}(t_{1})$. If we denote the number of states in
$\mathscr{X}_{i+1}$ that are generated from the $j$-th state from
$\mathscr{X}_{i}$ by $Q_{j,i}$ then we have the branching process recursion

$N_{i+1}=Q_{1,i}+Q_{2,i}+\cdots+Q_{N_{t},i}$

where it is clear that $Q_{j,i}\sim\mathsf{Bin}(s,c)$, that is, a Binomial
distributed random variable with probability $c$. Thus, we have
$\mathbb{E}[Q_{j,i}]=sc=1$ and $\mathrm{Var}(Q_{j,i})=sc(1-c)=1-c$. By
standard branching process arguments [20] we have $\mathbb{E}[N_{i}]=1$ and
$\mathrm{Var}(N_{i})=(i-1)(1-c)$ for $1<i<L$. Hence, for the unbiased
estimator $W=\frac{|\mathscr{X}_{L}|}{s^{L-1}}=\frac{N_{L}}{s^{L-1}}$ we have
$\mathbb{E}[W]=\ell=\frac{1}{s^{L-1}}$ and
$\mathrm{Var}(W)=\frac{(L-1)(1-c)}{s^{2L-2}}$ with $\log(\ell)=(L-1)\log(s)$.

An estimator $\hat{\ell}$ of $\ell$ is _logarithmically efficient_ [26] if the
following condition holds:

$\limsup_{\ell\downarrow
0}\left|\frac{\log(\mathrm{Var}(\hat{\ell}))}{\log(\ell^{2})}\right|\geq 1.$

For the logarithmic efficiency criterion we have

$\begin{split}\lim_{\ell\downarrow
0}\left|\frac{\log(\mathrm{Var}(W))}{\log(\ell^{2})}\right|&=\lim_{\ell\downarrow
0}\left|\frac{\log(L-1)+\log(1-c)-(2L-2)\log(s)}{(2L-2)\log(s)}\right|\\\
&=\lim_{s\uparrow\infty}\left|\frac{\log(L-1)+\log(1-c)-(2L-2)\log(s)}{(2L-2)\log(s)}\right|\\\
&=1.\end{split}$

Therefore, under the idealized assumption the estimator $W$ is logarithmically
efficient. Note that the simulation effort, starting from $t_{1}$, is a random
variable $s\sum_{i=1}^{L}N_{i}$ with expected value $s(L-1)$. The expected
_relative time variance product_ [26] is thus given by

$\begin{split}\frac{\mathrm{Var}(W)}{\ell^{2}}s(L-1)&=(L-1)(1-c)s(L-1)\\\
&=(L-1)^{2}(s-1)\\\
&=\left(\frac{\log\ell}{\log(s)}\right)^{2}(s-1),\end{split}$

which is minimized as a function of $s$ for $s>1$ at $s=4.92155363$ or $s=5$
when constrained on the integers.

###  A.3 Adaptive Dynamic Splitting

The description of the dynamic splitting method described up to this point has
allowed arbitrary choices for the intermediate time levels
$t_{i}\in\\{t_{1},\cdots,t_{L-1}\\}$. We now describe a pilot algorithm to
select optimal values for the intermediate time levels $\\{t_{i}\\}$. The
motivation behind this selection is to ensure that the conditional
probabilities
$\mathbb{P}(S(\boldsymbol{\Lambda}(t_{i}))\geq\gamma|S(\boldsymbol{\Lambda}(t_{i-1}))\geq\gamma)$
are not too small and not rare-event probabilities so that they can be easily
estimated with CMC. [7] notes that this formulation will in fact lead to
biased estimates of $\ell$ and complications in the computation of the
variance of $\hat{\ell}$. As a consequence, the relative error of the
estimator can be difficult to compute, thereby limiting comparisons and
benchmarks against other Monte Carlo estimators. Thus, the _Fixed Effort
Splitting_ algorithm is recommended as it will lead to unbiased estimates for
$\ell$ and readily available estimates of variance by running the algorithm
several times independently. The following algorithm implements a simple
procedure similar to numerical root-finding by utilizing the Gamma Bridge
Sampling algorithm.

Algorithm 10 : Adaptive Dynamic Splitting

0: start time $t_{s}$; end time $t_{e}$; endpoints
$\boldsymbol{\Lambda}(t_{s})\leq\boldsymbol{\Lambda}(t_{e})$; total sample at
each level $s$; time tolerance $\epsilon_{t}$; proportion tolerance
$\epsilon_{p}$

$\mathscr{T}\leftarrow\emptyset$

$\boldsymbol{\Lambda}(t_{l})\leftarrow\boldsymbol{\Lambda}(t_{s})$

$\boldsymbol{\Lambda}(t_{u})\leftarrow\boldsymbol{\Lambda}(t_{e})$

$t_{m}\leftarrow\frac{t_{s}+t_{e}}{2}$

for $j=1,\ldots,s$ do

Generate
$\boldsymbol{\Lambda}_{j}(t_{m})=(\Lambda_{1}(t_{m}),\ldots,\Lambda_{d}(t_{m}))^{\top}$
using Gamma Bridge Sampling with endpoints
$\boldsymbol{\Lambda}(t_{l})\leq\boldsymbol{\Lambda}(t_{u})$

$i\leftarrow 1$

while $t_{e}-t_{m}>\epsilon_{t}$ do

while
$|\frac{1}{s}\sum_{j=1}^{s}\mathbb{I}\\{S(\boldsymbol{\Lambda}_{j}(t_{m}))\geq\gamma\\}-\rho|>\epsilon_{p}$
do

if
$\frac{1}{s}\sum_{j=1}^{s}\mathbb{I}\\{S(\boldsymbol{\Lambda}(t_{m})\geq\gamma\\}-\rho>0$
then

$\boldsymbol{\Lambda}(t_{l})\leftarrow\boldsymbol{\Lambda}(t_{m})$

$t_{m}\leftarrow\frac{t_{m}+t_{u}}{2}$

else if
$\frac{1}{s}\sum_{j=1}^{s}\mathbb{I}\\{S(\boldsymbol{\Lambda}(t_{m})\geq\gamma\\}-\rho<0$
then

$\boldsymbol{\Lambda}(t_{u})\leftarrow\boldsymbol{\Lambda}(t_{m})$

$t_{m}\leftarrow\frac{t_{l}+t_{m}}{2}$

for $j=1,\ldots,s$ do

Generate
$\boldsymbol{\Lambda}_{j}(t_{m})=(\Lambda_{1}(t_{m}),\ldots,\Lambda_{d}(t_{m}))^{\top}$
using Gamma Bridge Sampling with endpoints
$\boldsymbol{\Lambda}(t_{l})\leq\boldsymbol{\Lambda}(t_{u})$

$t_{i}\leftarrow t_{m}$

add $t_{i}$ to $\mathscr{T}$

$i\leftarrow i+1$

return $\mathscr{T}$

Algorithm 11 : Gamma bridge sampling

0: Endpoints $\boldsymbol{\Lambda}(t_{l})\leq\boldsymbol{\Lambda}(t_{u})$ and
$t\in(t_{l},t_{u})$.

for $k=1,\ldots,d$ do

$B_{k}\sim\mathsf{Beta}(t-t_{l},t_{u}-t)$, independently

$\Lambda_{k}(t)\leftarrow\Lambda_{k}(t_{l})+(\Lambda_{k}(t_{u})-\Lambda_{k}(t_{l}))B_{k}$

return
$\boldsymbol{\Lambda}(t)=(\Lambda_{1}(t),\ldots,\Lambda_{d}(t))^{\top}$.

## References

  * [1] Søren Asmussen and Peter W Glynn.  Stochastic simulation: Algorithms and analysis, volume 57.  Springer Science & Business Media, 2007.
  * [2] Siu-Kui Au and James L. Beck.  Estimation of small failure probabilities in high dimensions by subset simulation.  Probabilistic Engineering Mechanics, 16(4):263–277, 2001.
  * [3] Ole E. Barndorff-Nielsen, Jan Pedersen, and Ken-iti Sato.  Multivariate subordination, self-decomposability and stability.  Advances in Applied Probability, pages 160–187, 2001.
  * [4] Achal Bassamboo, Sandeep Juneja, and Assaf Zeevi.  Portfolio credit risk with extremal dependence: Asymptotic analysis and efficient simulation.  Operations Research, 56(3):593–606, 2008.
  * [5] Christian Bluhm, Ludger Overbeck, and Christoph Wagner.  Introduction to credit risk modeling.  Crc Press, 2010.
  * [6] Zdravko I. Botev and Dirk P. Kroese.  An efficient algorithm for rare-event probability estimation, combinatorial optimization, and counting.  Methodology and Computing in Applied Probability, 10(4):471–505, 2008.
  * [7] Zdravko I. Botev and Dirk P. Kroese.  Efficient monte carlo simulation via the generalized splitting method.  Statistics and Computing, 22(1):1–16, 2012.
  * [8] Tim J Brereton, Dirk P Kroese, and Joshua C Chan.  Monte carlo methods for portfolio credit risk.  2013\.
  * [9] Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng.  Handbook of Markov Chain Monte Carlo.  CRC press, 2011.
  * [10] James Bucklew.  Introduction to rare event simulation.  Springer Science & Business Media, 2013.
  * [11] Joshua CC Chan and Dirk P Kroese.  Efficient estimation of large portfolio loss probabilities in t-copula models.  European Journal of Operational Research, 205(2):361–367, 2010\.
  * [12] Amir Dembo and Ofer Zeitouni.  Large deviations techniques and applications, volume 38.  Springer Science & Business Media, 2009.
  * [13] Rüdiger Frey and Alexander J McNeil.  Dependent defaults in models of portfolio credit risk.  Journal of Risk, 6:59–92, 2003.
  * [14] Paul Glasserman.  Monte Carlo methods in financial engineering, volume 53.  Springer Science & Business Media, 2003.
  * [15] Paul Glasserman, Philip Heidelberger, and Perwez Shahabuddin.  Asymptotically optimal importance sampling and stratification for pricing path-dependent options.  Mathematical finance, 9(2):117–152, 1999.
  * [16] Paul Glasserman, Philip Heidelberger, Perwez Shahabuddin, and Tim Zajic.  A large deviations perspective on the efficiency of multilevel splitting.  Automatic Control, IEEE Transactions on, 43(12):1666–1679, 1998\.
  * [17] Paul Glasserman, Philip Heidelberger, Perwez Shahabuddin, and Tim Zajic.  A look at multilevel splitting.  Springer, 1998.
  * [18] Paul Glasserman and Jingyi Li.  Importance sampling for portfolio credit risk.  Management science, 51(11):1643–1656, 2005.
  * [19] Peter W Glynn.  Importance sampling for monte carlo estimation of quantiles.  In Mathematical Methods in Stochastic Simulation and Experimental Design: Proceedings of the 2nd St. Petersburg Workshop on Simulation, pages 180–185, 1996.
  * [20] Theodore E Harris.  The theory of branching processes.  Courier Corporation, 2002.
  * [21] L. Jeff Hong, Zhaolin Hu, and Guangwu Liu.  Monte carlo methods for value-at-risk and conditional value-at-risk: A review.  ACM Transactions on Modeling and Comput. Simulation, 24(4):22:1–22:37, 2014.
  * [22] L. Jeff Hong and Guangwu Liu.  Monte carlo estimation of value-at-risk, conditional value-at-risk and their sensitivities.  In Proceedings of the Winter Simulation Conference, pages 95–107, 2011.
  * [23] Herman Kahn and Ted E. Harris.  Estimation of particle transmission by random sampling.  National Bureau of Standards applied mathematics series, 12:27–30, 1951.
  * [24] Wanmo Kang and Perwez Shahabuddin.  Fast simulation for multifactor portfolio credit risk in the t-copula model.  In Proceedings of the 37th conference on Winter simulation, pages 1859–1868, 2005.
  * [25] Dirk P. Kroese, Reuven Y. Rubinstein, and Peter W. Glynn.  The cross-entropy method for estimation.  Handbook of Statistics, edited by V. Govindaraju and CR Rao, 31, 2010.
  * [26] Dirk P. Kroese, Thomas Taimre, and Zdravko I. Botev.  Handbook of Monte Carlo Methods, volume 706.  John Wiley & Sons, 2011.
  * [27] David X Li.  On default correlation: A copula function approach.  Available at SSRN 187289, 1999.
  * [28] Reuven Y Rubinstein and Dirk P Kroese.  Simulation and the Monte Carlo method, volume 707.  John Wiley & Sons, 2011.
  * [29] Reuven Y Rubinstein and Dirk P Kroese.  The cross-entropy method: a unified approach to combinatorial optimization, Monte-Carlo simulation and machine learning.  Springer Science & Business Media, 2013.
  * [30] Alexandre Scott.  Extensions of the cross-entropy method with applications to diffusion processes and portfolio losses.  2015\.
  * [31] Alexandre Scott and Adam Metzler.  A general importance sampling algorithm for estimating portfolio loss probabilities in linear factor models.  Available at SSRN, 2015.
  * [32] Aad W Van der Vaart.  Asymptotic statistics, volume 3.  Cambridge university press, 2000.