The Federal Reserve Board eagle logo links to home page

Skip to: [Printable Version (PDF)] [Bibliography] [Footnotes]
Finance and Economics Discussion Series: 2008-21 Screen Reader version

Nested Simulation in Portfolio Risk Measurement*

Michael B. Gordy
Federal Reserve Board
Sandeep Juneja
Tata Institute of Fundamental Research
8 April 2008



Keywords: Nested simulation, loss distribution, value-at-risk, expected shortfall, jackknife estimator, dynamic allocation

Abstract:

Risk measurement for derivative portfolios almost invariably calls for nested simulation. In the outer step one draws realizations of all risk factors up to the horizon, and in the inner step one re-prices each instrument in the portfolio at the horizon conditional on the drawn risk factors. Practitioners may perceive the computational burden of such nested schemes to be unacceptable, and adopt a variety of second-best pricing techniques to avoid the inner simulation. In this paper, we question whether such short cuts are necessary. We show that a relatively small number of trials in the inner step can yield accurate estimates, and analyze how a fixed computational budget may be allocated to the inner and the outer step to minimize the mean square error of the resultant estimator. Finally, we introduce a jackknife procedure for bias reduction and a dynamic allocation scheme for improved efficiency.

JEL Codes: G32, C15



For a wide variety of derivative instruments, computational costs may pose a binding constraint on the choice of pricing model. The more realistic and flexible the model, the less likely that there will exist an analytical pricing formula, and so the more likely that simulation-based pricing algorithms will be required. For plain-vanilla options trading in fast-moving markets, simulation is prohibitively slow. Simple models with analytical solutions are typically employed with ad-hoc adjustments (such as local volatility surfaces) to obtain better fit to the cross-section of market prices. As such models capture underlying processes in crude fashion, they tend to require frequent recalibration and perform poorly in time-series forecasting. For path-dependent options (e.g, lookback options) and complex basket derivatives (e.g., CDO of ABS), simulation is almost unavoidable, though even here computational shortcuts may be adopted at the expense of bias.1

Risk-management applications introduce additional challenges. Time constraints are less pressing than in trading applications, but the computational task may be more formidable. When loss is measured on a mark-to-market basis, estimation via simulation of large loss probabilities or of risk-measures such as Value-at-Risk (VaR) calls for a nested procedure: In the outer step one draws realizations of all risk factors up to the horizon, and in the inner step one re-prices each position in the portfolio at the horizon conditional on the drawn risk factors. It has been widely assumed that simulation-based pricing algorithms would be infeasible in the inner step, because the inner step must be executed once for each trial in the outer step.

In this paper, we question whether inner step simulations must necessarily impose a large computational burden. We show that a relatively small number of trials in the inner step can yield accurate estimates for large loss probabilities and portfolio risk-measures such as Value-at-Risk and Expected Shortfall, particularly when the portfolio contains a large number of positions. Since an expectation is replaced by a noisy sample mean, the estimator is biased, and we are able to characterize this bias asymptotically. We analyze how a fixed and large computational budget may be allocated to the inner and the outer step to minimize the mean square error of the resultant estimator. We show how the jackknifing technique may be applied to reduce the bias in our estimator, and how this alters the optimal budget allocation. In addition, we introduce a dynamic allocation scheme for choosing the number of inner step trials as a function of the generated output. This technique can significantly reduce the computational effort to achieve a given level of accuracy.

The most studied application of nested simulation in the finance literature is the pricing of American options. An influential paper by Longstaff and Schwartz (2001) proposes a least-squares methodology in which a small number of inner step samples are used to estimate a parametric relationship between the state vector at the horizon (in this case, the stock price) and the continuation value of the option. This "LSM" estimator is applicable to a broad range of nested problems, so long as the dimension of the state vector is not too large and the relationship between state vector and continuation value is not too nonlinear. However, some care must be taken in the choice of basis functions, and in general it may be difficult to assess the associated bias (Glasserman, 2004, §8.6). Our methodology, by contrast, is well-suited to portfolios of high-dimensional and highly nonlinear instruments, can be applied to a variety of derivative types without customization, and has bias of known form.

Our optimization results for large loss probabilities and Value-at-Risk are similar to those of Lee (1998).2 Lee's analysis relies on a different and somewhat more intricate set of assumptions than ours, which are in the spirit of the sensitivity analysis of VaR by Gouriéroux et al. (2000) and the subsequent literature on "granularity adjustment" of credit VaR (Gordy, 2004; Martin and Wilde, 2002). The resulting asymptotic formulae, however, are the same.3 Our extension of this methodology to Expected Shortfall is new, as is our analysis of large portfolio asymptotics. Furthermore, so far as we are aware, we are the first to examine the performance of jackknife estimators and dynamic allocation schemes in a nested simulation setting.

In Section 1 we set out a very general modeling framework for a portfolio of financial instruments. We introduce the nested simulation methodology in Section 2. We characterize the bias in and variance of the simulation estimator, and analyze the optimal allocation of computational resources between the two stages that minimizes the mean square error of the resultant estimator. Numerical illustrations of our main results are provided in Section 3. In the last two sections, we propose some refinements to further improve computational performance of nested simulation. Simple jackknife methods for bias reduction are developed in Section 4. Our dynamic allocation scheme is introduced and examined in Section 5.


1 Model framework

Let  X_t be a vector of  m state variables that govern all prices. The vector  X_t might include interest rates, commodity prices, equity prices, and other underlying prices referenced by derivatives. Let  {\cal F}_t be the filtration generated by  X_t. For use in discounting future cash flows, we denote by  B_t(s) the value at time  s of $1 invested at time  t\leq s in a risk free money market account, i.e.,

\displaystyle B_t(s) = \exp\left(\int_t^s r(u)du\right).
If interest rates are stochastic, then  B_t(s) depends on  {\cal F}_s.

The portfolio consists of  K+1 positions. The price of position  k at time  t depends on  t,  {\cal F}_t, and the contractual terms of the instrument.4 Position 0 represents the sub-portfolio of instruments for which there exist analytical pricing functions. Without loss of generality, we treat this as a single composite instrument. Among the contractual terms for an instrument is its maturity. We assume maturity  T_k is finite for  k=1,\ldots,K. As in all risk measurement exercises, the portfolio is assumed to be held static over the model horizon.

Conditional on  {\cal F}_t, the cashflows up to time  t are nonstochastic functions of time that depend on the contractual terms. Let  C_k(t) be the cumulative cashflow for  k on  (0,t]. Note that increments to  C_k(t) can be positive or negative, and can arrive at discrete time intervals or continuously. The market value of each position is the present discounted expected value of its cashflows under the risk-neutral measure  Q:

