Importance Sampling in the Presence of PD-LGD Correlation

Abstract

This paper seeks to identify computationally efficient importance sampling (IS) algorithms for estimating large deviation probabilities for the loss on a portfolio of loans. Related literature typically assumes that realised losses on defaulted loans can be predicted with certainty, i.e., that loss given default (LGD) is non-random. In practice, however, LGD is impossible to predict and tends to be positively correlated with the default rate and the latter phenomenon is typically referred to as PD-LGD correlation (here PD refers to probability of default, which is often used synonymously with default rate). There is a large literature on modelling stochastic LGD and PD-LGD correlation, but there is a dearth of literature on using importance sampling to estimate large deviation probabilities in those models. Numerical evidence indicates that the proposed algorithms are extremely effective at reducing the computational burden associated with obtaining accurate estimates of large deviation probabilities across a wide variety of PD-LGD correlation models that have been proposed in the literature.

1. Introduction

This paper seeks to identify computationally efficient importance sampling (IS) algorithms for estimating large deviation probabilities for the loss on a portfolio of loans. Related literature assumes that realised losses on defaulted loans can be predicted with certainty, i.e., that loss given default (LGD) is non-random. In practice, however, LGD is impossible to predict and tends to be positively correlated with the default rate and the latter phenomenon is typically referred to as PD-LGD correlation (here PD refers to probability of default, which is often used synonymously with default rate). There is a large literature on modelling stochastic LGD and PD-LGD correlation, but there is a paucity of literature on using importance sampling to estimate large deviation probabilities in those models. This gap in the literature was brought to our attention by a risk management professional at a large Canadian financial institution, and filling that gap is the ultimate goal of this paper.

Problem Formulation and Related Literature

Consider a portfolio of N exposures of equal size. Let $L_1,L_2,\dots ,L_N$ denote the losses on the individual loans, expressed as a percentage of notional value. The percentage loss on the entire portfolio is:

$${\overline{L}}_N:=\frac{1}{N}\sum _{i=1}^N{L}_i.$$

(1)

We are interested in using IS to estimate large deviation probabilities of the form:

$$p_x:=\mathbb{P}({\overline{L}}_N\ge x),$$

(2)

where $x>>\mathbb{E}\left[L_i\right]=\mathbb{E}\left[{\overline{L}}_N\right]$ is some large, user-defined, threshold.

In practice the number of exposures is large (e.g., in the thousands) and prudent risk management requires one to assume that the individual losses are correlated. In practice, then, ${\overline{L}}_N$ is the average of a large number of correlated variables. As such, its probability distribution is highly intractable and Monte Carlo is the method of choice for approximating $p_x$. As the probability of interest is typically small (e.g., on the order of ${10}^{-3}$ or ${10}^{-4}$), the computational burden required to obtain an accurate estimate of $p_x$ using Monte Carlo can be prohibitive. For instance if $p_x$ is on the order of ${10}^{-3}$ and N is on the order of 1000 then, in the absence of any variance reduction techniques, the sample size required to reduce the estimator’s relative error1 to 10% is on the order of one hundred thousand. Since each realisation of ${\overline{L}}_N$ requires simulation of one thousand individual losses, a sample size of $\mathrm{100,000}$ requires one to generate one hundred million variables. If the desired degree of accuracy is reduced to 1%, the number of variables that must be generated increases to a staggering 10 billion.

Importance sampling (IS) is a variance reduction technique that has the potential to significantly reduce the computational burden associated with obtaining accurate estimates of large deviation probabilities. In the present context, effective IS algorithms have been identified for a variety of popular risk management models, but most are limited to the special case that loss given default (LGD) is non-random. The seminal paper in the area is (Glasserman and Li 2005), other papers include (Chan and Kroese 2010) and (Scott and Metzler 2015). It is well documented empirically, however, that portfolio-level LGD is not only stochastic, but positively correlated with the portfolio-level default rate as seen, for instance, in any of the studies listed in (Kupiec 2008) or (Frye and Jacobs 2012). This phenomenon is typically referred to as PD-LGD correlation. (Miu and Ozdemir 2006) show that ignoring PD-LGD correlation when it is in fact present can lead to material underestimates of portfolio risk measures.

There is a large literature on modelling PD-LGD correlation (Frye 2000); (Pykhtin 2003); (Miu and Ozdemir 2006); (Kupiec 2008); (Sen 2008); (Witzany 2011); (de Wit 2016); (Eckert et al. 2016); and others listed in (Frye and Jacobs 2012), but there is a much smaller literature on using IS to estimate large deviation probabilities in such models. To the best of our knowledge only (Deng et al. 2012) and (Jeon et al. 2017) have developed algorithms that allow for PD-LGD correlation (the former paper considers a dynamic intensity-based framework, the latter considers a static model with asymmetric and heavy-tailed risk factors). The present paper contributes to this nascent literature by developing algorithms that can be applied in a wide variety of PD-LGD correlation models that have been proposed in the literature, and are popular in practice.

The paper is structured as follows. Section 2 outlines important assumptions, notation, and terminology. Section 3 theoretically motivates the proposed algorithm in a general setting, and Section 4 discusses a few practical issues that arise when implementing the algorithm. Section 5 describes a general framework for PD-LGD correlation modelling that includes, as special cases, many of the models that have been developed in the literature and Section 6 describes how to implement the proposed algorithm in this general framework. Numerical results are presented and discussed in Section 7 and demonstrate that the proposed algorithms are extremely effective at reducing the computational burden required to obtain an accurate estimate of $p_x$.

2. Assumptions, Notation and Terminology

We assume that individual losses are of the form $L_i=\mathcal{L}(\mathbf{Z},{\mathbf{Y}}_i)$, where $\mathcal{L}$ is some deterministic function, $\mathbf{Z}=(Z_1,\dots ,Z_d)$ is a d-dimensional vector of systematic risk factors that affect all exposures, and ${\mathbf{Y}}_i$ is a vector of idiosyncratic risk factors that only affect exposure i. We assume that $\mathbf{Z},{\mathbf{Y}}_1,{\mathbf{Y}}_2,\dots$ are independent, and that the ${\mathbf{Y}}_i$ are identically distributed. The primary role of the systematic risk factors is to induce correlation among the individual exposures, and it is common to interpret the realised values of the systematic risk factors as determining the overall macroeconomic environment. It is worth noting that the we do not require the components of $\mathbf{Z}$ to be independent of one another, etc. for the components of ${\mathbf{Y}}_i$.

2.1. Large Portfolios and the Region of Interest

In a large portfolio, the influence of the idiosyncratic risk factors is negligible. Indeed, since individual losses are conditionally independent, given the realised values of the systematic risk factors, we have the almost sure limit:

$$\underset{N\to \infty }{lim}{\overline{L}}_N=\mu \left(\mathbf{Z}\right),$$

(3)

where

$$\mu \left(\mathbf{z}\right):=\mathbb{E}\left[L_i\right|\mathbf{Z}=\mathbf{z}]=\mathbb{E}\left[{\overline{L}}_N\right|\mathbf{Z}=\mathbf{z}].$$

(4)

Since $\mu \left(\mathbf{Z}\right)\approx {\overline{L}}_N$ for large N by Equation (3), the random variable $\mu \left(\mathbf{Z}\right)$ is often called the large portfolio approximation (LPA) to ${\overline{L}}_N$. The LPA is often used to formalise the intuitive notion that, in a large portfolio, all risk is systematic (i.e., idiosyncratic is “diversified away”). We define the region of interest as the set:

$$\{\mathbf{z}\in {\mathbb{R}}^d:\mu \left(\mathbf{z}\right)\ge x\}.$$

(5)

The region of interest is “responsible” for large deviations in the sense that:

$$\underset{N\to \infty }{lim}\mathbb{P}(\mu \left(\mathbf{Z}\right)\ge x|{\overline{L}}_N\ge x)=1$$

(6)

for most values2 of x. Together, Equations (3) and (6) suggest that for large portfolios, it is relatively more important to identify an effective IS distribution for the systematic risk factors, as compared to the idiosyncratic risk factors.

2.2. Systematic Risk Factors

We assume that $\mathbf{Z}$ is continuous and let $f\left(\mathbf{z}\right)$ denote its joint density. We assume that f is a member of an exponential family (see Bickel and Doksum 2001 for definitions and important properties) with natural sufficient statistic $S:{\mathbb{R}}^d\mapsto {\mathbb{R}}^p$. Any other member of the family can be put in the form:

$$f_{\lambda }\left(\mathbf{z}\right):=exp({\lambda }^{T}S\left(\mathbf{z}\right)-K\left(\lambda \right))·f\left(\mathbf{z}\right),$$

(7)

where $K(·)$ is the cumulant generating function (cgf) of $S\left(\mathbf{Z}\right)$ and $\lambda \in {\mathbb{R}}^p$ is such that $K\left(\lambda \right)$ is well-defined. The parameter $\lambda$ is called the natural parameter of the family in Equation (7). Appendix B embeds the Gaussian and multivariate t families into this general framework.

We will eventually be using densities of the form in Equation (7) as IS densities for the systematic risk factors. The associated IS weight is:

$$\frac{f\left(\mathbf{Z}\right)}{f_{\lambda }\left(\mathbf{Z}\right)}=exp(-{\lambda }^{T}S\left(\mathbf{Z}\right)+K\left(\lambda \right)),$$

(8)

and it will be important to know when the variance of the IS weight is finite. The following observation is readily verified.

Remark 1.

If $\mathbf{Z}\sim f_{\lambda }$, then Equation (8) has finite variance if and only if both $K\left(\lambda \right)$ and $K(-\lambda )$ are well defined.

A standard result in the theory of exponential families is that:

$$\nabla K\left(\lambda \right)={\mathbb{E}}_{\lambda }\left[S\left(\mathbf{Z}\right)\right],$$

(9)

where ∇ denotes gradient and ${\mathbb{E}}_{\lambda }$ denotes expectation with respect to the density $f_{\lambda }$.

2.3. Individual Losses

We assume that $L_i$ takes values in the unit interval. In general $L_i$ will have a point mass at zero (if it did not, the loan would not be prudent) and the conditional distribution of $L_i$, given that $L_i>0$, is called the (account-level) LGD distribution. We allow the LGD distribution to be arbitrary in the sense that it could be either discrete or continuous, or a mixture of both. This contrasts with the case of non-random LGD, where the LGD distribution is degenerate at a single point. We let ${\ell }_{max}\in (0,1]$ denote the supremum of the support of $L_i$. Individual losses will therefore never exceed ${\ell }_{max}$ but could take on values arbitrarily close (and possibly equal) to ${\ell }_{max}$.

Remark 2.

Despite the fact that $L_i$ is not a continuous variable, in what follows we will proceed as if it was and make repeated reference to its “density." This is done without loss of generality, and in an interest of simplifying the presentation and discussion. Nothing in the sequel requires $L_i$ to be a continuous variable, and everything carries over to the case where it is either discrete or continuous, or has both a discrete and continuous component.

For $\mathbf{z}\in {\mathbb{R}}^d$ we let $g(\ell |\mathbf{z})$ denote the conditional density of $L_i$, given that $\mathbf{Z}=\mathbf{z}$. We assume that the support of $g(·|\mathbf{z})$ is identical to the unconditional support, in particular it does not depend on the value of $\mathbf{z}$. Note that $\mu \left(\mathbf{z}\right)$ is the mean of $g(·|\mathbf{z})$.

In practice (i.e., for all of the PD-LGD correlation models listed in the introduction) $g(·|\mathbf{z})$ is not a member of an established parametric family, and direct simulation from $g(·|\mathbf{z})$ using a standard technique such as inverse transform or rejection sampling is not straightforward. Simulation from $g(·|\mathbf{z})$ is most easily accomplished by simulating the idiosyncratic risk factors, ${\mathbf{Y}}_i$, from their density, say $\eta \left(\mathbf{y}\right)$, and then setting $L_i=\mathcal{L}(\mathbf{z},{\mathbf{Y}}_i)$. In other words, in order to simulate from $g(·|\mathbf{z})$ we make use of the fact that $L_i=\mathcal{L}(\mathbf{z},{\mathbf{Y}}_i)$ is a drawing from $g(·|\mathbf{z})$ whenever ${\mathbf{Y}}_i$ is a drawing from $\eta (·)$.

For $\theta \in \mathbb{R}$ and $\mathbf{z}\in {\mathbb{R}}^d$ we let:

$$k(\theta ,\mathbf{z}):=log\left(\mathbb{E}[exp\left(\theta L_i\right)|\mathbf{Z}=\mathbf{z}]\right)$$

and

$$k^{\prime }(\theta ,\mathbf{z}):=\frac{\partial k}{\partial \theta }(\theta ,\mathbf{z}).$$

Then $k(·,\mathbf{z})$ is the conditional cgf of $L_i$, given that $\mathbf{Z}=\mathbf{z}$, and $k^{\prime }(·,\mathbf{z})$ is its first derivative. In practice, neither $k(·,\mathbf{z})$ nor $k^{\prime }(·,\mathbf{z})$ is available in closed form. In the examples we consider later in the paper each can be expressed as a one-dimensional integral, but the numerical values of those integrals must be approximated using quadrature. This contrasts with the case of non-random LGD, where the conditional cgf can be computed in closed form3.

For $x\in (0,{\ell }_{max})$ and $\mathbf{z}\in {\mathbb{R}}^d$ we let $\widehat{\theta }(x,\mathbf{z})$ denote the unique solution to the equation $k^{\prime }(\theta ,\mathbf{z})=max(x,\mu \left(\mathbf{z}\right))$. We often suppress dependence on x and $\mathbf{z}$, and simply write $\widehat{\theta }$ instead of $\widehat{\theta }(x,\mathbf{z})$. That $\widehat{\theta }$ is well-defined follows immediately from the developments in Appendix A.1. Based on the discussion there we find that $\widehat{\theta }$ is zero whenever $\mathbf{z}$ lies in the region of interest, and is strictly positive otherwise.