\displaystyle V_k(t) = E^Q\left\lbrack \int_t^{T_k} \frac{dC_k(s)}{B_t(s)} \bigg\vert {\cal F}_t \right\rbrack
The valuations are expressed in currency units, so there is no need for portfolio weights.5 We set  V_k(t)=0 whenever  t\geq T_k.

The present time is normalized to 0 and the model horizon is  H. "Loss" is defined as the difference between current value and discounted future value at the horizon, adjusting for interim cashflows. Portfolio loss is

\displaystyle Y = \sum_{k=0}^K \left(V_k(0) - \frac{1}{B_0(H)} \left(V_k(H) + \int_0^H B_t(H)\, dC_k(t)\right)\right).
The implicit assumption here is that interim cashflows are reinvested in the money market until time  H, but other conventions are easily accommodated.


2 Simulation framework

We now develop notation related to the simulation process. The simulation is nested: There is an "outer step" in which we draw histories up to the horizon  H. For each trial in the outer step, there is an "inner step" simulation needed for repricing at the horizon.

Let  L be the number of trials in the outer step. In each of these trials, we

  1. Draw a single path  X_t for  t=(0,H] under the physical measure. Let  \xi represent the filtration  {\cal F}_H that is generated by this path.
  2. Evaluate the accrued value at  H of the interim cashflows.
  3. Evaluate the price of each position at  H.
    • Closed-form price for instrument 0.
    • Simulation with  N_k "inner step" trials for remaining positions  k=1,\ldots,K. These paths are simulated under the risk-neutral measure.
  4. Discount back to time 0 to get our loss  Y(\xi).

Observe that the full dependence structure across the portfolio is captured in the period up to the model horizon. Inner step simulations, in contrast, are run independently across positions. This is because the value of position  k at time  H is simply a conditional expectation (given  {\cal F}_H and under the risk-neutral measure) of its own subsequent cash flows, and does not depend on future cash flows of other positions. Intuition might suggest that it would be more efficient from a simulation perspective to run inner step simulations simultaneously across all positions in order to reduce the total number of sampled paths of  X_t on  (H,\max\{T_k\}]. However, if we use the same samples of  X_t across inner step simulations, pricing errors are no longer independent across the positions, and so do not diversify away as effectively at the portfolio level. Furthermore, when the positions are repriced independently, to reprice position  k we need only draw joint paths for the elements of  X_t that influence that instrument. This may greatly reduce the memory footprint of the simulation, in particular when the number of state variables ( m) is large and when some of the maturities  \{T_k\} are very long relative to the horizon  H.

We have assumed that initial prices  V_k(0) are already known and can be taken as constants in our algorithm. Of course, this can be relaxed.

In the following three subsections, we discuss estimation of large loss probabilities (§2.1), Value-at-Risk (§2.2), and Expected Shortfall (§2.3). For simplicity, we impose a single value of  N across all positions (i.e.,  N_k=N for  k=1,...,K). This restriction is relaxed in Section 2.4. In Section 2.5, we consider the asymptotic behavior of the optimal allocation of computational resources as the portfolio size grows large. Last, in Section 2.6, we elaborate on the trade-offs associated with simultaneous repricing.


2.1 Estimating the probability of large losses

We first consider the problem of efficient estimation of  \alpha=P(Y(\xi)>u) via simulation for a given  u. If for each generated  \xi, the mark-to-market values of each position were known, the associated  Y(\xi) would be known and simulation would involve generating i.i.d. samples  Y(\xi_1), Y(\xi_2), \ldots, Y(\xi_L) and taking the average

\displaystyle \frac{1}{L} \sum_{i=1}^L \ensuremath{1[Y(\xi_i)>u]}
as an estimator of  \alpha. However, the mark-to-market value of each position is not known and is instead estimated via the inner step simulations.

Within the inner step simulation for repricing position  k, each trial gives an unbiased (but very noisy) estimate of  Y(\xi). Let  \ensuremath{{\zeta_{ki}}}(\xi) denote the zero-mean pricing error associated with the  i^{\rm th} such sample for position  k, let  Z_i(\xi) denote the portfolio pricing error for the  i^{\rm th} inner step sample, and finally let

\displaystyle \bar{Z}^\ensuremath{{\scriptscriptstyle N}}(\xi)= \frac{1}{N} \sum_{i=1}^N Z_i(\xi) = \frac{1}{N} \sum_{i=1}^N \sum_{k=1}^K \ensuremath{{\zeta_{ki}}}(\xi)
be the zero-mean average pricing error for the portfolio as a whole. In place of  Y(\xi), we take as its surrogate  {\tilde Y}(\xi)\equiv Y(\xi) + \bar{Z}^\ensuremath{{\scriptscriptstyle N}}(\xi) as an estimate of loss in the portfolio. By the law of large numbers,
\displaystyle \bar{Z}^\ensuremath{{\scriptscriptstyle N}}(\xi) \rightarrow 0 \quad a.s.
as  N \rightarrow \infty for any fixed  K (assuming that  \ensuremath{{\operatorname E}\lbrack \vert\ensuremath{{\zeta_{ki}}}\vert\rbrack}< \infty for each  k). The estimator for  P(Y(\xi)>u) then involves generating i.i.d. samples  ({\tilde Y}_1(\xi_1), \ldots, {\tilde Y}_L(\xi_L)) via outer and inner step simulation and taking the average
\displaystyle \ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}}= \frac{1}{L}\sum_{\ell=1}^L \ensuremath{1[{\tilde Y}_\ell(\xi_\ell)>u]}.

We now examine the mean square error of  \ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}}. Let  \alpha_\ensuremath{{\scriptscriptstyle N}} denote  P({\tilde Y}(\xi)>u). The mean square error of the estimator  \ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}} separates into

\displaystyle \ensuremath{{\operatorname E}\lbrack (\ensuremath{{\hat\alpha_\en... ...ptstyle N}})^2\rbrack}+ (\alpha_\ensuremath{{\scriptscriptstyle N}}-\alpha)^2.
Further, note that
\displaystyle \ensuremath{{\operatorname E}\lbrack (\ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}}- \alpha_\ensuremath{{\scriptscriptstyle N}})^2\rbrack} = \ensuremath{{\operatorname V}\lbrack \ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}}\rbrack} = \frac{\alpha_\ensuremath{{\scriptscriptstyle N}}(1-\alpha_\ensuremath{{\scriptscriptstyle N}})}{L}
where for any random variable  X,  \ensuremath{{\operatorname V}\lbrack X\rbrack} denotes its variance. Suppose that the computational effort to generate one sample of  \xi (outer step simulation) on an average equals  \gamma_0>0, and to generate an inner step simulation sample on the average equals  \gamma_1. Then, average effort per iteration of simulation equals  N \gamma_1 + \gamma_0 and average effort for  L iterations equals  L(N \gamma_1 + \gamma_0). By law of large numbers this is close to the actual effort when  L is large. We suppose that the overall computational budget is fixed at  \ensuremath{\Gamma}. We analyze the problem of choosing  (N,L) to minimize the mean square error of the estimator  \ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}} when  L(N \gamma_1 + \gamma_0)=\ensuremath{\Gamma} for large  \ensuremath{\Gamma}. Thus we consider the optimization problem
\displaystyle \min_{N,L \geq 0} \frac{\alpha_\ensuremath{{\scriptscriptstyle N}}(1-\alpha_\ensuremath{{\scriptscriptstyle N}})}{L} + (\alpha_\ensuremath{{\scriptscriptstyle N}}-\alpha)^2   subject to\displaystyle \quad L(N \gamma_1 + \gamma_0)=\ensuremath{\Gamma}, (1)

as  \ensuremath{\Gamma}\rightarrow \infty. Proposition 1 states our essential result for solving the optimization problem. We need some notation and a technical assumption for this.

Let  \tilde{Z}_\ensuremath{{\scriptscriptstyle N}}= \bar{Z}^\ensuremath{{\scriptscriptstyle N}}\cdot\sqrt{N} so that  \tilde{Z}_\ensuremath{{\scriptscriptstyle N}} has a non-trivial limit as  N \rightarrow \infty. Then  \alpha_\ensuremath{{\scriptscriptstyle N}}= P(Y+\tilde{Z}_\ensuremath{{\scriptscriptstyle N}}/\sqrt{N} >u). Our asymptotic analysis relies on Taylor series expansion of the joint density function  g_\ensuremath{{\scriptscriptstyle N}}(y,z) of  (Y, \tilde{Z}_\ensuremath{{\scriptscriptstyle N}}) and its partial derivatives. Assumption 1 ensures that higher order terms in such expansions can be ignored.

Assumption 1  
  1. The joint pdf  g_\ensuremath{{\scriptscriptstyle N}}(\cdot,\cdot) of  Y and  \tilde{Z}_\ensuremath{{\scriptscriptstyle N}} and its partial derivatives  \frac{\partial}{\partial y} g_\ensuremath{{\scriptscriptstyle N}}(y,z) and  \frac{\partial^2}{\partial y^2} g_\ensuremath{{\scriptscriptstyle N}}(y,z) exist for each  N and for all  (y,z).
  2. For  N \geq 1, there exist non-negative functions  p_{0,N}(\cdot),  p_{1,N}(\cdot) and  p_{2,N}(\cdot) such that
    \displaystyle g_\ensuremath{{\scriptscriptstyle N}}(y,z) \displaystyle \leq \displaystyle p_{0,N}(z),  
    \displaystyle \bigg\vert\frac{\partial}{\partial y} g_\ensuremath{{\scriptscriptstyle N}}(y,z)\bigg\vert \displaystyle \leq \displaystyle p_{1,N}(z)  
    \displaystyle \bigg\vert\frac{\partial^2}{\partial y^2} g_\ensuremath{{\scriptscriptstyle N}}(y,z)\bigg\vert \displaystyle \leq \displaystyle p_{2,N}(z)  

    for all  y,z. In addition,
    \displaystyle \sup_N \int_{-\infty}^{\infty}\vert z\vert^rp_{i,N}(z)dz < \infty
    for  i =0,1,2, and  0 \leq r \leq 4.

This assumption may be expected to be true in a large portfolio where there are at least a few positions that have a sufficiently smooth payoff. Alternatively, this assumption may be satisfied by perturbing  Y and  \tilde{Z}_N through adding to both of them mean zero, variance  \epsilon independent Gaussian random variables, also independent of  Y and  \tilde{Z}_N. For small  \epsilon this has a negligible impact on the tail measures. Then, if

\displaystyle \sup_N \int_{-\infty}^{\infty}z^4dF_{\tilde{Z}_N}(z) < \infty, (2)

where  F_{\tilde{Z}_N}(\cdot) denotes the distribution function of  \tilde{Z}_N, Assumption 1 can be seen to hold. To see this, let  Y_{\epsilon} and  \tilde{Z}_{N,\epsilon} denote the random variables obtained by perturbing  Y and  \tilde{Z}_N as described above. Then, the joint pdf of  (Y_{\epsilon},\tilde{Z}_{N,\epsilon}) equals
\displaystyle g_{\ensuremath{{\scriptscriptstyle N}},\epsilon}(y,z) = \int_{\Re^2}\phi_{\epsilon}(y-a) \phi_{\epsilon}(z-b) dG_\ensuremath{{\scriptscriptstyle N}}(a,b),
where  \phi_{\epsilon}(\cdot) denotes the pdf of Gaussian random variable with mean zero and variance  \epsilon, and  G_\ensuremath{{\scriptscriptstyle N}}(\cdot,\cdot) denotes the joint distribution function of  (Y, \tilde{Z}_\ensuremath{{\scriptscriptstyle N}}). Since,  \phi_{\epsilon}(\cdot) and its first two derivatives are bounded, it is easy to see that equation (2) implies that Assumption 1 holds for  (Y_{\epsilon},\tilde{Z}_{\ensuremath{{\scriptscriptstyle N}},\epsilon}).

Assumption 1 is sufficient to deliver a useful convergence property. Here and henceforth, let  f and  F denote the density and cumulative distribution function for  Y, and let  \ensuremath{\tilde{f}_\ensuremath{{\scriptscriptstyle N}}} and  \tilde{F}_\ensuremath{{\scriptscriptstyle N}} denote the density and cumulative distribution function for  \ensuremath{{\tilde Y}}. Now let  y_\ensuremath{{\scriptscriptstyle N}} be some sequence of real numbers that converges to a real number  y. In Appendix A.1, we prove the following lemma:

Lemma 1
Under Assumption 1, if  y_\ensuremath{{\scriptscriptstyle N}}\rightarrow y, then  \ensuremath{\tilde{f}_\ensuremath{{\scriptscriptstyle N}}}(y_\ensuremath{{\scriptscriptstyle N}}) \rightarrow f(y) and  \ensuremath{\tilde{f}_\ensuremath{{\scriptscriptstyle N}}}'(y_\ensuremath{{\scriptscriptstyle N}}) \rightarrow f'(y) as  N \rightarrow \infty.