Remark 3.

In practice, the value of $\widehat{\theta }$ cannot be computed in closed form and must be approximated using a numerical root-finding algorithm. Since each evaluation of the function $k^{\prime }(·,\mathbf{z})$ requires quadrature, computing $\widehat{\theta }$ is straightforward but relatively time consuming. This contrasts with the case of non-random LGD, where $\widehat{\theta }$ can be computed in closed form at essentially no cost.

For $\mathbf{z}\in {\mathbb{R}}^d$ we let $q(·,\mathbf{z})$ denote the Legendre transform of $k(·,\mathbf{z})$ over $[0,\infty )$. That is,

$$q(x,\mathbf{z}):=\underset{\theta \ge 0}{max}(\theta x-k(\theta ,\mathbf{z}))=\widehat{\theta }x-k(\widehat{\theta },\mathbf{z}).$$

(10)

That $\widehat{\theta }$ is the uniquely defined point at which the function $\theta \mapsto \theta x-k(\theta ,\mathbf{z})$ attains its maximum on $[0,\infty )$ follows from the developments in Appendix A.2. Based on the discussion there, we find that both $\widehat{\theta }$ and q are equal to zero whenever $\mathbf{z}$ lies in the region of interest, and that both are strictly positive otherwise.

2.4. Conditional Tail Probabilities

Given the realised values of the systematic risk factors, individual losses are independent. Large deviations theory can therefore provide useful insights into the large-N behaviour of the tail probability $\mathbb{P}({\overline{L}}_N\ge x|\mathbf{Z}=\mathbf{z})$. For instance, Chernoff’s bound yields the estimate:

$$\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\le exp(-Nq(x,\mathbf{z})),$$

(11)

and Cramér’s (large deviation) theorem yields the limit:

$$\underset{N\to \infty }{lim}\frac{log\left(\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\right)}{N}=-q(x,\mathbf{z}).$$

(12)

Together these results are often used to justify the approximation:

$$\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\approx exp(-Nq(x,\mathbf{z})),$$

(13)

which will be used repeatedly throughout the paper. The approximation in Equation (13) is often called the large deviation approximation (LDA) to the tail probability $\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})$. Note that since $q(x,\mathbf{z})=0$ whenever $\mu \left(\mathbf{z}\right)\ge x$, the LDA suggests that $\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\approx 1$ whenever $\mathbf{z}$ lies in the region of interest.

2.5. Conditional Densities

Let $\mathbf{L}=(L_1,\dots ,L_N)$, noting that $\mathbf{L}$ takes values in ${[0,{\ell }_{max}]}^N$. For $\mathbf{z}\in {\mathbb{R}}^d$ and $\ell =({\ell }_1,\dots ,{\ell }_N)\in {[0,{\ell }_{max}]}^N$, we let $h_x(\mathbf{z},\ell )$ denote the conditional density of $(\mathbf{Z},\mathbf{L})$, given that ${\overline{L}}_N>x$. Then $h_x$ is given by:

$$h_x(\mathbf{z},\ell )=\frac{f\left(\mathbf{z}\right)·{\prod }_{i=1}^{N}g\left({\ell }_i\right|\mathbf{z})}{p_x}·{\mathbf{1}}_{\{\ell \in A_{N,x}\}},$$

(14)

where $A_{N,x}$ is the set of points $\ell \in {[0,{\ell }_{max}]}^N$ for which $N^{-1}{\sum }_{i=1}^N{\ell }_i>x$.

We let $f_x\left(\mathbf{z}\right)$ denote the conditional density of the systematic risk factors, given that ${\overline{L}}_N>x$, noting that:

$$f_x\left(\mathbf{z}\right)=\frac{\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})}{\mathbb{P}({\overline{L}}_N\ge x)}·f\left(\mathbf{z}\right).$$

(15)

In the examples we consider the mean of $f_x$ tends to lie inside, but close to the boundary of, the region of interest. And relative to the unconditional density f, the conditional density $f_x$ tends to be much more concentrated about its mean.

Finally, we let $g_x(\ell |\mathbf{z})$ denote the conditional density of an individual loss, given that $\mathbf{Z}=\mathbf{z}$ and ${\overline{L}}_N>x$, noting that:

$$g_x(\ell |\mathbf{z})=\frac{\mathbb{P}({\overline{L}}_{N-1}>x+\frac{x-\ell }{N-1}|\mathbf{Z}=\mathbf{z})}{\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})}·g(\ell |\mathbf{z}).$$

(16)

If the realised value of $\mathbf{z}$ lies inside the region of interest, the conditional density $g_x(·|\mathbf{z})$ tends to resemble the unconditional density $g(·|\mathbf{z})$. Intuitively, for such values of $\mathbf{z}$ the LDA informs that the event $\{{\overline{L}}_N>x\}$ is very likely, and conditioning on its occurrence is not overly informative. If the realised value of $\mathbf{z}$ does not lie in the region of interest then $g_x(·|\mathbf{z})$ tends to resemble the exponentially tilted version of $g(·|\mathbf{z})$ whose mean is exactly x. See Appendix A.3 for more details.

Neither $h_x,f_x$, nor $g_x$ are numerically tractable, but as we will soon see they do serve as useful benchmarks against which to compare candidate IS densities. In addition, it is worth noting here that the representations of Equations (15) and (16) lend themselves to numerical approximation via the LDA in Equation (13).

3. Proposed Algorithm

In practice, the most common approach to estimating $p_x$ via Monte Carlo simulation in this framework is summarised in Algorithm 1 below.

Algorithm 1 Standard Monte Carlo Algorithm for Estimating $p_x$

1: Simulate M i.i.d. copies of the systematic risk factors. Think of these as different economic scenarios and denote the simulated values by ${\mathbf{z}}_1,\dots ,{\mathbf{z}}_M$.   2: For each scenario m: (a) Simulate the idiosyncratic risk factors for each exposure. Denote the simulated values ${\mathbf{y}}_{1,m},\dots ,{\mathbf{y}}_{N,m}$. (b) Set ${\ell }_{i,m}=\mathcal{L}({\mathbf{z}}_m,{\mathbf{y}}_{i,m})$ for each exposures i, and ${\overline{\ell }}_m=\frac{1}{N}{\sum }_{i=1}^N{\ell }_{i,m}$. 3: Return ${\widehat{p}}_x=\frac{1}{M}{\sum }_{m=1}^M{\mathbf{1}}_{\{{\overline{\ell }}_m>x\}}$.

Algorithm 1 consists of two stages. In the first stage one simulates the systematic risk factors, and in the second stage one simulates the idiosyncratic risk factors for each exposure. Mathematically, the first stage induces independence among the individual exposures, so that the second stage amounts to simulating a large number of i.i.d. variables. Intuitively, it is useful to think of the first stage as determining the prevailing macroeconomic environment, which fixes economy-wide quantities such as default and loss-given-default rates. The second stage of the algorithm overlays idiosyncratic noise on top of economy-wide rates, to arrive at the default and loss-given-default rates for a particular portfolio.

Relative error is the preferred measure of accuracy for estimators of rare event probabilities. The relative error of the estimator ${\widehat{p}}_x$ in Algorithm 1 is:

$$\frac{1}{\sqrt{M}}\sqrt{\frac{1-p_x}{p_x}},$$

and the sample size required to ensure the relative error does not exceed some predetermined threshold $ϵ$ is:

$$M\left(ϵ\right)=\frac{1}{{ϵ}^2}\frac{1-p_x}{p_x}.$$

(17)

The number of variables that must be generated in order to achieve the desired degree of accuracy $ϵ$ is therefore $(N+d)·M\left(ϵ\right)$, which grows without bound as $p_x\to 0$. For instance if $p_x={10}^{-3}$, $N={10}^3$, $d=2$, and $ϵ=5·{10}^{-2}$ then the number of variables that must be generated is approximately four hundred million, which is an enormous computational burden for a modest degree of accuracy. In the next section we discuss general principles for selecting an IS algorithm that can reduce the computational burden required to obtain an accurate estimate of $p_x$.

3.1. General Principles

For practical reasons, we insist that our IS procedure retains conditional independence of individual losses, given the realised value of the systematic risk factors. This is important because it allows us to reduce the problem of simulating a large number of dependent variables to the (much) more computationally efficient problem of simulating a large number of independent variables.

In the first stage we simulate the systematic risk factors from the IS density $f_{IS}\left(\mathbf{z}\right)$. The IS weight associated with this first stage is therefore:

$${\mathsf{\Lambda }}_1\left(\mathbf{z}\right):=\frac{f\left(\mathbf{z}\right)}{f_{IS}\left(\mathbf{z}\right)}.$$

In the second stage we simulate the individual losses as i.i.d. drawings from the density $g_{IS}(\ell |\mathbf{z})$. The IS weight associated with this second stage is:

$${\mathsf{\Lambda }}_2(\mathbf{z},\ell )=\prod _{i=1}^N\frac{g\left({\ell }_i\right|\mathbf{z})}{g_{IS}\left({\ell }_i\right|\mathbf{z})},$$

and the IS density from which we sample $(\mathbf{Z},\mathbf{L})$ is therefore of the form:

$$h_{IS}(\mathbf{z},\ell )=f_{IS}\left(\mathbf{z}\right)·\prod _{i=1}^N{g}_{IS}\left({\ell }_i\right|\mathbf{z}).$$

(18)

The so-described algorithm, with as-yet unspecified IS densities, is summarised in Algorithm 2.

Algorithm 2 IS Algorithm for Estimating $p_x$

1: Simulate M i.i.d. copies of the systematic risk factors from the density $f_{IS}\left(\mathbf{z}\right)$. Think of these as different economic scenarios and denote the simulated values by ${\mathbf{z}}_1,\dots ,{\mathbf{z}}_M$.   2: For each scenario m: (a) Independently simulate ${\ell }_{1,m},{\ell }_{2,m},\dots ,{\ell }_{N,m}$ from the density $g_{IS}(·|{\mathbf{z}}_m)$.   (b) Set ${\overline{\ell }}_m=\frac{1}{N}{\sum }_{i=1}^N{\ell }_{i,m}$. 3: Return ${\widehat{p}}_x=\frac{1}{M}{\sum }_{m=1}^M{\mathsf{\Lambda }}_1\left({\mathbf{z}}_m\right)·{\mathsf{\Lambda }}_2({\mathbf{z}}_m,{\ell }_m)·{\mathbf{1}}_{\{{\overline{\ell }}_m>x\}}$, where ${\ell }_m=({\ell }_{1,m},\dots ,{\ell }_{N,m})$.

It is important to note that in the second stage, we will not be simulating individual losses directly from the (conditional) IS density $g_{IS}$. Rather, we will simulate the idiosyncratic risk factors ${\mathbf{Y}}_i$ in such a way as to ensure that for a given value of $\mathbf{z}$, the variable $L_i=\mathcal{L}(\mathbf{z},{\mathbf{Y}}_i)$ has the desired density $g_{IS}$. Focusing on the “indirect" IS density of $L_i$, as opposed to “direct" IS density of ${\mathbf{Y}}_i$, allows us to identify a much more effective second stage algorithm4.

The estimator ${\widehat{p}}_x$ produced by Algorithm 2 is demonstrably unbiased and its variance is:

$${\mathbb{E}}_{IS}[{(\mathsf{\Lambda }(\mathbf{Z},\mathbf{L})·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}-p_x)}^2]=p_x^2·{\mathbb{E}}_{IS}[{({\mathsf{\Lambda }}_x(\mathbf{Z},\mathbf{L})·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}-1)}^2],$$

(19)

where ${\mathbb{E}}_{IS}$ denotes expectation under the IS distribution, $\mathsf{\Lambda }(\mathbf{z},\ell ):={\mathsf{\Lambda }}_1\left(\mathbf{z}\right)·{\mathsf{\Lambda }}_2(\mathbf{z},\ell )$ and

$${\mathsf{\Lambda }}_x(\mathbf{z},\ell ):=\frac{\mathsf{\Lambda }(\mathbf{z},\ell )}{p_x}.$$

Note that ${\mathsf{\Lambda }}_x$ is the ratio of (i) the IS density in Equation (18) to (ii) the conditional density in Equation (14). The estimator’s squared relative error can then be decomposed as:

$${\mathbb{E}}_{IS}[{({\mathsf{\Lambda }}_x(\mathbf{Z},\mathbf{L})-1)}^2·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}]+[1-{\mathbb{P}}_{IS}({\overline{L}}_N>x)],$$

(20)

where ${\mathbb{P}}_{IS}$ denotes probability under the IS distribution.

Inspecting Equation (20) we see that an effective IS density should (i) assign a high probability to the event of interest and (ii) should resemble the conditional density in Equation (14) as closely as possible, in the sense that the ratio ${\mathsf{\Lambda }}_x$ should deviate as little as possible from unity. Clearly, an estimator that satisfies (ii) should also satisfy (i), since $h_x$ assigns probability one to the event that ${\overline{L}}_N>x$. The task now is to identify a density of the form in Equation (18) that resembles the ideal density in Equation (14), in some sense.

3.2. Identifying the Ideal IS Densities

Our measure of similarity is Kullback–Leibler divergence (KLD), or divergence for short. See Chatterjee and Diaconis (2018) for a general discussion of the merits of minimum divergence as a criteria for identifying effective IS distributions. We begin by writing:

$$\frac{h_x(\mathbf{z},\ell )}{h_{IS}(\mathbf{z},\ell )}=\frac{f_x\left(\mathbf{z}\right)}{f_{IS}\left(\mathbf{z}\right)}·\frac{{\tilde{g}}_x(\ell |\mathbf{z})}{{\tilde{g}}_{IS}(\ell |\mathbf{z})},$$