We now approximate  \alpha_\ensuremath{{\scriptscriptstyle N}} in orders of  1/N. We define the function

\displaystyle \Theta(u) = \frac{1}{2} f(u) \ensuremath{{\operatorname E}\lbrack \sigma^2_{\xi}\vert Y=u\rbrack}, (3)

where  \sigma^2_{\xi} denotes the conditional variance of  Z_i (conditioned on  \xi). Our approximation is given by:6
Proposition 1   Under Assumption 1,
\displaystyle \alpha_\ensuremath{{\scriptscriptstyle N}}= \alpha + \theta_u/N + O_N(1/N^{3/2})
where  \theta_u=-\Theta'(u).
Proof is in Appendix A.2.

For the distributions and large loss levels one might expect to appear in practice, the bias will be upwards (i.e.,  \theta_u>0). By construction, the distribution of  Y(\xi)+\bar{Z}^\ensuremath{{\scriptscriptstyle N}}(\xi) differs from the distribution of  Y(\xi) by a mean-preserving spread, in the sense of Rothschild and Stiglitz (1970). Unless the two distributions have an infinite number of crossings, there will exist a  u^* such that  \ensuremath{{\hat\alpha_\ensuremath{{\scriptscriptstyle L\negthinspace{,}\negthinspace{N}}}}}\geq\alpha for all  u>u^*.

Applying Proposition 1, the objective function reduces to finding  N that minimizes

\displaystyle \frac{N\gamma_1 + \gamma_0}{\ensuremath{\Gamma}} \left (\ensuremath{\alpha}(1-\ensuremath{\alpha}) + O_N(1/N)\right ) + \frac{\theta_u^2}{N^2} + O_N(1/N^{5/2}).

It is easy to see that an optimal  N for this has the form

\displaystyle N^* = \left (\frac{2 \theta_u^2}{\ensuremath{\alpha}(1-\ensuremath{\alpha}) \gamma_1} \right )^{1/3} \ensuremath{\Gamma}^{1/3} +o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{1/3}). (4)

Therefore optimal  L has the form
\displaystyle L^* = \left( \frac{\ensuremath{\alpha}(1-\ensuremath{\alpha})} {2 \gamma_1^2 \theta_u^2} \right)^{1/3}\ensuremath{\Gamma}^{2/3} +o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{2/3}), (5)

and the mean square error at optimal  N^* equals
\displaystyle 3\left(\frac{\theta_u \ensuremath{\alpha}(1-\ensuremath{\alpha}) \gamma_1}{2 \ensuremath{\Gamma}}\right )^{2/3} + o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{-2/3}).

For large computational budgets, we see that  L^* grows with the square of  N^*. Thus, marginal increments to  \ensuremath{\Gamma} are allocated mainly to the outer step. It is easy to intuitively see the imbalance between  N^* and  L^*. Note that when  N and  L are of the same order  \sqrt{\ensuremath{\Gamma}}, the squared bias term contributes much less to the mean square error compared to the variance term. By increasing  L at the expense of  N we reduce variance till it matches up in contribution to the squared bias term.


2.2 Estimating Value-at-Risk

We now consider the problem of efficient estimation of Value-at-Risk for  Y. For a target insolvency probability  \ensuremath{\alpha}, VaR is the value  y_\ensuremath{\alpha} given by

\displaystyle y_\ensuremath{\alpha}= \ensuremath{{\rm VaR}_\ensuremath{\alpha}\lbrack Y\rbrack} = \inf\{y:P(Y\leq y)\geq 1-\ensuremath{\alpha}\}.

Under Assumption 1,  Y is a continuous random variable so that  P(Y \geq y_\ensuremath{\alpha})= \ensuremath{\alpha}. As before, our nested simulation generates samples  ({\tilde Y}_1(\xi_1), \ldots, {\tilde Y}_L(\xi_L)) where  {\tilde Y}(\xi)\equiv Y(\xi) + \bar{Z}^{N}(\xi). We sort these draws as  \tilde{Y}_{[1]} \geq \ldots \geq \tilde{Y}_{[L]}, so that  \tilde{Y}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}} provides an estimate of  y_\ensuremath{\alpha}, where  \lceil a \rceil denotes the integer ceiling of the real number  a. Our interest is in characterizing the mean square error  \ensuremath{{\operatorname E}\lbrack (\tilde{Y}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}}-y_\ensuremath{\alpha})^2\rbrack} and then minimizing it. As before, we decompose MSE into variance and squared bias:

\displaystyle \ensuremath{{\operatorname E}\lbrack (\tilde{Y}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}}-y_\ensuremath{\alpha})^2\rbrack} = \ensuremath{{\operatorname V}\lbrack \tilde{Y}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}}\rbrack} + \ensuremath{{\operatorname E}\lbrack \tilde{Y}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}}-y_\ensuremath{\alpha}\rbrack}^2.

To approximate bias and variance, we use the following result:

Proposition 2   Under Assumption 1,
\displaystyle \ensuremath{{\operatorname E}\lbrack \ensuremath{{\tilde Y}}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}}\rbrack}-y_\ensuremath{\alpha} = \frac{\theta_\ensuremath{\alpha}}{Nf(y_\ensuremath{\alpha})} +o_N(1/N) + O_L(1/L) + o_N(1)O_L(1/L),
where  \theta_\ensuremath{\alpha}=-\Theta'(y_\ensuremath{\alpha}) and
\displaystyle \ensuremath{{\operatorname V}\lbrack \ensuremath{{\tilde Y}}_\ensuremath{{\lceil \ensuremath{\alpha}L \rceil}}\rbrack} = \frac{\ensuremath{\alpha}(1-\ensuremath{\alpha})} {(L+2)f(y_\ensuremath{\alpha})^2} + O_L(1/L^2) + o_N(1)O_L(1/L).
A result parallel to the bias approximation is used in the literature on "granularity adjustment" of credit VaR to adjust asymptotic approximations of VaR for undiversified idiosyncratic risk (Gordy, 2004; Martin and Wilde, 2002). To avoid lengthy technical digressions, our statement of the proposition and its derivation in Appendix A.3 abstract from certain mild but cumbersome regularity conditions; see the appendix for details.

Our budget allocation problem reduces to minimizing the mean square error