(21)

where for fixed $\mathbf{z}$,

$${\tilde{g}}_x(\ell |\mathbf{z})=\frac{{\prod }_{i=1}^{N}g\left({\ell }_i\right|\mathbf{z})}{\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})}·{\mathbf{1}}_{\{\ell \in A_{N,x}\}}$$

is the joint density of N independent variables having marginal density $g(·|\mathbf{z})$, conditioned on their average value exceeding the threshold x, and

$${\tilde{g}}_{IS}(\ell |\mathbf{z})=\prod _{i=1}^N{g}_{IS}\left({\ell }_i\right|\mathbf{z}).$$

is the joint density of N independent variables having marginal density $g_{IS}(·|\mathbf{z})$.

Using Equation (21) it is straightforward to decompose the divergence of $h_{IS}$ from $h_x$ as:

$$D(h_{IS}\left|\right|h_x)=D(f_{IS}\left|\right|f_x)+\mathbb{E}\left[\left.D\left({\tilde{g}}_{IS}(·|\mathbf{Z})\right|\left|{\tilde{g}}_x(·|\mathbf{Z})\right)\right|{\overline{L}}_N>x\right],$$

(22)

where $D\left(\xi \right|\left|\eta \right)$ denotes the divergence of the density $\xi$ from the density $\eta$. The first term in Equation (22) is the divergence of $f_{IS}$ from $f_x$, and is therefore minimised by setting $f_{IS}=f_x$. In other words, the best possible IS density for the systematic risk factors (according to the criteria of minimum divergence) is the conditional density $f_x$. The second term in Equation (22) is the average divergence of ${\tilde{g}}_{IS}(·|\mathbf{z})$ from ${\tilde{g}}_x(·|\mathbf{z})$, averaged over all possible realisations of the systematic risk factors and conditioned on portfolio loss exceeding the threshold. Based on the developments in Appendix A.5, for fixed $\mathbf{z}\in {\mathbb{R}}^d$ the divergence of ${\tilde{g}}_{IS}(·|\mathbf{z})$ from ${\tilde{g}}_x(·|\mathbf{z})$ is minimised by setting $g_{IS}(·|\mathbf{z})=g_x(·|\mathbf{z})$. The average divergence in Equation (22) is, therefore, also minimised by setting $g_{IS}(·|\mathbf{z})=g_x(·|\mathbf{z})$ for every $\mathbf{z}\in {\mathbb{R}}^d$.

Remark 4.

Among all densities of the form in Equation (18), the one that most resembles the ideal density $h_x$ (in the sense of minimum divergence) is the density:

$${\widehat{h}}_x(\mathbf{z},\ell ):=f_x\left(\mathbf{z}\right)·\prod _{i=1}^N{g}_x\left({\ell }_i\right|\mathbf{z}),\mathbf{z}\in {\mathbb{R}}^d,\ell \in {[0,{\ell }_{max}]}^N.$$

In other words, ${\widehat{h}}_x$ is the best possible IS density (among the class Equation (18) and according to the criteria of minimum divergence) from which to simulate $(\mathbf{Z},\mathbf{L})$.

It is worth noting that the IS density ${\widehat{h}}_x$ “gets marginal behaviour correct”, in the sense that the marginal distribution of the systematic risk factors, as well as the marginal distribution of an individual loss, is the same under ${\widehat{h}}_x$ as it is under the ideal density $h_x$. The dependence structure of individual losses is different under ${\widehat{h}}_x$ and $h_x$—this is the price that we must pay for insisting on conditional independence (i.e., computational efficiency).

3.3. Approximating the Ideal IS Densities

Simulating directly from ${\widehat{h}}_x$ requires an ability to simulate directly from $f_x$ and $g_x$. Unfortunately, neither $f_x$ nor $g_x$ is numerically tractable (witness the unknown quantities in Equations (15) and (16)), and it does not appear that either is amenable to direct simulation. Our next task is to identify tractable densities that resemble $f_x$ and $g_x$.

3.3.1. Systematic Risk Factors

As a tractable approximation to $f_x$, we suggest using that member of the parametric family in Equation (7) that most resembles $f_x$ in the sense of minimum divergence. Using Equations (7) and (15) we get that:

$$log\left(\frac{f_x\left(\mathbf{z}\right)}{f_{\lambda }\left(\mathbf{z}\right)}\right)=-{\lambda }^{T}S\left(\mathbf{z}\right)+K\left(\lambda \right)+log\left(\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\right)-log\left(p_x\right),$$

whence the divergence of $f_{\lambda }$ from $f_x$ is:

$$D(f_{\lambda }\left|\right|f_x)=-{\lambda }^T\mathbb{E}\left[S\left(\mathbf{Z}\right)\right|{\overline{L}}_N>x]+K\left(\lambda \right)+\mathbb{E}[log\left(\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\right)|{\overline{L}}_N>x]-log\left(p_x\right).$$

(23)

As a cgf, $K(·)$ is strictly convex. As such, Equation (23) attains its unique minimum at that value of $\lambda$ such that:

$$\nabla K\left(\lambda \right)=\mathbb{E}\left[S\left(\mathbf{Z}\right)\right|{\overline{L}}_N>x],$$

(24)

which, in light of Equation (9), is equivalent to:

$${\mathbb{E}}_{\lambda }\left[S\left(\mathbf{Z}\right)\right]=\mathbb{E}\left[S\left(\mathbf{Z}\right)\right|{\overline{L}}_N>x].$$

(25)

Intuitively, we suggest using that value of the IS parameter $\lambda$ for which the mean of $S\left(\mathbf{Z}\right)$ under the IS density matches the conditional mean of $S\left(\mathbf{Z}\right)$, given that portfolio losses exceed the threshold. In what follows we let ${\widehat{\lambda }}_x$ denote that suggested value of the IS parameter $\lambda$, i.e., that value of $\lambda$ that solves Equation (24).

Remark 5.

The first-stage IS weight associated with the so-described density is:

$${\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)=exp(-{\widehat{\lambda }}_x^{T}S\left(\mathbf{Z}\right)+K\left({\widehat{\lambda }}_x\right)).$$

(26)

It is entirely possible—and quite common in the examples we consider in this paper—that $K(-{\widehat{\lambda }}_x)$ is not well-defined, in which case Equation (26) has infinite variance under $f_{{\widehat{\lambda }}_x}$ (recall Remark 1). At first glance it might seem absurd to consider IS densities whose associated weights have infinite variance, but as we discuss in Section 4.2 it is straightforward to circumvent this issue by trimming large first-stage IS weights5.

It remains to develop a tractable approximation to the right hand side of Equation (24), so that we can approximate the value of ${\widehat{\lambda }}_x$. To this end we write the natural sufficient statistic as $S\left(\mathbf{z}\right)=(S_1\left(\mathbf{z}\right),\dots ,S_p\left(\mathbf{z}\right))$ and note that:

$$\mathbb{E}\left[S_i\left(\mathbf{Z}\right)\right|{\overline{L}}_N>x]=\frac{\mathbb{E}[S_i\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}]}{\mathbb{P}({\overline{L}}_N>x)}=\frac{\mathbb{E}[S_i\left(\mathbf{Z}\right)·\mathbb{P}({\overline{L}}_N>x|\mathbf{Z})]}{\mathbb{E}\left[\mathbb{P}({\overline{L}}_N>x|\mathbf{Z})\right]}.$$

Next, we use the LDA in Equation (13) to get:

$$\mathbb{E}\left[S_i\left(\mathbf{Z}\right)\right|{\overline{L}}_N>x]\approx \frac{\mathbb{E}[S_i\left(\mathbf{Z}\right)·exp(-Nq(x,\mathbf{Z}))]}{\mathbb{E}[exp(-Nq(x,\mathbf{Z})\left)\right]}.$$

(27)

As it only involves the systematic risk factors (and not the large number of idiosyncratic risk factors), the expectation on the right hand side of Equation (27) is amenable to either quadrature or Monte Carlo simulation.

3.3.2. Individual Losses

We encourage the reader unfamiliar with exponential tilts to consult Appendix A.3, before reading the remainder of this section. Our approximation to $g_x(\ell |\mathbf{z})$ is obtained by using the LDA of Equation (13) to approximate both conditional probabilities appearing in Equation (16) (see Appendix A.4 for details). The resulting approximation is:

$${\widehat{g}}_x(\ell |\mathbf{z}):=exp(\widehat{\theta }\ell -k(\widehat{\theta },\mathbf{z}))·g(\ell |\mathbf{z}),$$

(28)

where we recall that $\widehat{\theta }$ is defined and discussed in Section 2.3. If the realised values of the systematic risk factors obtained in the first stage lie in the region of interest then $\widehat{\theta }=0$ and ${\widehat{g}}_x$ is identical to g. Otherwise, $\widehat{\theta }$ is strictly positive and ${\widehat{g}}_x$ is the exponentially tilted version of g whose mean is x. Intuitively, we can interpret ${\widehat{g}}_x$ as that density that most resembles (in the sense of minimum divergence) $g_x$, among all densities whose mean is at least x, and the numerical value of $\widehat{\theta }$ as the degree to which the density $g(·|\mathbf{z})$ must be deformed, in order to produce a density whose mean is at least x.

Remark 6.

The mean of Equation (28) is $max\left(\mu \right(\mathbf{z}),x)$. The implication is that the event of interest is not a rare event under the proposed IS algorithm. Indeed,

$${\mathbb{E}}_{IS}\left[L_i\right]={\mathbb{E}}_{IS}\left[{\mathbb{E}}_{IS}\left[L_i\right|\mathbf{Z}]\right]={\mathbb{E}}_{f_{\widehat{\lambda }}}\left[{\mathbb{E}}_{{\widehat{g}}_x}\left[L_i\right|\mathbf{Z}]\right]={\mathbb{E}}_{f_{\widehat{\lambda }}}[max(x,\mu \left(\mathbf{Z}\right))]\ge x,$$

which implies that ${lim}_{N\to \infty }{\mathbb{P}}_{IS}({\overline{L}}_N>x)=1$.

The second-stage IS weight associated with Equation (28) is:

$${\mathsf{\Lambda }}_2(\mathbf{Z},\mathbf{L})=\prod _{i=1}^{N}exp(-\widehat{\theta }{L}_i+k(\widehat{\theta },\mathbf{z}))=exp(-N[\widehat{\theta }{\overline{L}}_N-k(\widehat{\theta },\mathbf{Z})]).$$

Since the second stage weight depends only on $\mathbf{Z}$ and ${\overline{L}}_N$, we will often write ${\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)$ instead of ${\mathsf{\Lambda }}_2(\mathbf{Z},\mathbf{L})$. In order to assess the stability of the second-stage IS weight, we note that:

$$exp(-N[\widehat{\theta }{\overline{L}}_N-k(\widehat{\theta },\mathbf{Z})])=exp(-\widehat{\theta }N[{\overline{L}}_N-x])·exp(-Nq(x,\mathbf{Z})).$$

If $\mathbf{Z}$ lies in the region of interest then $\widehat{\theta }=q=0$, whence ${\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)=1$ whatever the value of ${\overline{L}}_N$. Otherwise, both $\widehat{\theta }$ and q are strictly positive, which implies that ${\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)<1$ whenever ${\overline{L}}_N>x$. The net result of this discussion is that:

$${\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)\le 1\mathrm{whenever}{\overline{L}}_N>x.$$

(29)

The implication is that large, unstable, IS weights in the second stage will never be a problem.

If the realised value of $\mathbf{z}$ does lie in the region of interest then ${\widehat{g}}_x$ and g are identical, and simulation from g is straightforward. Our final task is to determine how to sample from Equation (28) in the case where $\mathbf{z}$ does not lie in the region of interest. One approach would be to identify a family of densities $\{{\eta }_{\mathbf{z}}\left(\mathbf{y}\right):\mathbf{z}\in {\mathbb{R}}^d\}$ such that $L_i=\mathcal{L}(\mathbf{z},{\mathbf{Y}}_i)$ is a draw from ${\widehat{g}}_x(·|\mathbf{z})$ whenever ${\mathbf{Y}}_i$ is a draw from ${\eta }_{\mathbf{z}}(·)$, but this approach appears to be overly complicated. A simpler approach is to sample from Equation (28) using rejection sampling with g as the proposal density. To this end, we note that for fixed $\mathbf{z}$, the ratio of ${\widehat{g}}_x$ to g is $exp(\widehat{\theta }\ell -k(\widehat{\theta },\mathbf{z}))$, which is bounded and strictly increasing on $[0,{\ell }_{max}]$. The best possible (i.e., smallest) rejection constant is therefore:

$$\widehat{c}=\widehat{c}(x,\mathbf{z}):=exp(\widehat{\theta }{\ell }_{max}-k(\widehat{\theta },\mathbf{z})),$$

(30)

and the algorithm for sampling from ${\widehat{g}}_x$ would proceed as follows. First, sample ${\mathbf{Y}}_i$ from its actual density and set ${\widehat{L}}_i=\mathcal{L}(\mathbf{z},{\mathbf{Y}}_i)$. Then generate a random number U, uniformly distributed on $[0,1]$ and independent of ${\mathbf{Y}}_i$. If,

$$U\le \frac{{\widehat{g}}_x\left({\widehat{L}}_i\right|\mathbf{z})}{\widehat{c}·g\left({\widehat{L}}_i\right|\mathbf{z})}=exp(-\widehat{\theta }({\ell }_{max}-{\widehat{L}}_i)),$$

set $L_i={\widehat{L}}_i$ and proceed to the next exposure. Otherwise return to the first step and sample another pair $({\mathbf{Y}}_i,U)$.

3.4. Summary and Intuition

The proposed algorithm is summarised in Algorithm 3 below. The initial step is to approximate the value of the first-stage IS parameter, ${\widehat{\lambda }}_x$. In our numerical examples we use a small pilot simulation (10% of the sample size that we eventually use to estimate $p_x$) and the approximation of Equation (27) in order to estimate ${\widehat{\lambda }}_x$.

Having computed ${\widehat{\lambda }}_x$, the first stage of the algorithm proceeds by simulating independent realisations of the systematic risk factors from the density $f_{{\widehat{\lambda }}_x}$, and computing the associated first-stage weights of Equation (26). Recall that we can interpret these realisations as corresponding to different economic scenarios. Intuitively, sampling from $f_{{\widehat{\lambda }}_x}$ instead of f increases the proportion of adverse scenarios that are generated in the first stage. In the examples we consider, $f_{{\widehat{\lambda }}_x}$ concentrates most of its mass near the boundary of the region of interest, and the effect is to concentrate the distribution of $\mu \left(\mathbf{Z}\right)$ near x.

In the second stage, one first checks whether or not the realised values of the systematic risk factors lie inside the region of interest. If they do then the event of interest is no longer rare and there is no need to apply further IS in the second stage. Otherwise, if we “miss” the region of interest in the first stage, we “correct” this mistake by applying an exponential tilt to the conditional distribution of individual losses. Specifically, we transfer mass from the left tail of g to the right tail, in order to ba x.

Algorithm 3 Proposed IS Algorithm for Estimating $p_x$

1: Compute $\widehat{\lambda }$ using a small pilot simulation.   2: Simulate M i.i.d. copies of the systematic risk factors from $f_{\widehat{\lambda }}\left(\mathbf{z}\right)$ and compute the corresponding first-stage IS weights. Denote the realised values of the factors by ${\mathbf{z}}_1,\dots ,{\mathbf{z}}_M$ and the associated IS weights by ${\mathsf{\Lambda }}_1\left({\mathbf{z}}_1\right),\dots ,{\mathsf{\Lambda }}_1\left({\mathbf{z}}_M\right)$.   3: For each scenario m, determine whether or not ${\mathbf{z}}_m$ lies in the region of interest (i.e., whether or not $\mu \left({\mathbf{z}}_m\right)\ge x$). If it does lie in the region, proceed as follows: (a) Simulate the idiosyncratic risk factors for each exposure. Denote the simulated values by ${\mathbf{y}}_{1,m},\dots ,{\mathbf{y}}_{N,m}$.   (b) Set ${\ell }_{i,m}=\mathcal{L}({\mathbf{z}}_m,{\mathbf{y}}_{i,m})$, ${\overline{\ell }}_m=\frac{1}{N}{\sum }_{i=1}^N{\ell }_{i,m}$ and ${\mathsf{\Lambda }}_2({\mathbf{z}}_m,{\overline{\ell }}_m)=1$.     Otherwise, proceed as follows: (a) Compute $\widehat{\theta }=\widehat{\theta }(x,{\mathbf{z}}_m)$, $\widehat{k}=k(\widehat{\theta },{\mathbf{z}}_m)$ and $\widehat{c}=exp(\widehat{\theta }{\ell }_{max}-\widehat{k})$. For each exposure i: (i) Simulate the exposure’s idiosyncratic risk factor (denote the realised value by ${\widehat{\mathbf{y}}}_{i,m}$) and set ${\widehat{\ell }}_{i,m}=\mathcal{L}({\mathbf{z}}_m,{\mathbf{y}}_{i,m})$.   (ii) Simulate a random number drawn uniformly from the unit interval (denote the realised value by u) and determine whether or not $u\le exp(-\widehat{\theta }({\ell }_{max}-{\widehat{\ell }}_{i,m}))$. If it is, set ${\ell }_{i,m}={\widehat{\ell }}_{i,m}$ and proceed to the next exposure. Otherwise, return to step (i). (b) Set ${\overline{\ell }}_m=\frac{1}{N}{\sum }_{i=1}^N{\ell }_{i,m}$ and ${\mathsf{\Lambda }}_2({\mathbf{z}}_m,{\overline{\ell }}_m)=exp(-N[\widehat{\theta }{\overline{\ell }}_m-\widehat{k}])$ 4: Return ${\widehat{p}}_x=\frac{1}{M}{\sum }_{m=1}^M{\mathsf{\Lambda }}_1\left({\mathbf{z}}_m\right)·{\mathsf{\Lambda }}_2({\mathbf{z}}_m,{\overline{\ell }}_m)·{\mathbf{1}}_{\{{\overline{\ell }}_m>x\}}$.

4. Practical Considerations

In this section we discuss some of the practical issues that arise when implementing the proposed methodology.

4.1. One- and Two-Stage Estimators

The rejection sampling procedure employed in the second stage of the proposed algorithm involves repeated evaluation of $\widehat{\theta }$, which requires a non-trivial amount of computational time time. In addition, rejection sampling in general requires relatively complicated code. As such, it is worth considering a simpler algorithm that only applies importance sampling in the first stage, and is therefore easier to implement and faster to run.

In what follows we will distinguish between one- and two-stage IS algorithms. A one-stage algorithm only applies IS in the first stage and samples $(\mathbf{Z},\mathbf{L})$ from the IS density:

$$h_{1\mathrm{S}}(\mathbf{z},\ell ):=f_{{\widehat{\lambda }}_x}\left(\mathbf{z}\right)·\prod _{i=1}^{N}g\left({\ell }_i\right|\mathbf{z}).$$

(31)

The associated IS weight is ${\mathsf{\Lambda }}_1\left(\mathbf{z}\right)$ and the one-stage algorithm is summarised in Algorithm 4 below. Note the simplicity of Algorithm 4, relative to Algorithm 3. The two-stage algorithm applies IS in both the first stage and the second stage, sampling ($\mathbf{Z},\mathbf{L})$ from the IS density:

$$h_{2\mathrm{S}}(\mathbf{z},\ell ):=f_{{\widehat{\lambda }}_x}\left(\mathbf{z}\right)·\prod _{i=1}^N{\widehat{g}}_x\left({\ell }_i\right|\mathbf{z}).$$

(32)

The associated IS weight is ${\mathsf{\Lambda }}_1\left(\mathbf{z}\right)·{\mathsf{\Lambda }}_2(\mathbf{z},{\overline{\ell }}_N)$, and the two-stage algorithm was summarised previously in Algorithm 3.

Algorithm 4 Proposed One-Stage IS Algorithm for Estimating $p_x$

1: Compute ${\widehat{\lambda }}_x$ using a small pilot simulation.   2: Simulate M i.i.d. copies of the systematic risk factors from $f_{\widehat{\lambda }}\left(\mathbf{z}\right)$ and compute the corresponding first-stage IS weights. Denote the realised values of the factors by ${\mathbf{z}}_1,\dots ,{\mathbf{z}}_M$ and the associated IS weights by ${\mathsf{\Lambda }}_1\left({\mathbf{z}}_1\right),\dots ,{\mathsf{\Lambda }}_1\left({\mathbf{z}}_M\right)$.   3: For each scenario m: (a) Simulate the idiosyncratic risk factors for each exposure. Denote the simulated values by ${\mathbf{y}}_{1,m},\dots ,{\mathbf{y}}_{N,m}$.   (b) Set ${\ell }_{i,m}=\mathcal{L}({\mathbf{z}}_m,{\mathbf{y}}_{i,m})$ and ${\overline{\ell }}_m=\frac{1}{N}{\sum }_{i=1}^N{\ell }_{i,m}$. 4: Return ${\widehat{p}}_x=\frac{1}{M}{\sum }_{m=1}^M{\mathsf{\Lambda }}_1\left({\mathbf{z}}_m\right)·{\mathbf{1}}_{\{{\overline{\ell }}_m>x\}}$.

Although it is simpler to implement and faster to run, the one-stage algorithm is less accurate than the two-stage algorithm. More precisely, the two-stage estimator never has larger variance than the one-stage estimator. To see this, first let ${\mathbb{E}}_{1\mathrm{S}}$ denote expectation under the one-stage IS density $h_{1\mathrm{S}}(\mathbf{z},\ell )$ given in Equation (31). Then the variance of the one-stage estimator is:

$${\displaystyle \frac{{\mathbb{E}}_{1\mathrm{S}}[{({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}})}^2]-p_x^2}{M}},$$

where M denotes sample size. And if we let ${\mathbb{E}}_{2\mathrm{S}}$ denote expectation under the two-stage IS density $h_{2\mathrm{S}}(\mathbf{z},\ell )$ given in Equation (32) then the variance of the two-stage estimator is:

$${\displaystyle \frac{{\mathbb{E}}_{2\mathrm{S}}[{({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}})}^2]-p_x^2}{M}}.$$

In order to compare variances it suffices to compare the second moments appearing above under the actual density $h(\mathbf{z},\ell )$, and we let $\mathbb{E}$ denote expectation with respect to this density. To this end we note that:

$${\mathbb{E}}_{1\mathrm{S}}[{({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}})}^2]=\mathbb{E}[{\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}}]$$

and

$${\mathbb{E}}_{2\mathrm{S}}[{({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}})}^2]=\mathbb{E}[{\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)]·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}}]$$

In light of Equation (29) we get that:

$${\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}\le 1·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}={\mathbf{1}}_{\{{\overline{L}}_N>x\}},$$

(33)

whence

$$\begin{array}{lll}{\mathbb{E}}_{2\mathrm{S}}[{({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}})}^2]\hfill & =& \mathbb{E}[{\mathsf{\Lambda }}_1(\mathbf{Z})·{\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}}]\hfill \\ & \le & \mathbb{E}[{\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}}]\hfill \\ & =& {\mathbb{E}}_{1\mathrm{S}}[{({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N\ge x\}})}^2].\hfill \end{array}$$

The two-stage estimator will therefore never have larger variance than the the one-stage estimator.

4.2. Large First-Stage Weights

In the examples that we consider in this paper, the systematic risk factors are Gaussian. When selecting their IS density, one could either (i) shift their means and leave their variances (and correlations) unchanged or (ii) shift their means and adjust their variances (and correlations). In general the latter approach will lead to a much better approximation to the ideal density $f_x$, but could lead to an IS weight that has infinite variance. By contrast, the former approach will always lead to an IS weight with finite variance, but could lead to a poor approximation of the ideal density. At first glance it might seem absurd to consider IS densities whose weights are so unstable as to have infinite variance, but we have found that adjusting the variances of the systematic risk factors can lead to more effective estimators, in terms of both statistical accuracy and run time (see Section 6.1 for more details), provided one stabilises the resulting IS weights in some way. In the remainder of this section we describe a simple stabilisation technique that leads to a computable upper bound on the associated bias (an alternative would be to stabilize unruly IS weights via truncation, as discussed in Ionides (2008)).

Returning now to the general case, suppose that the first-stage IS parameter, ${\widehat{\lambda }}_x$, is such that the first-stage IS weight, ${\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)$, has infinite variance. We trim large first-stage weights by fixing a set $A\subset {\mathbb{R}}^d$ such that ${\mathsf{\Lambda }}_1(·)$ is bounded over A, and discarding those simulations for which $\mathbf{Z}\notin A$. Specifically, the last line of Algorithm 3 would be altered to return the trimmed estimate:

$${\widehat{p}}_x=\frac{1}{M}\sum _{m=1}^M{\mathsf{\Lambda }}_1\left({\mathbf{z}}_m\right)·{\mathsf{\Lambda }}_2(\mathbf{z},{\overline{\ell }}_m)·{\mathbf{1}}_{\{{\overline{\ell }}_m>x\}}·{\mathbf{1}}_{\{{\mathbf{z}}_m\in A\}},$$

etc. for Algorithm 4. The variance of the so-trimmed estimator is necessarily finite (recall that ${\mathsf{\Lambda }}_2(\mathbf{z},\overline{\ell })\le 1$ if $\overline{\ell }>x$), and its bias is:

$${\mathbb{E}}_{2\mathrm{S}}[{\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathsf{\Lambda }}_2(\mathbf{Z},{\overline{L}}_N)·{\mathbf{1}}_{\{{\overline{L}}_N>x\}}·{\mathbf{1}}_{\{\mathbf{Z}\notin A\}}]=\mathbb{E}[{\mathbf{1}}_{\{{\overline{L}}_N>x\}}·{\mathbf{1}}_{\{\mathbf{Z}\notin A\}}]=\mathbb{E}[\mathbb{P}({\overline{L}}_N>x|\mathbf{Z})·{\mathbf{1}}_{\{\mathbf{Z}\notin A\}}],$$

where we have used the tower property (conditioning on $\mathbf{Z}$) to obtain the last equality. Using Chernoff’s bound in Equation (11) we get that:

$$\mathbb{E}[\mathbb{P}({\overline{L}}_N>x|\mathbf{Z})·{\mathbf{1}}_{\{\mathbf{Z}\notin A\}}]\le \mathbb{E}[exp(-Nq(x,\mathbf{Z}))·{\mathbf{1}}_{\{\mathbf{Z}\notin A\}}].$$

(34)

As it only depends on the small number of systematic risk factors, and not the large number of idiosyncratic risk factors, the right-hand side of Equation (34) is a tractable upper bound on the bias committed by trimming large (first-stage) IS weights. This upper bound can be used to assess whether or not the bias associated with a given set A is acceptable.

4.3. Large Rejection Constants