\displaystyle \frac{\ensuremath{\alpha}(1-\ensuremath{\alpha})} {(L+2)f(y_\ensuremath{\alpha})^2} + \frac{\theta_\ensuremath{\alpha}^2}{N^2f(y_\ensuremath{\alpha})^2} +o_N(1/N^2) + o_N(1)O_L(1/L)+ O_L(1/L^2),
subject to  L(N \gamma_1 + \gamma_0)=\ensuremath{\Gamma}. It is easy to see that the optimal solution is
\displaystyle N^* = \left (\frac{2 \theta_\ensuremath{\alpha}^2}{\ensuremath{\alpha}(1-\ensuremath{\alpha}) \gamma_1} \right )^{1/3} \ensuremath{\Gamma}^{1/3} +o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{1/3}).
and
\displaystyle L^* = \left(\frac{\ensuremath{\alpha}(1-\ensuremath{\alpha})} {2 \gamma_1 ^2 \theta_\ensuremath{\alpha}^2}\right)^{1/3}\ensuremath{\Gamma}^{2/3} +o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{2/3}).
These values are identical up to terms of size  o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{2/3}) to the optimal values for estimating  P(Y>u) derived in the previous section when  u=y_\ensuremath{\alpha}. The mean square error at optimal  N^* equals
\displaystyle \frac{3}{f(y_\ensuremath{\alpha})^2}\left(\frac{\theta \ensuremath{\alpha}(1-\ensuremath{\alpha}) \gamma_1}{2 \ensuremath{\Gamma}}\right )^{2/3} + o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{-2/3}).


2.3 Estimating Expected Shortfall

Although Value-at-Risk is ubiquitous in industry practice, it is well understood that it has significant theoretical and practical shortcomings. It ignores the distribution of losses beyond the target quantile, so may give incentives to build portfolios that are highly sensitive to extreme tail events. More formally, Value-at-Risk fails to satisfy the sub-addivity property, so a merger of two portfolios can yield VaR greater than the sum of the two stand-alone VaRs. For this reason, Value-at-Risk is not a coherent risk-measure, in the sense of Artzner et al. (1999).

As an alternative to VaR, Acerbi and Tasche (2002) propose using generalized Expected Shortfall ("ES"), defined by

\displaystyle \ensuremath{{\rm ES}_\ensuremath{\alpha}\lbrack Y\rbrack}=\ensuremath{\alpha}^{-1}\left( \ensuremath{{\operatorname E}\lbrack Y\cdot\ensuremath{1[Y\geq y_\ensuremath{\alpha}]}\rbrack}+y_\ensuremath{\alpha}(1-\ensuremath{\alpha}-\Pr(Y<y_\ensuremath{\alpha}))\right), (6)

The first term is often used as the definition of Expected Shortfall for continuous variables. It is also known as "tail conditional expectations." The second term is a correction for mass at the quantile  y_\ensuremath{\alpha}. In our setting,  Y and  \ensuremath{{\tilde Y}} are continuous in distribution, so
\displaystyle \ensuremath{{\rm ES}_\ensuremath{\alpha}\lbrack Y\rbrack} \displaystyle = \displaystyle \frac{1}{\alpha}\ensuremath{{\operatorname E}\lbrack {Y}\cdot\ensuremath{1[{Y} > {y}_{\alpha}]}\rbrack}  
\displaystyle \ensuremath{{\rm ES}_\ensuremath{\alpha}\lbrack \ensuremath{{\tilde Y}}\rbrack} \displaystyle = \displaystyle \frac{1}{\alpha}\ensuremath{{\operatorname E}\lbrack {\ensuremath{{\tilde Y}}}\cdot\ensuremath{1[{\ensuremath{{\tilde Y}}} > \tilde{y}_{\alpha}]}\rbrack}  

Acerbi and Tasche (2002) show that ES is coherent and equivalent to the "conditional VaR" (CVaR) measure of Rockafellar and Uryasev (2002).

We begin with the more general problem of optimally allocating a computational budget to efficiently estimate  \Upsilon(u) = \ensuremath{{\operatorname E}\lbrack Y\cdot\ensuremath{1[Y>u]}\rbrack} for arbitrary  u. This is easier than the problem of estimating  \ensuremath{{\operatorname E}\lbrack Y\cdot \ensuremath{1[{Y} > {y}_{\alpha}]}\rbrack} since here  u is specified while in the latter case  {y}_{\alpha} is estimated. We return later to analyze the bias associated with the estimate of  \ensuremath{{\operatorname E}\lbrack Y\cdot\ensuremath{1[Y > {y}_{\alpha}]}\rbrack}.

Again, our sample output from the simulation to estimate  {Y\cdot \ensuremath{1[{Y}>u]}} equals  \ensuremath{{\tilde Y}}\cdot\ensuremath{1[\ensuremath{{\tilde Y}}>u]}. Let  \Upsilon_\ensuremath{{\scriptscriptstyle N}}(u) denote  \ensuremath{{\operatorname E}\lbrack \ensuremath{{\tilde Y}}\cdot\ensuremath{1[\ensuremath{{\tilde Y}}>u]}\rbrack}. The following proposition evaluates the bias associated with this term.