The smaller the $\widehat{c}$, the more efficient is the rejection sampling algorithm employed in the second stage. Indeed the average number of proposals that must be generated in order to obtain one realisation from ${\widehat{g}}_x$ is $1/\widehat{c}$. In the examples we consider in this paper, $\widehat{c}$ is (essentially) a decreasing function $\mu \left(\mathbf{z}\right)$, such that $\widehat{c}\to 1$ as $\mu \left(\mathbf{z}\right)\to x$ and $\widehat{c}\to \infty$ as $\mu \left(\mathbf{z}\right)\to 0$ (see Figure 1). The second-stage rejection algorithm is therefore quite efficient when $\mu \left(\mathbf{z}\right)\approx x$ and quite inefficient when $\mu \left(\mathbf{z}\right)\approx 0$. Now, the IS density for the first-stage risk factors is such that the distribution of $\mu \left(\mathbf{Z}\right)$ concentrates most of its mass near x (where $\widehat{c}$ is a reasonable size), but it is still theoretically possible to obtain a realisation of the systematic risk factors for which $\mu \left(\mathbf{z}\right)$ is very small and $\widehat{c}$ is unacceptably large (e.g., ${10}^4$). In such situations the algorithm effectively grinds to a halt, as one endlessly generates proposed losses that have no realistic chance of being accepted. It is extremely unlikely that one does obtain such a scenario under the first-stage IS distribution, but it is still important to protect oneself against this unlikely event. To this end we suggest fixing some maximum acceptable rejection constant $c_{max}$, and only applying the second stage IS to those first-stage realizations for which $\mu \left(\mathbf{z}\right)<x$ and $\widehat{c}\le c_{max}$. In other words, even if the realised values of the systematic risk factors lie outside the region of interest, we avoid applying the second stage if the associated rejection constant exceeds the predefined threshold.

4.4. Computing $\widehat{\theta }$

Repeated evaluations of $\widehat{\theta }(x,·)$ are necessary when computing ${\widehat{\lambda }}_x$ at the outset of the algorithm, as well as during the second stage of the two-stage algorithm. Recall that in order to compute $\widehat{\theta }(x,\mathbf{z})$ “exactly” one must numerically solve the equation $k^{\prime }(\theta ,\mathbf{z})=x$, which requires a non-trivial amount of CPU time. As each evaluation of $\widehat{\theta }$ is relatively costly, repeated evaluation would, in the absence of any further approximation (over and above that inherent in numerical root-finding), account for the vast majority of the algorithm’s total run time.

In order to reduce the amount of time spent evaluating $\widehat{\theta }$ we fit a low degree polynomial to the function $\widehat{\theta }(x,·)$ that can be evaluated extremely quickly, considerably reducing total run time. Specifically, suppose that we must compute $\widehat{\theta }(x,{\mathbf{z}}_n)$ for each of n points ${\mathbf{z}}_1,\dots ,{\mathbf{z}}_n$ (either the sample points from the pilot simulation, or the first-stage realisations that did not land in the region of interest). We identify a small set $C\subset {\mathbb{R}}^d$ that contains each of the n points, construct a mesh of $m<<n$ points in C, evaluate $\widehat{\theta }$ exactly at each mesh point, and then fit a fifth degree polynomial to the resulting data. Letting $\overline{\theta }(x,·)$ denote the resulting polynomial, we then evaluate $\overline{\theta }(x,{\mathbf{z}}_1),\dots ,\overline{\theta }(x,{\mathbf{z}}_n)$ instead of $\widehat{\theta }(x,{\mathbf{z}}_1),\dots ,\widehat{\theta }(x,{\mathbf{z}}_n)$. If m is substantially smaller than n, then the reduction in CPU time is considerable.

5. PD-LGD Correlation Framework

All of the PD-LGD correlation models listed in the introduction are special cases of the following general framework—an observation that, to the best of our knowledge, has not been made in the literature. The systematic risk factors take the form $\mathbf{Z}=(Z_D,Z_L)$, where $Z_D$ and $Z_L$ are bivariate normal with standard normal margins and correlation ${\rho }_S$. Idiosyncratic risk factors take the form ${\mathbf{Y}}_i=(Y_{i,D},Y_{i,L})$, where $Y_{i,D}$ and $Y_{i,L}$ are bivariate normal with standard normal margins and correlation ${\rho }_I$.

Associated with each exposure is a default driver $X_{i,D}$ and a loss driver $X_{i,L}$, defined as follows:

$$\begin{array}{c}\hfill X_{i,D}={\alpha }_D{Z}_D+\sqrt{1-{\alpha }_D^2}{Y}_{i,D},\end{array}$$

(35)

$$\begin{array}{c}\hfill X_{i,L}={\alpha }_L{Z}_L+\sqrt{1-{\alpha }_L^2}{Y}_{i,L}.\end{array}$$

(36)

The factor loadings ${\alpha }_D$ and ${\alpha }_L$ are constants taking values in the unit interval, and dictate the relative importance of systematic risk versus idiosyncratic risk. The correlation between default drivers of distinct exposures is ${\rho }_D:={\alpha }_D^2$ and the correlation between loss drivers of distinct exposures is ${\rho }_L:={\alpha }_L^2$. The correlation between the default and potential loss drivers of a particular exposure is:

$${\rho }_{DL}:={\alpha }_D{\alpha }_L{\rho }_S+\sqrt{1-{\alpha }_D^2}\sqrt{1-{\alpha }_L^2}{\rho }_I,$$

which can be positive or negative (or zero). Note that if ${\rho }_S$ and ${\rho }_I$ have the same sign then, since both factor loadings are positive, ${\rho }_{DL}$ inherits this common sign.

The realised loss on exposure i is $L_i={\mathcal{D}}_i·{\mathcal{L}}_i$, where:

$${\mathcal{D}}_i={\mathbf{1}}_{\{X_{i,D}\le {\mathsf{\Phi }}^{-1}\left(P\right)\}}$$

is the default indicator associated with exposure i and

$${\mathcal{L}}_i=h\left(X_{i,L}\right)$$

is called the potential loss (our terminology) associated with exposure i. Here P denotes the common default probability of all exposures and h is some function from $\mathbb{R}$ to $[0,{\ell }_{max}]$. It is useful (but not necessary) to think of potential loss as ${\mathcal{L}}_i=max(0,1-{\mathcal{C}}_i)$, where ${\mathcal{C}}_i$ is the value of the collateral pledged to exposure i expressed as a fraction of the loan’s notional value.

Models in this framework are characterised by (i) the correlation structure of the risk factors, specifically restrictions on the values of ${\rho }_I$ and ${\rho }_S$, and (ii) the marginal distribution of potential loss. For instance:

• Frye (2000) assumes perfect systematic correlation (${\rho }_S=1$) and zero idiosyncratic correlation (${\rho }_I=0$);
• Pykhtin (2003) assumes perfect systematic correlation (${\rho }_S=1$) but allows for arbitrary idiosyncratic correlation (${\rho }_I$ unrestricted);
• Witzany (2011) allows for arbitrary systematic correlation (${\rho }_S$ unrestricted) but insists on zero idiosyncratic correlation (${\rho }_I=0$);
• Miu and Ozdemir (2006) allow for arbitrary systematic correlation (${\rho }_S$ unrestricted) and arbitrary idiosyncratic correlation (${\rho }_I$ unrestricted).

Note that if ${\rho }_S=\pm 1$ then the systematic risk factor is effectively one-dimensional. Indeed if ${\rho }_S=1$ then $\mathbf{Z}=(Z,Z)$ from some standard Gaussian variable Z, and if ${\rho }_S=-1$ then $\mathbf{Z}=(Z,-Z)$. We refer to the case $|{\rho }_S|=1$ as the one-factor case, and the case $|{\rho }_S|<1$ as the two-factor case. In the one-factor case we use Z, and not $\mathbf{Z}$, to denote the systematic risk factor. The first two models listed above are one-factor models, the last two are two-factor models.

The marginal distribution of potential loss is determined by the specification of the function h. For instance:

• Frye (2000) specifies $h\left(x\right)=max(0,1-a(1+bx\left)\right)$ for constants $a\in \mathbb{R}$ and $b>0$. Potential loss takes values in $[0,\infty )$. Its density has a point mass at zero and is proportional to a Gaussian density on $(0,\infty )$. Since $L_i$ is not constrained to lie in the unit interval, this specification violates the assumptions made in Section 2.3;
• Pykhtin (2003) specifies $h\left(x\right)=max(0,1-e^{a+bx})$ for constants $a\in \mathbb{R}$ and $b>0$. Potential loss takes values in $[0,1)$. Its density has a point mass at zero, and is proportional to a shifted lognormal density over $(0,1)$;
• Witzany (2011) and Miu and Ozdemir (2006) both specify $h\left(x\right)=B_{a,b}^{-1}\left(\mathsf{\Phi }\left(x\right)\right)$, where $a,b>0$ and $B_{a,b}$ denotes the cdf of the beta distribution with parameters a and b. Potential loss takes values in $(0,1)$. It is a continuous variable and follows a beta distribution.

The sign of ${\rho }_{DL}$ and the nature of the function h (increasing or decreasing) will in general determine the sign of the relationship between ${\mathcal{D}}_i$ and ${\mathcal{L}}_i$. If ${\rho }_{DL}>0$ then the relationship will be positive [negative] provided h is decreasing [increasing], and vice versa if ${\rho }_{DL}<0$.

5.1. Computing $\mu \left(\mathbf{z}\right)$

Here vectors $\mathbf{z}\in {\mathbb{R}}^2$ take the form $\mathbf{z}={(z_D,z_L)}^T$. In order to obtain an expression for $\mu \left(\mathbf{z}\right)=\mathbb{E}\left[L_i\right|\mathbf{Z}=\mathbf{z}]$, we begin with the observation that:

$$\mathbb{E}\left[L_i\right|\mathbf{Z}]=\mathbb{E}\left[{\mathcal{L}}_i{\mathcal{D}}_i\right|\mathbf{Z}]=\mathbb{E}\left[{\mathcal{L}}_i\mathbb{E}\left[{\mathcal{D}}_i\right|{\mathbf{X}}_{i,L},\mathbf{Z}]\right|\mathbf{Z}]=\mathbb{E}\left[{\mathcal{L}}_i\mathbb{P}({\mathcal{D}}_i=1|X_{i,L},\mathbf{Z})\right|\mathbf{Z}].$$

Thus,

$$\mu \left(\mathbf{z}\right)={\int }_{\mathbb{R}}h\left(x_L\right)·\mathsf{\Phi }(d,m(x_L,\mathbf{z}),v)·\varphi (x_L,{\alpha }_L{z}_L,1-{\alpha }_L^2)d{x}_L,$$

(37)

where

$$m(x_L,\mathbf{z}):={\alpha }_D{z}_D+{\rho }_I·\sqrt{\frac{1-{\alpha }_D^2}{1-{\alpha }_L^2}}·(x_L-{\alpha }_L{z}_L)$$

and

$$v=v(x_L,\mathbf{z}):=(1-{\alpha }_D^2)(1-{\rho }_I^2)$$

are the conditional mean and variance of $X_{i,D}$, respectively, given that $(X_{i,L},\mathbf{Z})=(x_L,\mathbf{z})$. In general $\mu \left(\mathbf{z}\right)$ must be evaluated using quadrature, and doing so is straightforward6. On average (across parameter values and points $\mathbf{z}\in {\mathbb{R}}^2)$ a single evaluation of $\mu (·)$ requires approximately one millisecond. In the one-factor case with ${\rho }_S=1$ [${\rho }_S=-1$] the expression for $\mu \left(z\right)=\mathbb{E}\left[L_i\right|Z=z]$ is obtained by plugging $\mathbf{z}=(z,z)$ [$\mathbf{z}=(z,-z)$] into Equation (37).

5.2. Computing $k(\theta ,\mathbf{z})$ and $\widehat{\theta }(x,\mathbf{z})$

Here again, vectors $\mathbf{z}\in {\mathbb{R}}^2$ take the form $\mathbf{z}={(z_D,z_L)}^T$. In order to derive an expression for $k(\theta ,\mathbf{z})$ we begin with the observation that:

$$e^{\theta L_i}=\mathbf{1}({\mathcal{D}}_i=0)+e^{\theta {\mathcal{L}}_i}\mathbf{1}({\mathcal{D}}_i>0)=1+(e^{\theta {\mathcal{L}}_i}-1)·\mathbf{1}({\mathcal{D}}_i>0),$$

and since $k(\theta ,\mathbf{z})=log\left(\mathbb{E}\left[e^{\theta L_i}\right|\mathbf{Z}=\mathbf{z}]\right)$, we get that:

$$k(\theta ,\mathbf{z})=log\left(1+{\int }_{\mathbb{R}}(e^{\theta h\left(x_L\right)}-1)·\mathsf{\Phi }(d,m(x_L,\mathbf{z}),v)·\varphi (x_L,{\alpha }_L{z}_L,1-{\alpha }_L^2)d{x}_L\right),$$

(38)

where $m(x_L,\mathbf{z})$ and v are given in the previous section. In the one-factor case with ${\rho }_S=1$ [${\rho }_S=-1$] the expression for $k(\theta ,z)=log\left(\mathbb{E}[exp\left(\theta L_i\right)|Z=z]\right)$ is obtained by plugging $\mathbf{z}=(z,z)$ [$\mathbf{z}=(z,-z)$] into Equation (38). As with $\mu \left(\mathbf{z}\right)$, $k(\theta ,\mathbf{z})$ must in general be evaluated using quadrature, which is straightforward. The time required for a single evaluation of $k(\theta ,·)$ is comparable to that required for a single evaluation of $\mu (·)$.

In order to compute $\widehat{\theta }$ we must solve the equation $k^{\prime }(\theta ,\mathbf{z})=x$ with respect to $\theta$. Differentiating Equation (38) we get:

$$k^{\prime }(\theta ,\mathbf{z})=\frac{\partial k(\theta ,\mathbf{z})}{\partial \theta }=\frac{{\int }_{\mathbb{R}}h\left(x_L\right)·e^{\theta h\left(x_L\right)}·\mathsf{\Phi }(d,m(x_L,\mathbf{z}),v)·\varphi (x_L,{\alpha }_L{z}_L,1-{\alpha }_L^2)d{x}_L}{exp\left(k\right(\theta ,\mathbf{z}\left)\right)},$$

(39)

which is straightforward to compute using quadrature. A single evaluation of $k^{\prime }(\theta ,\mathbf{z})$ requires approximately twice as much time as a single evaluation of $k(\theta ,\mathbf{z})$. As the root of $k^{\prime }(\theta ,\mathbf{z})=x$ must be evaluated numerically, evaluating $\widehat{\theta }$ is much more time consuming than evaluating k or $k^{\prime }$. Across parameter values and points $\mathbf{z}\in {\mathbb{R}}^2$, and using $\theta =0$ as an initial guess, the average time required for a single evaluation7 of $\widehat{\theta }(x,·)$ is slightly less than one tenth of one second.

The right panel of Figure 1 illustrates the relationship between expected losses and the rejection constant employed in the second stage, $\widehat{c}=exp(\widehat{\theta }-k(\widehat{\theta },\mathbf{z}))$. We see that $\widehat{c}$ is essentially a decreasing function of $\mu \left(\mathbf{z}\right)$, such that $\widehat{c}\to 1$ as $\mu \left(\mathbf{z}\right)\to x$ and $\widehat{c}\to \infty$ as $\mu \left(\mathbf{z}\right)\to 0$. The left panel of Figure 1 illustrates the graph of the LDA approximation $\mathbb{P}({\overline{L}}_N>x|\mathbf{Z}=\mathbf{z})\approx exp(-Nq(x,\mathbf{z}))$. The approximation is identically equal to one inside the region of interest, and decays to zero very rapidly outside the region. In other words, most of the variability in the function $q(x,·)$ occurs along, and just outside, the boundary of the region of interest.

5.3. Exploring the Parameter Space

The model contains five parameters, in addition to any parameters associated with the transformation h. We are ultimately interested in how well the proposed algorithms perform across a wide range of different parameter sets. As such, in our numerical experiments we will randomly select a large number of parameter sets according to the procedure described below, and assess the algorithms’ performance for each parameter set.

• Generate the default probability P uniformly between 0% and 10%, and generate each of the correlations ${\rho }_D={\alpha }_D^2$ and ${\rho }_L={\alpha }_L^2$ uniformly between 0% and 50%;
• In the one-factor model, generate ${\rho }_S$ uniformly on $\{-1,1\}$, i.e., ${\rho }_S$ takes on the value $-1$ or $+1$ with equal probability. If ${\rho }_S=1$ we generate ${\rho }_I$ uniformly between 0% and 100%, and if ${\rho }_S=-1$ we generate ${\rho }_I$ uniformly between $-100\%$ and $0\%$. This allows us to control the sign of ${\rho }_{DL}$, which we must do in order to ensure a positive relationship between default and potential loss. In the two-factor model we randomly generated ${\rho }_S$ uniformly on $[-1,1]$. If ${\rho }_S$ is positive, randomly generate ${\rho }_I$ uniformly on $[0,1]$, otherwise randomly generate ${\rho }_I$ uniformly on $[-1,0]$;
• We choose the transformation $h(·)$ to ensure that (i) potential loss is beta distributed and (ii) there is a positive relationship between default and loss. The paramters a and b of the beta distribution are generated independently from an exponential distribution with unit mean. If ${\rho }_{DL}<0$ we set $h\left(x\right)=B_{a,b}^{-1}\left(\mathsf{\Phi }\left(x\right)\right)$ and if ${\rho }_{DL}>0$ we set $h\left(x\right)=B_{a,b}^{-1}\left(\mathsf{\Phi }(-x)\right)$, where $B_{a,b}(·)$ is the cumulative distribution function for the beta distribution with parameters a and b. Note that under these restrictions, in the one-factor model the expected loss function $\mu \left(z\right)$ is monotone decreasing.

In order to ensure that we are considering cases of practical interest, we randomise the portfolio size and loss threshold as follows.

• Generate the number of exposures randomly between 10 and 5000;
• In the one-factor model we generate the threshold x by setting $x=\mu \left({\mathsf{\Phi }}^{-1}\left({10}^{-q}\right)\right)$, where q is uniformly distributed on $[1,5]$. The LPA suggests that

  $$p_x=\mathbb{P}({\overline{L}}_N>x)\approx \mathbb{P}(\mu \left(Z\right)>x)=\mathbb{P}(Z<{\mu }^{-1}\left(x\right))={10}^{-q}.$$

  This means that $log\left(p_x\right)$, the order of magnitude of the probability of interested, is approximately uniformly distributed on $[-5,-1]$. In the two-factor model we set $x=\mu \left({\mathbf{z}}_q\right)$, where ${\mathbf{z}}_q=({\mathsf{\Phi }}^{-1}\left(q\right),{\rho }_S{\mathsf{\Phi }}^{-1}\left(q\right))$ and q is uniformly distributed on $[-5,-1]$.

6. Implementation

In this section we discuss our implementation of the algorithm proposed in Section 3 in the general framework outlined in Section 5. As the general framework encompasses many of the PD-LGD correlations that have been proposed in the literature, this section effectively discusses implementation of the proposed algorithm across a wide variety of models that are used in practice.

6.1. Selecting the IS Density for the Systematic Risk Factors

The systematic risk factors here are Gaussian. When constructing their IS density we could either shift their means and leave their variances (and correlations) unchanged, or shift their means and adjust their variances (and correlations). Recall that the ultimate goal is to choose an IS density that closely resembles the ideal density $f_x$ given in Equation (15). As illustrated8 in Figure 2, the ideal density $f_x$ tends to be very tightly concentrated about its mean, and adjusting the variance of the systematic risk factors leads to a much better approximation to the ideal density for “typical values” of the ideal density. The left tail of the ideal density is, however, heavier than the variance-adjusted IS density, an issue that can be resolved by trimming large IS weights.

The downside to adjusting the variance of the systematic risk factors is that it can lead to first-stage IS weights with infinite variance, but numerical evidence suggests that this issue can be mitigated by trimming large weights. Indeed, numerical experiments9 suggest that adjusting variance and trimming large weights leads to substantially more accurate estimators of $p_x$. Intuitively, it is more important for the IS density to mimic the behaviour of the ideal density over its “typical range”, as opposed to faithfully representing its tail behaviour. In addition to improving statistical accuracy, adjusting variance has the added benefit of making the second stage of the algorithm more computationally efficient in terms of run time. Indeed, as discussed in more detail in Section 6.3, adjusting variance tends to increase the proportion of first-stage simulations that land in the region of interest (thereby reducing the number of times the rejection sampling algorithm must be employed in the second stage) and reduces the average size of the rejection constants employed in the second stage (thereby making the rejection algorithm more effective whenever it must be employed).

6.2. First Stage

In this section we explain how to efficiently approximate the parameters of the optimal IS density for the systematic risk factors, in both the one- and two-factor models. We also explain how we trim large IS weights, and demonstrate that the resulting bias is negligible.

6.2.1. Computing Parameters in the Two-Factor Model

In the two-factor model the systematic risk factors are bivariate Gaussian with zero mean vector and covariance matrix:

$$\mathsf{\Sigma }=\left[\begin{array}{cc}1& {\rho }_S\\ {\rho }_S& 1\end{array}\right].$$

The mean vector and covariance matrix that satisfy the criteria of Equation (25) are:10

$${\mu }_{IS}:=\mathbb{E}\left[\mathbf{Z}\right|{\overline{L}}_N>x]$$

(41)

and

$${\mathsf{\Sigma }}_{IS}:=\mathbb{E}\left[(\mathbf{Z}-{\mu }_{IS}){(\mathbf{Z}-{\mu }_{IS})}^T\right|{\overline{L}}_N>x],$$

(42)

respectively. In order to approximate the suggested mean vector and covariance matrix we use Equation (27) to get:

$${\mu }_{IS}\approx \frac{\mathbb{E}[exp(-Nq(x,\mathbf{Z}))·\mathbf{Z}]}{\mathbb{E}[exp(-Nq(x,\mathbf{Z})\left)\right]}$$

(43)

and

$${\mathsf{\Sigma }}_{IS}\approx \frac{\mathbb{E}[exp(-Nq(x,\mathbf{Z}))·(\mathbf{Z}-{\mu }_{IS}){(\mathbf{Z}-{\mu }_{IS})}^T]}{\mathbb{E}[exp(-Nq(x,\mathbf{Z})\left)\right]}.$$

(44)

The expected values appearing on the right-hand sides of Equations (43) and (44) are both amenable to simulation, and we use a small pilot simulation of size $M_p<<M$ to approximate them. In our numerical examples, the size of the pilot simulation is 10% of the sample size that is eventually used to estimate $p_x$.

In order to implement the approximation we must first simulate the systematic risk factors and then compute $q(x,\mathbf{z})$ for each sample point $\mathbf{z}$. The most natural way to proceed is to (i) sample the systematic risk factors from their actual distribution (bivariate Gaussian with zero mean vector and covariance matrix $\mathsf{\Sigma }$) and (ii) numerically solve the equation $k^{\prime }(\theta ,\mathbf{z})=x$ in order to compute compute $\widehat{\theta }(x,\mathbf{z})$ for each pilot sample point $\mathbf{z}$ that lies outside the region of interest. In our experience this leads to unacceptably inefficient estimators, in terms of both (i) statistical accuracy and (ii) computational time. We deal with each issue in turn.

As most of the variation in $q(x,·)$ occurs just outside the boundary of the region of interest (recall the right panel of Figure 1), we suggest using an IS distribution for the pilot simulation that is centered on the boundary of the region. Specifically, we suggest using that point on the boundary at which the density of the systematic risk factors attains its maximum value (i.e., the most likely point on the boundary):

$${\mathbf{z}}_x:=argmin\{{\mathbf{z}}^T{\mathsf{\Sigma }}^{-1}\mathbf{z}:\mu \left(\mathbf{z}\right)=x\}.$$

(45)

The non-linear minimisation problem appearing above is easily and rapidly solved using standard techniques. We used fmincon function in Matlab.

As ${\mathbf{z}}_x$ lies on the boundary of the region of interest, roughly half the pilot sample will lie outside the region. In Section 5.2 we noted that it takes nearly one tenth of one second to numerically solve the equation $k^{\prime }(\theta ,\mathbf{z})=x$. As such, if we are to compute $\widehat{\theta }$ exactly (i.e., by numerically solving the indicated equation) for each sample point that lies outside the region of interest, the total time required (in seconds) to estimate the first-stage IS parameters will be at least $M_p/20$. In our numerical examples we use a pilot sample size of $M_p=1000$, which means that it would take nearly one full minute to compute the first-stage IS parameters. This discussion suggests that reducing the number of times we must numerically solve the equation $k^{\prime }(\theta ,\mathbf{z})=x$ could lead to a dramatic reduction in computational time.

We suggest fitting a low degree polynomial to the function $\widehat{\theta }(x,·)$, over a small region in ${\mathbb{R}}^2$ that contains all of the pilot sample points that lie outside the region of interest. Specifically, we determine the smallest rectangle that contains all of the pilot sample points, and discretize the rectangle using a mesh of $n_g^2$ points, equally spaced in each direction. Next, we identify those mesh points that lie outside the region of interest and compute $\widehat{\theta }(x,\mathbf{z})$ exactly (i.e., by solving $k^{\prime }(\theta ,\mathbf{z})=x$ numerically) for each such point. Finally, we fit a polynomial to the resulting $(\mathbf{z},\widehat{\theta }(x,\mathbf{z}))$ pairs and call the resulting function $\overline{\theta }(x,·)$. Numerical evidence indicates at using a fifth-degree polynomial and a mesh with ${15}^2=225$ points leads to a sufficiently accurate approximation to $\widehat{\theta }(x,·)$ over the indicated range (the intersection of (i) the smallest rectangle that contains all sample points and (ii) the complement of the region of interest). Note that $\overline{\theta }$ could be an extremely inaccurate approximation to $\widehat{\theta }$ outside this range, but that is not a concern because we will never need to evaluate it there.

It remains to compute $q(x,\mathbf{z})$ for each of the pilot points $\mathbf{z}$. For those points $\mathbf{z}$ that lie inside the region of interest, we set $q(x,\mathbf{z})=0$. For those points that lie inside the region, we set $q(x,\mathbf{z})=\overline{\theta }x-k(\overline{\theta },\mathbf{z})$, where $\overline{\theta }=\overline{\theta }(x,\mathbf{z})$. Evaluating $\overline{\theta }(x,·)$ requires essentially no computational time (it is a polynomial), and if the mesh size and degree are chosen appropriately the difference between $\widehat{\theta }$ and $\overline{\theta }$ is very small. In total, the suggested procedure reduces the number of evaluations of $\widehat{\theta }$ from $M_p/2$ to $n_g/2$, for a percentage reduction of $n_g^2/M_p$. In our numerical examples we use $n_g=15$ and $M_p=1000$, which corresponds to a reduction of 75% in computational time.

To summarise, we estimate the optimal first-stage IS parameters as follows. First, we compute ${\mathbf{z}}_x$. Second, we draw a random sample of size $M_p$ from the Gaussian distribution with mean vector ${\mathbf{z}}_x$ and covariance matrix $\mathsf{\Sigma }$. Third, we construct $\overline{\theta }(x,·)$, the polynomial approximation to $\widehat{\theta }(x,·)$, as described in the previous paragraph. Fourth, for those sample points $\mathbf{z}$ that lie outside the region of interest we compute $q(x,\mathbf{z})$ using $\overline{\theta }$ instead $\widehat{\theta }$. The estimates of the optimal first-stage IS parameters are then:

$${\widehat{\mu }}_{IS}=\frac{{\sum }_{m=1}^{M_p}w\left({\mathbf{Z}}_m\right)exp(-Nq(x,{\mathbf{Z}}_m)){\mathbf{Z}}_m}{{\sum }_{m=1}^{M_p}w\left({\mathbf{Z}}_m\right)exp(-Nq(x,{\mathbf{Z}}_m))}$$

and

$${\widehat{\mathsf{\Sigma }}}_{IS}=\frac{{\sum }_{m=1}^{M_p}w\left({\mathbf{Z}}_m\right)exp(-Nq(x,{\mathbf{Z}}_m))({\mathbf{Z}}_m-{\widehat{\mu }}_{IS}){({\mathbf{Z}}_m-{\widehat{\mu }}_{IS})}^T}{{\sum }_{m=1}^{M_p}w\left({\mathbf{Z}}_m\right)exp(-Nq(x,{\mathbf{Z}}_m))},$$

where ${\mathbf{Z}}_1,\dots ,{\mathbf{Z}}_{M_p}$ is the random sample and

$$w\left(\mathbf{z}\right)=\frac{\varphi (\mathbf{z};\mathbf{0},\mathsf{\Sigma })}{\varphi (\mathbf{z};{\mathbf{z}}_x,\mathsf{\Sigma })}$$

is the IS weight associated with shifting the mean of the systematic risk factors from $\mathbf{0}$ to ${\mathbf{z}}_x$. The upper left panel of Figure 3 illustrates a typical situation where the mean of the IS distribution lies “just inside” the region of interest.

6.2.2. Computing Parameters in the One-Factor Model

The procedure described in the previous section specialises in the one-factor case as follows. First, under the parameter restrictions outlined in Section 5.3, the expected loss function $\mu \left(z\right)$ is a strictly decreasing function of z. As such, the region of interest is the semi-infinite interval $(-\infty ,z_x)$, where $z_x:={\mu }^{-1}\left(x\right)$, and its boundary is the single point $z_x$. In general $z_x$ must be computed numerically, which is straightforward. Second, we draw a random sample of size $M_p$ from the Gaussian distribution with mean $z_x$ and unit variance. Third, the polynomial approximation to $\widehat{\theta }$ is constructed by evaluating $\widehat{\theta }$ exactly (i.e., by numerically solving the equation $k^{\prime }(\theta ,z)=x$) at each of $n_g$ equally-spaced points z in the interval $[z_{-},z_{+}]$, where $z_{-}$ and $z_{+}$ are the largest and smallest values obtained in the pilot simulation, respectively, and then fitting a polynomial to the resulting $(z,\widehat{\theta }(x,z))$ pairs. Fourth, we evaluate $q(x,z)$ for each pilot sample point z as follows—if z lies inside the region of interest we set $q(x,z)=0$, otherwise we compute $q(x,z)$ by replacing the exact value $\widehat{\theta }(x,z)$ with the approximate value $\overline{\theta }(x,z)$, where $\overline{\theta }$ is the polynomial constructed in the previous step. Note that a single evaluation of $\overline{\theta }$ requires far less computational time than a single evaluation of $\widehat{\theta }$. Finally, the approximations to the first-stage IS parameters are:

$${\widehat{\mu }}_{IS}=\frac{{\sum }_{m=1}^{M_p}w\left(Z_m\right)exp(-Nq(x,Z_m))Z_m}{{\sum }_{m=1}^{M_p}w\left(Z_m\right)exp(-Nq(x,Z_m))}$$

and

$${\widehat{\sigma }}_{IS}^2=\frac{{\sum }_{m=1}^{M_p}w\left(Z_m\right)exp(-Nq(x,Z_m)){(Z_m-{\widehat{\mu }}_{IS})}^2}{{\sum }_{m=1}^{M_p}w\left(Z_m\right)exp(-Nq(x,Z_m))},$$

where $Z_1,\dots ,Z_{M_p}$ is the random sample and

$$w\left(z\right)=\frac{\varphi (z;0,1)}{\varphi (z;z_x,1)}.$$

is the IS weight associated with shifting the mean of the systematic risk factor from 0 to $z_x$.

6.2.3. Trimming Large Weights

In the one-factor model the first-stage IS weight will have infinite variance whenever ${\sigma }_{IS}^2<0.5$ (see Remark A1 in Appendix B). In a sample of 100 parameter sets, randomly selected according to the procedure in Section 5.3, the largest realised value of ${\sigma }_{IS}^2$ was $0.38$, and the mean and median were $0.11$ and $0.09$, respectively. It appears, then, that the first-stage IS weight in the one-factor model will have infinite variance in all cases of practical interest. We trim large weights as described in Section 4.2, using the set:

$$A=\left\{z\in \mathbb{R}:|z-{\widehat{\mu }}_{IS}|\le C{\widehat{\sigma }}_{IS}\right\}$$

for some constant C. In the numerical examples that follow we use $C=4$, in which case we expect to trim less than $0.01\%$ of the entire sample. Specialising Equation (34) to the present context, we get that an upper bound on the associated bias is given by:

$${\int }_{A^C}exp(-Nq(x,z))\varphi \left(z\right)dz,$$

(46)

which is straightforward (albeit slow) to compute using quadrature. Figure 4 illustrates the relationship between the probability of interest $p_x$ and the upper bound of Equation (46) for the 100 randomly generated parameter sets, and clearly demonstrates that the bias associated with our trimming procedure is negligible. For instance, for probabilities on the order of ${10}^{-3}$ the bias is no larger than ${10}^{-5}$, or 1% of the quantity of interest.

In the two-factor model the first-stage IS weight will have infinite variance whenever $\det (2{\mathsf{\Sigma }}^{-1}-{\mathsf{\Sigma }}_{IS}^{-1})<0$. In a random sample of 100 parameter sets, this condition occurred 96 times. As in the one-factor model, then, the first-stage IS weight in the two-factor model can be expected to have infinite variance in most cases of practical interest. We trim large weights using the set:

$$A=\left\{\mathbf{z}\in {\mathbb{R}}^2:{(\mathbf{z}-{\widehat{\mu }}_{IS})}^T{\widehat{\mathsf{\Sigma }}}_{IS}^{-1}(\mathbf{z}-{\widehat{\mu }}_{IS})|\le C^2\right\}$$

for some constant C, and use $C=4$ in the numerical examples that follow.

6.3. Second Stage

The first stage of the algorithm consists of (i) computing the first-stage IS parameters, (ii) simulating a random sample of size M from the systematic risk factors’ IS distribution, and (iii) computing the associated IS weights, trimming large weights appropriately. Having completed these tasks, the next step is to simulate individual losses in the second stage. In the remainder of this section we let $\mathbf{z}=(z_D,z_L)$ denote a generic realisation of the systematic risk factors obtained in the first stage.

6.3.1. Approximating $\widehat{\theta }$

Before generating any individual losses first construct the polynomial approximation to $\widehat{\theta }$, using the same procedure described in Section 6.2.1. The basic idea is to fit a relatively low degree polynomial to the surface of $\widehat{\theta }(x,·)$, over a small region that contains all of the first-stage sample points. The values of $\mathbf{z}$ obtained in the pilot sample are invariably different from those obtained in the first stage, so it is essential that the polynomial is refit to account for this fact. In what follows we use $\overline{\theta }$ to approximate $\widehat{\theta }$ whenever the numerical value of $\widehat{\theta }$ is required, but since the difference between the two is small we do not distinguish between the two (i.e., we write $\widehat{\theta }$ in this document, but use $\overline{\theta }$ in our code).

6.3.2. Sampling Individual Losses

In this section we describe how to sample individual losses in the two-factor model. The procedure carries over in an obvious way to the one-factor model, so we do not discuss that case explicitly.

If $\mathbf{z}$ lies inside the region of interest then the second stage is straightforward. For a given exposure i, we first simulate the exposure’s idiosyncratic risk factors ${\mathbf{Y}}_i=(Y_{i,D},Y_{i,L})$, from the bivariate normal distribution with standard normal margins and correlation ${\rho }_I$. Next, we set:

$$(X_{i,D},X_{i,L})=({\alpha }_D{z}_D+\sqrt{1-{\alpha }_D^2}{Y}_{i,D},{\alpha }_L{z}_L+\sqrt{1-{\alpha }_L^2}{Y}_{i,L}).$$

If $X_{i,D}>{\mathsf{\Phi }}^{-1}\left(P\right)$ then the exposure did not default and we set $L_i=0$ and proceed to the next exposure. Otherwise the exposure did default, in which case we must compute $h\left(X_{i,L}\right)$, set ${\ell }_i=h\left(x_{i,L}\right)$ and then proceed to the next exposure. Note that we only evaluate h for defaulted exposures—this is important since evaluating h requires numerical inversion of the beta cdf, which is relatively slow. Having computed the individual losses associated with each exposure, we then compute the average loss $\overline{\ell }=N^{-1}{\sum }_{i=1}^N{\ell }_i$ and set ${\mathsf{\Lambda }}_2(\mathbf{z},\overline{\ell })=1$.

If $\mathbf{z}$ lies outside the region of interest we must compute $\widehat{\theta }$, $k\left(\widehat{\theta }\right)$ and $\widehat{c}$, which we do approximately using the polynomial approximation $\overline{\theta }$. We then sample from ${\widehat{g}}_x(·|\mathbf{z})$ as follows. First simulate the idiosyncratic risk factors ${\mathbf{Y}}_i=(Y_{i,D},Y_{i,L})$ from the bivariate normal distribution with standard normal margins and correlation ${\rho }_I$. Also generate a random number U, independent of ${\mathbf{Y}}_i$. Then set:

$$(X_{i,D},X_{i,L})=({\alpha }_D{z}_D+\sqrt{1-{\alpha }_D^2}{Y}_{i,D},{\alpha }_L{z}_L+\sqrt{1-{\alpha }_L^2}{Y}_{i,L}).$$

If the exposure did not default we set ${\widehat{L}}_i=0$, otherwise we compute h and set ${\widehat{L}}_i=h\left(X_{i,L}\right)$. Next we check whether or not

$$U\le \frac{1}{\widehat{c}}·\frac{{\widehat{g}}_x\left({\widehat{L}}_i\right|\mathbf{z})}{g\left({\widehat{L}}_i\right|\mathbf{z})}=exp(-\widehat{\theta }({\ell }_{max}-{\widehat{L}}_i))$$

(47)

then accept ${\widehat{L}}_i$ as a drawing from ${\widehat{g}}_x$, that is, set $L_i={\widehat{L}}_i$ and proceed to exposure i. Otherwise, draw another random number U and set of idiosyncratic factors. Once we have sampled the individual losses associated with each exposure we compute the average loss $\overline{\ell }=N^{-1}{\sum }_{i=1}^N{\ell }_i$ and set ${\mathsf{\Lambda }}_2(\mathbf{z},\overline{\ell })=exp(-N[\widehat{\theta }\overline{\ell }-k(\widehat{\theta },\mathbf{z})])$, using the polynomial approximation to estimate the value of $\widehat{\theta }$.

6.3.3. Efficiency of the Second Stage

The frequency with which the rejection sampling algorithm must be applied in the second stage is governed by ${\mathbb{P}}_{IS}(\mu \left(\mathbf{Z}\right)>x)$. The left panel of Figure 5 illustrates the empirical distribution of this probability across 100 randomly selected parameter sets. The distribution is concentrated towards small values (the median fraction is 27%) but does have a relatively thick right tail (the mean fraction is 35%). In some cases—particularly when the value of the parameter ${\rho }_D$ is close to zero, in which case individual losses are very nearly independent and systematic risk is largely irrelevant—the vast majority of first-stage simulations require further IS in the second stage.

The efficiency of the rejection sampling algorithm, when it must be applied, is governed by the conditional distribution of $\widehat{c}=\widehat{c}(x,\mathbf{Z})$ given that $\mu \left(\mathbf{Z}\right)<x$. For each of the 100 parameter sets we estimate ${\mathbb{E}}_{IS}\left[\widehat{c}(x,\mathbf{Z})\right|\mu \left(\mathbf{Z}\right)<x]$, which determines the average size of the rejection constant for a given set of parameters, by computing the associated value of $\widehat{c}$ for each first-stage realisation that lies outside the region of interest and then averaging the resulting values. The right panel of Figure 5 illustrates the results, and we note that the mean and median of the data presented there are $1.09$ and $1.17$, respectively. The figure clearly indicates that the rejection sampling algorithm can be expected to be quite efficient, whenever it must be applied.

The distributions of ${\mathbb{P}}_{IS}(\mu \left(\mathbf{Z}\right)<x)$ and ${\mathbb{E}}_{IS}\left[\widehat{c}(x,\mathbf{Z})\right|\mu \left(\mathbf{Z}\right)<x]$ across parameters depend heavily on whether or not we adjust the variance of the systematic risk factors in the first stage. When we do not adjust variance, the mean and median of ${\mathbb{P}}_{IS}(\mu \left(\mathbf{Z}\right)<x)$ (across 100 randomly selected parameter sets) rise to 49% and 45% (as compared to 35% and 27% when we do adjust variance), and the mean and median of ${\mathbb{E}}_{IS}\left[\widehat{c}(x,\mathbf{Z})\right|\mu \left(\mathbf{Z}\right)<x]$ rise to $18.6$ and $1.8$, respectively (as compared to 1.17 and 1.09 when we do adjust variance).

Remark 7.

If we do not adjust the variance of the systematic risk factors in the first stage, then (i) the rejection sampling algorithm must be applied more frequently and (ii) is less efficient whenever it must be applied. As such, adjusting the variance of the systematic risk factors reduces the total time required to implement the two-stage algorithm.

The intuition behind this fact is as follows. First recall that the mean of the systematic risk factors tends to lie just inside the region of interest (recall Figure 3). In such cases the effect of reducing the variance of the systematic risk factors is to concentrate the distribution of $\mathbf{Z}$ just inside the boundary of the region of interest. Not only will this ensure that more first-stage realisations lie inside the region of interest (thereby reducing the fraction of points that require further IS in the second stage), it will also ensure that those realisations that lie outside the region (i.e., for which $\mu \left(\mathbf{z}\right)<x$) do not lie “that far” outside the region (i.e., that $\mu \left(\mathbf{z}\right)$ is not “that much less” than x), which in turn ensures that the typical size of $\widehat{c}$ is relatively close to one (recall the left panel of Figure 1).

7. Performance Evaluation

In this section we investigate the proposed algorithms’ performance in terms of statistical accuracy, computational time, and overall. Unless otherwise mentioned, we use a pilot sample size of $M_p=1000$ to estimate the first-stage IS parameters and a sample size of M = 10,000 to estimate the probability of interest ($p_x$). We use the value $C=4$ to trim large first-stage IS weights, and a value of $c_{max}=10$ to trim large rejection constants.

7.1. Statistical Accuracy

The standard error of any estimator that we consider is of the form ${\nu }_x/\sqrt{M}$ for some constant ${\nu }_x$ that depends on the algorithm used and the model parameters. For instance, for the one-stage estimator in the two-factor case we have ${\nu }_x={\mathrm{SD}}_{1\mathrm{S}}({\mathsf{\Lambda }}_1\left(\mathbf{Z}\right)·{\mathbf{1}}_{\{{\overline{L}}_N>x\}})$, where ${\mathrm{SD}}_{1\mathrm{S}}$ denotes standard deviation under the one-stage IS density of Equation (31). Note that in the absence of IS we have ${\nu }_x=\sqrt{p_x(1-p_x)}\sim p_x^{0.5}$ as $p_x\to 0$.

Figure 6 illustrates the relationship between ${\nu }_x$ and $p_x$ using 100 randomly selected parameters sets, for the two-stage algorithm and in the two-factor case. Importantly, we see that (i) ${\nu }_x$ seems to be a function of $p_x$ (i.e., it only depends on model parameters through $p_x$) and (ii) for small probabilities the functional relationship appears to be of the form ${\nu }_x=a{p}_x^b$ for constants a and b. These features are also present in the case of the one-stage estimator, as well as for both estimators in the one-factor model. The numerical values of a and b are easily estimated using the line of best fit (on the logarithmic scale), and the estimated values for both the one- and two-factor cases are summarised in

Table 1. This table reports fitted values of the relationship ${\nu }_x\approx a{p}_x^b$ for each estimator (one- and two-stage) and each model (one- and two-factor). Values of a and b are obtained by determining the line of best fit on the logarithmic scale (i.e., the line appearing in Figure 6). Note that in the absence of IS we would have ${\nu }_x=\sqrt{p_x(1-p_x)}\approx p_x^{0.5}$.

	One-Stage Algorithm	Two-Stage Algorithm
One-Factor Model	$0.91p_x^{0.98}$	$0.81p_x^{0.99}$
Two-Factor Model	$0.98p_x^{0.98}$	$0.81p_x^{0.98}$

. Of particular note is the fact that the value of b is extremely close to one in every case.

Of particular interest in the rare event context is an estimator’s relative error, defined as the ratio of its standard error to the true value of the quantity being estimated. For any of the estimators that we consider, the component of relative error that does not depend on sample size is ${\nu }_x/p_x\approx a{p}_x^{b-1}$. In the absence of IS we have $b-1=-0.5$, in which case relative error grows rapidly as $p_x\to 0$ (i.e., ${\nu }_x\to 0$ but ${\nu }_x/p_x\to \infty$ as $p_x\to 0$). By contrast, $b\approx 1$ for any of our IS estimators, in which case there is weak dependence of relative error on $p_x$. The minimum sample size required to ensure that an estimator’s relative error does not exceed the threshold $ϵ$ is $v_x^2/{\left(p_xϵ\right)}^2\approx a^2{p}_x^{2(b-1)}{ϵ}^{-2}$. In the absence of IS we have $b\approx 0.5$, in which case the sample size (and therefore computational burden) required to achieve a given degree of accuracy increases rapidly as $p_x\to 0$. By contrast, for all of our IS estimators we have $b\approx 1$, in which case the minimum sample size (and computational burden) is nearly independent of $p_x$.

Our ultimate goal is to reduce the computational burden associated with estimating $p_x$, in situations where $p_x$ is small. To see how effective the proposed algorithms are in this regard, note that the sample size required to achieve a given degree of accuracy using the proposed algorithm, relative to that required to achieve the same degree of accuracy in the absence of IS, is approximately

$$\frac{a^2{p}_x^{2(b-1)}{ϵ}^{-2}}{p_x^{-1}{ϵ}^{-2}}=a^2{p}_x^{2b-1},$$

which does not depend on $ϵ$. Since $a<1$ and $b>0.5$ (recall Table 1), we have that $a^2{p}_x^{2b-1}<p_x$.

Remark 8.

The relative sample size required to achieve a given degree of accuracy using the proposed algorithm, relative to that required in the absence of IS, is not larger than the probability of interest. For example, if the probability of interest is approximately 1%, then the proposed algorithm requires a sample size that is less than 1% of what would be required in the absence of IS (regardless of the desired degree of accuracy). And if the probability of interest is 0.1%, then the proposed algorithm requires a sample size that is less than 0.1% of what would be required in the absence of IS. In other words, the proposed algorithm is extremely effective at reducing the sample size required to achieve a given degree of accuracy.

It is also insightful to compare the efficiency of the two-stage estimator, relative to the one-stage estimator. In the one-factor case, the minimum sample size required using the two-stage algorithm, relative to that required using the one-stage algorithm, is approximately:

$$\frac{0.66p_x^{-0.02}{ϵ}^{-2}}{0.83p_x^{-0.04}{ϵ}^{-2}}=0.80p_x^{0.02}.$$

As $p_x$ ranges from $1\%$ to $0.01\%$ the estimated relative sample size ranges from $0.73$ to $0.67$. In the two-factor case, the relative sample size is approximately 0.69, regardless of the value of $p_x$.

Remark 9.

In both the one- and two-factor models, the two-stage algorithm is more efficient than the one-stage algorithm, in the sense that it requires a smaller sample size in order to achieve a given degree of accuracy. Indeed, in cases of practical interest (probabilities in the range of 1% to 0.01%) the minimum sample size required to achieve a given degree of accuracy using the two-stage algorithm is roughly 70% of what would be required using the one-stage algorithm.

7.2. Computational Time

Figure 7 illustrates the relationship between sample size (M) and run time (total time required to estimate $p_x$ using a particular algorithm), for one randomly selected set of parameters. Across both models and algorithms, the relationship is almost perfectly linear. In the absence of IS the intercept is zero (i.e., run time is directly proportional to sample size), whereas the intercepts are non-zero for the IS algorithms. The non-zero intercepts are due to the overhead associated with (i) computing the first-stage IS parameters, which accounts for almost all of the difference between the intercepts of the solid (no IS) and dashed (one-stage IS) intercepts, and (ii) computing the second-stage polynomial approximation to $\widehat{\theta }$, which accounts for almost all of the difference between the intercepts of the the dashed (one-stage IS) and dash-dot (two-stage IS) lines. It is also worth noting that a given increase in sample size will have a greater impact on the run times for the IS algorithms than it will on the standard algorithm. This is because we only calculate $h\left(X_{i,L}\right)$ for defaulted exposures (evaluating $h(·)$ is slow because it requires numerical inversion of the beta distribution function), and the default rate is higher under the IS distribution.

Across 100 randomly generated parameter sets, portfolio size (N) is most highly correlated with run time and the relationship is roughly linear.

Table 2. This table reports summary statistics—in seconds, and across 100 randomly selected parameter sets—for total run time (first three columns), time required to estimate the first-stage IS parameters (fourth column) and time required to fit the second-stage polynomial approximation to $\widehat{\theta }$ (final column).

Average Run Times
	No IS	One-Stage IS	Two-Stage IS	${\mathbf{\mu }}_{\mathbf{IS}},{\mathsf{\Sigma }}_{\mathbf{IS}}$	$\widehat{\mathbf{\theta }}$
One Factor	7.3	25.6	33.7	1.5	0.8
Two Factor	7.4	39.0	55.5	14.3	8.9

reports summary statistics on run times, across algorithms and models.

7.3. Overall Performance

Recall that the ultimate goal of this paper is to reduce the computational burden associated with estimating $p_x$, when $p_x$ is small. The computational burden associated with a particular algorithm is a function of both its statistical accuracy and total run time. We have seen that the proposed algorithms are substantially more accurate, but require considerably more run time. In this section we demonstrate that the benefit of increased accuracy is well worth the cost of additional run time, by considering the amount of time required by a particular algorithm in order to achieve a given degree of accuracy (as measured by relative error).

To begin, let $t\left(M\right)$ denote the total run time required by a particular algorithm to estimate $p_x$ using a sample of size M. As illustrated in Figure 7 we have $t\left(M\right)\approx c+dM$ for constants c and d that depend on the underlying model parameters (particularly portfolio size, N) as well as the algorithm being used. In Section 7.1 we saw that the minimum sample size required to ensure that the estimator’s relative error does not exceed the threshold $ϵ$ isL

$$M\left(ϵ\right)\approx a^2{p}_x^{2(b-1)}{ϵ}^{-2},$$

for constants a and b depending on the underlying model (one- or two-factor) and algorithm being used. Thus, if $T\left(ϵ\right)$ denotes the total CPU time required to ensure that the estimator’s relative error does not exceed $ϵ$, we have:

$$T\left(ϵ\right)\approx c+d{a}^2{p}_x^{2(b-1)}{ϵ}^{-2}.$$

(48)

Table 3. This table reports the time (in seconds) required to achieve a given degree of accuracy (computed using Equation (48)) for several values of $p_x$ and $ϵ$, for the parameter values corresponding to the left panel of Figure 7. Note that this is for the one-factor model. Values of c and d are obtained from the lines of best fit appearing in the left panel of Figure 7, values of a and b are obtained from Table 1.

No IS					One-Stage IS (Two-Stage IS)
	${\mathit{p}}_{\mathbf{x}}$	1%	0.1%	0.01%		${\mathit{p}}_{\mathbf{x}}$	1%	0.1%	0.01%
$\mathbf{ϵ}$					$\mathbf{ϵ}$
$10\%$		6	60	600	$10\%$		1.2 (2.3)	1.2 (2.3)	1.3 (2.4)
$5\%$		24	240	2400	$5\%$		1.8 (2.8)	1.9 (2.9)	1.9 (2.9)
$1\%$		600	6000	60,000	$1\%$		20.0 (18.8)	21.8 (19.6)	23.8 (20.4)

contains sample calculations for several different values of $p_x$ and $ϵ$, using the data appearing in the left panel of Figure 7 to estimate c and d and the values of a and b implicitly reported in Table 1. The results reported in the table are representative of those obtained using different parameter sets. It is clear that the proposed algorithms can substantially reduce the computational burden associated with accurate estimation of small probabilities. For instance, if the probability of interest is on the order of $0.1\%$ then either of the proposed algorithms can achieve 5% accuracy within 2–3 s, as compared to 4 min (80 times longer) in the absence of IS.

The two-stage estimator is statistically more accurate (Section 7.1) but computationally more expensive (Section 7.2) than the one-stage estimator. It is important to determine whether or not the benefit of increased accuracy outweighs the cost of increased computational time. Table 3 suggests that, in some cases at least, implementing the second stage is indeed worth the effort, in the sense that it can achieve the same degree of accuracy in less time.

Figure 8 illustrates the overall efficiency of the proposed algorithms, as a function of the desired degree of accuracy. Specifically, the left panel illustrates the ratio of (i) the total CPU time required to ensure the standard estimator’s relative error does not exceed a given threshold to (ii) the total time required by the proposed algorithms, for a randomly selected set of parameter values in the one-factor model. The right panel illustrates the same ratio for a randomly selected set of parameters in the two-factor model.

In the one-factor model, it would take hundreds of times longer to obtain an estimate of $p_x$ whose relative error is less than 10%, and thousands of times longer to obtain an estimate whose relative error is less than 1%. The figure also suggests that, since it requires less run time to obtain very accurate estimates, the two-stage algorithm is preferable to the one-stage algorithm in the one-factor model. In the two-factor model—where estimating IS parameters and fitting the second-stage polynomial approximation to $\widehat{\theta }$ is more time consuming—the proposed algorithms are hundreds of times more efficient than the standard algorithm. In addition, it appears that the one-stage algorithm is preferable to the two-stage algorithm in this case. Although the numerical values discussed here are specific to the parameter set used to produce the figure, they are representative of other parameter sets. In other words, the behaviour illustrated in Figure 8 is representative of the general framework overall.

8. Concluding Remarks

This paper developed an importance sampling (IS) algorithm for estimating large deviation probabilities for the loss on a portfolio of loans. In contrast to existing literature, we allowed loss given default to be stochastic and correlate with the default rate. The proposed algorithm proceeded in two stages. In the first stage one generates systematic risk factors from an IS distribution that is designed to increase the rate at which adverse macroeconomic scenarios are generated. In the second stage one checks whether or not the simulated macro environment is sufficiently adverse—if it is then no further IS is applied and idiosyncratic risk factors are drawn from their actual (conditional) probability distribution, if it is not then one indirectly applies IS to the conditional distribution of the idiosyncratic risk factors. Numerical evidence indicated that the proposed algorithm could be thousands of times more efficient than algorithms that did not employ any variance reduction techniques, across a wide variety of PD-LGD correlation models that are used in practice.