Proposition 3   Under Assumption 1,
\displaystyle \Upsilon_\ensuremath{{\scriptscriptstyle N}}(u) = \Upsilon(u) + \ensuremath{{(\Theta(u)-u\Theta'(u))}}/N + O_N(1/N^{3/2})
The proof is given in Appendix A.4.

Using similar analysis to the proof of Proposition 3, we can establish that

\displaystyle \ensuremath{{\operatorname E}\lbrack \ensuremath{{\tilde Y}}^2\cdot\ensuremath{1[\ensuremath{{\tilde Y}}>u]}\rbrack}- \ensuremath{{\operatorname E}\lbrack Y^2\cdot\ensuremath{1[Y>u]}\rbrack}= O_N(1/N),
and therefore that
\displaystyle \ensuremath{{\operatorname V}\lbrack \ensuremath{{\tilde Y}}\cdot\ensuremath{1[\ensuremath{{\tilde Y}}>u]}\rbrack} =\ensuremath{{\operatorname V}\lbrack Y \cdot\ensuremath{1[Y>u]}\rbrack}+O_N(1/N). (7)

Thus, analogous to Section 2.1, we consider the optimization problem
\displaystyle \min_{N,L \geq 0} \frac{\ensuremath{{\operatorname V}\lbrack \ensuremath{{\tilde Y}}\cdot\ensuremath{1[\ensuremath{{\tilde Y}}>u]}\rbrack}}{L} + (\Upsilon_\ensuremath{{\scriptscriptstyle N}}(u)-\Upsilon(u))^2   subject to\displaystyle \quad L(N \gamma_1 + \gamma_0)=\ensuremath{\Gamma},
as  \ensuremath{\Gamma}\rightarrow \infty.

Applying Proposition 3 and (7), the objective function reduces to finding  N that minimizes

\displaystyle \frac{V[Y\cdot\ensuremath{1[Y>u]}]+O_N(1/N) }{L} + \frac{\ensuremath{{(\Theta(u)-u\Theta'(u))}}^2}{N^2}+ O_{N}(1/N^{5/2}).
It is easy to see that an optimal  N for this has the form
\displaystyle N^* = \left (\frac{2 \ensuremath{{(\Theta(u)-u\Theta'(u))}}^2}{\ensuremath{{\operatorname V}\lbrack Y\cdot\ensuremath{1[Y>u]}\rbrack} \gamma_1} \right )^{1/3} \ensuremath{\Gamma}^{1/3} +o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{1/3}), (8)

and optimal  L has the form
\displaystyle L^* = \left( \frac{\ensuremath{{\operatorname V}\lbrack Y\cdot\ensuremath{1[Y>u]}\rbrack}} {2 \gamma_1^2 \ensuremath{{(\Theta(u)-u\Theta'(u))}}^2} \right)^{1/3}\ensuremath{\Gamma}^{2/3} +o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{2/3}). (9)

The mean square error at optimal  N^* equals
\displaystyle 3 \left(\frac{\ensuremath{{(\Theta(u)-u\Theta'(u))}}\ensuremath{{\operatorname V}\lbrack Y\cdot\ensuremath{1[Y>u]}\rbrack} \gamma_1}{2 \ensuremath{\Gamma}}\right )^{2/3} + o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{-2/3}).

We return now to the problem of the bias of  \ensuremath{{\rm ES}_\ensuremath{\alpha}\lbrack \ensuremath{{\tilde Y}}\rbrack}. We can write the difference between the Expected Shortfall of random variables  \ensuremath{{\tilde Y}} and  Y as

\begin{multline} \ensuremath{{\operatorname E}\lbrack \ensuremath{{\tilde Y}}\vert\ensuremath{{\tilde Y}}>\tilde{y}_{\alpha}\rbrack} - \ensuremath{{\operatorname E}\lbrack {Y}\vert{Y}>{y}_{\alpha}\rbrack} = \frac{1}{\alpha}(\Upsilon_\ensuremath{{\scriptscriptstyle N}}(\tilde{y}_{\alpha})- \Upsilon(y_{\alpha}))\ = \frac{1}{\alpha}\left( (\Upsilon_\ensuremath{{\scriptscriptstyle N}}(\tilde{y}_{\alpha})- \Upsilon_\ensuremath{{\scriptscriptstyle N}}({y}_{\alpha})) + (\Upsilon_\ensuremath{{\scriptscriptstyle N}}({y}_{\alpha}) - \Upsilon(y_{\alpha})) \right) \end{multline}

From Proposition 3, we have
\displaystyle \Upsilon_\ensuremath{{\scriptscriptstyle N}}({y}_{\alpha}) - \Upsilon(y_{\alpha}) =\frac{1 }{N}\left(\Theta(y_{\alpha})-y_{\alpha}\Theta'(y_{\alpha})\right) + O_N(1/N^{3/2}). (10)

Now from the mean value theorem
\displaystyle \Upsilon_\ensuremath{{\scriptscriptstyle N}}(\tilde{y}_{\alpha})- \Upsilon_\ensuremath{{\scriptscriptstyle N}}({y}_{\alpha}) = (\tilde{y}_{\alpha}-{y}_{\alpha}){\Upsilon_\ensuremath{{\scriptscriptstyle N}}}'(\breve{y})
where  \breve{y} lies between  {y}_{\alpha} and  \tilde{y}_{\alpha}. Note that  {\Upsilon_\ensuremath{{\scriptscriptstyle N}}}'({y})= -y \ensuremath{\tilde{f}_\ensuremath{{\scriptscriptstyle N}}}(y) and that  \breve{y} \rightarrow y_{\alpha} as  N \rightarrow \infty. From Lemma 1 it follows that  \ensuremath{\tilde{f}_\ensuremath{{\scriptscriptstyle N}}}(\breve{y}) \rightarrow f(y). Therefore,
\displaystyle \Upsilon_\ensuremath{{\scriptscriptstyle N}}(\tilde{y}_{\alpha})- \Upsilon_\ensuremath{{\scriptscriptstyle N}}({y}_{\alpha}) = -(\tilde{y}_{\alpha}-{y}_{\alpha})(y_{\alpha}f(y_{\alpha}) +o_N(1)) = - \frac{y_\ensuremath{\alpha}\theta_\ensuremath{\alpha}}{N} +o_N(1/N). (11)

where the last equality follows from Proposition 2. By substituting equations (11) and (12) into (10), we arrive at
\displaystyle \ensuremath{{\rm ES}_\ensuremath{\alpha}\lbrack \ensuremath{{\tilde Y}}\rbrack}- \ensuremath{{\rm ES}_\ensuremath{\alpha}\lbrack Y\rbrack} = \Theta(y_{\alpha})/N + o_N(1/N). (12)

A similar result is noted by Martin and Tasche (2007) and Gordy (2004).


2.4 Optimal allocation within the inner step

In this subsection we relax the restriction that  N_k is equal across  k. We focus on estimation of large loss probabilities. Similar analysis would allow us to vary  N_k across positions in estimating VaR and ES.

We redefine  N as the total number of inner step simulations. This aggregate  N is to be divided up among the positions  k=1,\ldots,K by allocating  p_k N simulations for position  k where each  p_k \geq 0 and  \sum_{k \leq K}p_k=1. Suppose that the average effort to generate a single such inner step simulation for  k is  \gamma_{1,k}. Then, total inner step simulation effort equals  N \gamma_1 where

\displaystyle \gamma_1 = \sum_{k=1}^K p_k\gamma_{1,k}.
For a single trial of the outer step simulation, the inner step simulation generates the estimate
\displaystyle \ensuremath{{\tilde Y}}(\xi) = Y(\xi) + \sum_{k=1}^K \frac{1}{p_k N}\sum_{i=1}^{p_k N} \ensuremath{{\zeta_{ki}}}(\xi).
Here we ignore minor technicalities associated with  p_k N not being an integer.

The analysis to compute the mean square error proceeds exactly as in Section 2.1. The resultant  \theta_u in this setting is

\displaystyle \theta_u = -\frac{1}{2} \frac{d}{du}f(u) \sum_{k=1}^K \frac{1}{p_k} \ensuremath{{\operatorname E}\lbrack \sigma^2_{k,\xi}\vert Y=u\rbrack} = \sum_{k=1}^K \frac{\theta_{u,k}}{p_k} (13)

where  \sigma^2_{k,\xi} denotes the variance of  \ensuremath{{\zeta_{k,\cdot}}}(\xi) conditioned on  \xi, and where we define
\displaystyle \theta_{u,k} = -\frac{1}{2} \frac{d}{du}f(u) \ensuremath{{\operatorname E}\lbrack \sigma^2_{k,\xi}\vert Y=u\rbrack}.

Recall from Section 2.1 that the mean square error at optimal  N^* equals

\displaystyle 3\left(\frac{\theta_u \alpha(1-\alpha) \gamma_1}{2\ensuremath{\Gamma}}\right)^{2/3} + o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{-2/3}).
Holding  N^* fixed, we now consider the problem of determining an approximation to optimal  \{p_k \geq 0: k \leq K\}. As  u and hence  \ensuremath{\alpha} are fixed, it is reasonable to ignore the  o_{\ensuremath{\Gamma}}(\ensuremath{\Gamma}^{-2/3}) residual and simply minimize the product  \theta_u \gamma_1, which we can write as
\displaystyle \theta_u \gamma_1 = \left(\sum_{k=1}^K \frac{\theta_{u,k}}{p_k}\right) \times\left(\sum_{j=1}^K p_j \gamma_{1,j}\right).
Since the terms  p_k and  p_j appear as ratios in the objective, the constraint  \sum_{k \leq K}p_k=1 simply involves normalizing any solution of the unconstrained problem. From the first order conditions, we can easily verify that the solution  (p_k^*: k \leq K) is
\displaystyle p_k^* = \frac{\sqrt{\theta_{u,k}/\gamma_{1,k}}} {\sum_{j=1}^K\sqrt{\theta_{u,j}/\gamma_{1,j}}}.
This is intuitive as one expects that higher computation resources should be allocated to a position with higher contribution to bias and lower computational effort. This is captured by  \theta_{u,k} in the numerator and  \gamma_{1,k} in the denominator.


2.5 Large portfolio asymptotics

Intuition suggests that as the portfolio size increases, the optimal number of inner loops needed becomes small, even falling to  N^*=1 for a sufficiently large portfolio. We formalize this intuition by considering an asymptotic framework where both the portfolio size and the computational budget increases to infinity. To avoid cumbersome notation and tedious technical arguments, we focus on the case of a portfolio of exchangeable (i.e., statistically homogeneous) positions. The arguments given are somewhat heuristic to give the flavor of analysis involved while avoiding the cumbersome and lengthy notation and assumptions needed to make it completely rigorous.

Consider an infinite sequence of exchangeable position indexed by  k, and let  W_k be the loss on position  k. Let  \bar{Y}^\ensuremath{{\scriptscriptstyle K}} be the average loss per position on a portfolio consisting of the first  K positions in the sequence, i.e.,

\displaystyle \bar{Y}^\ensuremath{{\scriptscriptstyle K}}= \frac{1}{K} \sum_{k=1}^K W_k.
We assume that Assumption 1 holds for the individual  W_k and their respective pricing errors  \ensuremath{{\zeta_{ki}}}.

As before, instead of observing  W_k(\xi), we generate  N inner step samples  W_k(\xi)+\ensuremath{{\zeta_{ki}}}(\xi) for  i=1, \ldots, N, so that our simulation provides an unbiased estimator for the probability

\displaystyle P(\bar{Y}^\ensuremath{{\scriptscriptstyle K}}+ \bar{Z}^{\ensuremath{{\scriptscriptstyle K,N}}}>u),
where
\displaystyle \bar{Z}^{\ensuremath{{\scriptscriptstyle K,N}}} = \frac{1}{KN} \sum_{i=1}^N\sum_{k=1}^K \ensuremath{{\zeta_{ki}}}.
We approximate the bias  P(\bar{Y}^\ensuremath{{\scriptscriptstyle K}}+ \bar{Z}^{\ensuremath{{\scriptscriptstyle K,N}}}>u)-P(\bar{Y}^\ensuremath{{\scriptscriptstyle K}}>u) as  K and  N grow large. Applying the same arguments as in proof of Proposition 1, we can show
\displaystyle P(\bar{Y}^\ensuremath{{\scriptscriptstyle K}}+ \bar{Z}^{\ensuremath{{\scriptscriptstyle K,N}}}>u) = P(\bar{Y}^\ensuremath{{\scriptscriptstyle K}}>u) + \theta_u^\ensuremath{{\scriptscriptstyle K}}/(KN) +o_{KN}(1/(KN)),
where  \theta_u^\ensuremath{{\scriptscriptstyle K}}= -\Theta_\ensuremath{{\scriptscriptstyle K}}'(u) and
\displaystyle \Theta_\ensuremath{{\scriptscriptstyle K}}(u) = \frac{1}{2} \frac{d}{du}\bar{f}_\ensuremath{{\scriptscriptstyle K}}(u) \ensuremath{{\operatorname E}\lbrack \tilde{\sigma}^2_{\xi}\vert\bar{Y}^\ensuremath{{\scriptscriptstyle K}}=u\rbrack}
where  \bar{f}_\ensuremath{{\scriptscriptstyle K}}(\cdot) denotes the pdf of  \bar{Y}^\ensuremath{{\scriptscriptstyle K}}, and  \tilde{\sigma}^2_{\xi} = \ensuremath{{\operatorname V}\lbrack \zeta_{\cdot,\cdot}\vert\xi\rbrack}.

By the law of large numbers,  \vert \bar{Y}^\ensuremath{{\scriptscriptstyle K}}-\ensuremath{{\operatorname E}\lbrack W_1\vert\xi\rbrack}\vert \rightarrow 0, almost surely, so the cdf  \bar{F}_\ensuremath{{\scriptscriptstyle K}} converges to a non-degenerate limiting distribution  \bar{F}, which is the distribution of  \ensuremath{{\operatorname E}\lbrack W_1\vert\xi\rbrack}. Similarly,  \ensuremath{{\operatorname E}\lbrack \tilde{\sigma}^2_{\xi}\vert\bar{Y}^\ensuremath{{\scriptscriptstyle K}}=u\rbrack} converges to  \ensuremath{{\operatorname E}\lbrack \tilde{\sigma}^2_{\xi}\vert\ensuremath{{\operatorname E}\lbrack W_1\vert\xi\rbrack}=u\rbrack}. Under suitable regularity conditions,  \theta_u^\ensuremath{{\scriptscriptstyle K}}\rightarrow\bar{\theta}_u=-\bar{\Theta}'(u), where

\displaystyle \bar{\Theta}(u) = \frac{1}{2} \frac{d}{du}\bar{f}(u) \ensuremath{{\operatorname E}\lbrack \tilde{\sigma}^2_{\xi}\vert\ensuremath{{\operatorname E}\lbrack W_1\vert\xi\rbrack}=u\rbrack}.
Therefore, the bias has the form  \bar{\theta}_u/(KN) + o_{KN}(1/KN).

We assume that the computational budget  \ensuremath{\Gamma}= \chi K^{\beta} for  \chi>0 and  \beta \geq 1. The value of  \beta captures the size of computational budget available relative to the time taken to generate a single inner loop sample. Note that if  \beta<1 then asymptotically even a single sample cannot be generated.

Recall that  m denotes the number of underlying state variables that control the prices of  K positions. Suppose that the computational effort to generate one sample of outer step simulation on average equals  \psi(m,K) for some function  \psi, and to generate an inner step simulation sample on the average equals  K \gamma_1 for a constant  \gamma_1>0. Average effort per outer step trial then equals  KN\gamma_1 +\psi(m,K) and average effort for  L such trials equals  L(KN\gamma_1 +\psi(m,K)). We analyze the order of magnitude of  N and  L that minimize the resultant mean square error of the estimator.

The mean square error of the estimator equals

\displaystyle \frac{\alpha(1-\alpha)}{L} + \frac{\bar{\theta}_u^2}{K^2N^2}
plus terms that are relatively negligible for large values of  K. Noting that  L= \frac{\chi K^{\beta-1}}{N \gamma_1 + \psi(m,K)/K}, it is easily seen that the value of  N \geq 1 that minimizes the dominant terms in the mean square error equals
\displaystyle N^* = \max \left(1, \left (\frac{2 \bar{\theta}_u^2 \chi}{\alpha (1- \alpha)\gamma_1} \right )^{1/3}K^{\beta/3-1} \right)
In particular, if  \beta <3 then,  N=1 for  K sufficiently large. Intuitively, this means that if the portfolio has a large number of positions and the computational budget is limited, then,  N may be kept equal to 1 irrespective of the form of the fixed-cost function  \psi(m,K). In this case  L^* is of order  \min(K^{\beta-1}, K^{\beta}/\psi(m,K)). Only when  \beta>3 does  N grow with  K.


2.6 Simultaneous repricing

Up to this point, we have stipulated that the inner step samples for each position in the portfolio are generated independently (conditional on  \xi) across positions. In application to derivative portfolios, there may be factors common to many positions (e.g., the prices of underlying securities) and it may be computationally efficient to generate these factors once for all positions rather than generating them independently for each position. While this reduces the computational effort required to generate a single sample of each position, it induces dependence across positions in the generated samples. If the dependence is such that the sum of the resultant noise from each position has lesser variance than if these samples were generated independently, as might be the case when there are many offsetting positions, then the former is a preferred method. However, typically the noises generated may have positive dependence and that may enhance the variance of the resulting samples thereby increasing the total number of samples required to achieve specified accuracy.

We now make this idea precise in a very simple setting. Consider the case where we want to find the expectation of  \sum_{i=1}^K X_i via simulation. Suppose that average computational effort needed to generate a sample of  \sum_{i=1}^K X_i by generating independent samples of  (X_1, \ldots, X_K) equals  K \gamma for some constant  \gamma>0. Let  \ensuremath{{\operatorname V}\lbrack X\rbrack} denote the variance of these  X_i's (to keep the discussion simple we assume that all rv have the same variance). Then, the computational effort required to get a specified accuracy is proportional to the variance of the sample  K\cdot \ensuremath{{\operatorname V}\lbrack X\rbrack} times the expected effort required to generate a single sample  K \gamma (see Glynn and Whitt, 1992), i.e.,  K^2 \gamma\ensuremath{{\operatorname V}\lbrack X\rbrack}. We refer to this measure as the simulation efficiency.

Now consider the case where we generate  \sum_{i=1}^K X_i by generating dependent samples of  (X_1, \ldots, X_K). Suppose that the computational effort to generate these samples on average equals  K\gamma (1-\beta) for some  \beta<1. Further suppose that the correlation between any two random variables  X_i and  X_j for  i \ne j is  \rho \in (0,1). Then the variance of  \sum_{i=1}^K X_i equals  K(1+ \rho(K-1)) \ensuremath{{\operatorname V}\lbrack X\rbrack}. So the simulation efficiency equals  K^2 \gamma (1-\beta)(1+ \rho(K-1))\ensuremath{{\operatorname V}\lbrack X\rbrack}. We therefore prefer to draw dependent samples whenever

\displaystyle \beta> 1- \frac{1}{(1+ \rho(K-1))}.
Unless  1-\beta=o_{\ensuremath{{\scriptscriptstyle K}}}(1/K), we will prefer to draw independent samples for any  K sufficiently large. This broadly indicates the benefits of independent samples in many finance settings. For the remainder of this paper, we reinstate the stipulation of independent repricing.


3 Numerical examples

We illustrate our results with a parametric example. Distributions for loss  Y and the pricing errors  Z are specified to ensure that the bias and variance of our simulation estimators are in closed-form. While the example is highly stylized, it allows us to compare our asymptotically optimal  (N^*,L^*) to the exact optimal solution under a finite computational budget. We have used simulation to perform similar exercises on the somewhat more realistic example of a portfolio of equity options. All our conclusions are robust.

Consider a homogeneous portfolio of  K positions. Let the state variable,  X_t, represent a single-dimensional market risk factor, and assume  X_H\sim{\cal N}(0,1). Let  \epsilon_k be the idiosyncratic component to the return on position  k at the horizon, so that the loss on position  k is  (X_H+\epsilon_k) per unit of exposure. We assume that the  \epsilon_k are i.i.d.   {\cal N}(0,\nu^2). To facilitate comparative statics on  K, we scale exposure sizes by  1/K.

The exact distribution for portfolio loss  Y is  {\cal N}(0,1+\nu^2/K). We assume that the position-level inner step pricing errors  \ensuremath{{\zeta_{k,\cdot}}} are i.i.d.   {\cal N}(0,\eta^2) per unit of exposure, so that the portfolio pricing error has variance  \sigma^2=\eta^2/K across inner step trials. This implies that the simulated loss variable