Nested Simulation in Portfolio Risk Measurement

Finance and Economics Discussion Series Divisions of Research & Statistics and Monetary Affairs Federal Reserve Board, Washington, D.C. Nested Simulation in Portfolio Risk Measurement Michael B. Gordy and Sandeep Juneja 2008-21 NOTE: Staff working papers in the Finance and Economics Discussion Series (FEDS) are preliminary materials circulated to stimulate discussion and critical comment. The analysis and conclusions set forth are those of the authors and do not indicate concurrence by other members of the research staff or the Board of Governors. References in publications to the Finance and Economics Discussion Series (other than acknowledgement) should be cleared with the author(s) to protect the tentative character of these papers.

Nested Simulation in Portfolio Risk Measurement∗ Michael B. Gordy Sandeep Juneja Federal Reserve Board Tata Institute of Fundamental Research 8 April 2008 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, wequestionwhethersuchshortcutsarenecessary. Weshowthatarelativelysmallnumberoftrials 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 Key words: nested simulation, loss distribution, value-at-risk, expected shortfall, jackknife estimator, dynamic allocation ∗WethankAaronBrownandDarrellDuffieforhelpfulsuggestions. Muchofthisworkwasconductedwhileboth authorswerevisitingtheIndianSchoolofBusiness. Theopinionsexpressedhereareourown,anddonotreflectthe views of the Board of Governors or its staff. Email: hmichael.gordy@frb.govi, hjuneja@tifr.res.ini.

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(suchaslocalvolatilitysurfaces)toobtainbetterfittothecross-sectionofmarketprices. 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-measuressuchasValue-at-Risk(VaR)callsforanestedprocedure: Intheoutersteponedraws 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 allocationschemeforchoosingthenumberofinnersteptrialsasafunctionofthegeneratedoutput. Thistechniquecansignificantlyreducethecomputationalefforttoachieveagivenlevelofaccuracy. 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 1ManyoftheideasinthispaperweredevelopedoriginallyinGordyandJuneja(2006)specificallyforapplication to portfolios of CDO tranches. 2Lee (1998) is an unpublished PhD thesis, which we encountered shortly before completion of this draft. Lee’s 1

than ours, which are in the spirit of the sensitivity analysis of VaR by Gouri´eroux et al. (2000) and the subsequent literature on “granularity adjustment” of credit VaR (Martin and Wilde, 2002; Gordy, 2004). 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. InSection1wesetoutaverygeneralmodelingframeworkforaportfoliooffinancialinstruments. 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 ofcomputational resources between the two stages that minimizes the mean square error of the resultant estimator. Numerical illustrationsofourmainresultsareprovidedinSection3. Inthelasttwosections, weproposesome refinements to further improve computational performance of nested simulation. Simple jackknife methodsforbiasreductionaredevelopedinSection4. Ourdynamicallocationschemeisintroduced and examined in Section 5. 1 Model framework Let X be a vector of m state variables that govern all prices. The vector X might include interest t t rates, commodity prices, equity prices, and other underlying prices referenced by derivatives. Let F be the filtration generated by X . For use in discounting future cash flows, we denote by B (s) t t t the value at time s of $1 invested at time t ≤ s in a risk free money market account, i.e., (cid:18)Z s (cid:19) B (s) = exp r(u)du . t t If interest rates are stochastic, then B (s) depends on F . t s The portfolio consists of K+1 positions. The price of position k at time t depends on t, F , and t 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 is finite for k = 1,...,K. As in all risk measurement exercises, the portfolio is k assumed to be held static over the model horizon. Conditional on F , the cashflows up to time t are nonstochastic functions of time that depend t onthecontractualterms. LetC (t)bethecumulativecashflowfork on(0,t]. Notethatincrements k to C (t) can be positive or negative, and can arrive at discrete time intervals or continuously. The k market value of each position is the present discounted expected value of its cashflows under the risk-neutral measure Q: V k (t) = EQ (cid:20)Z T k dC k (s) (cid:12) (cid:12) (cid:12)F t (cid:21) B (s) (cid:12) t t contribution anticipates the work of Gouri´eroux et al. (2000) and the literature on granularity adjustment, and appears to have been overlooked in the finance literature. 3Foraverydifferentapproachtooptimizationofcomputationalbudgetinanestedsimulation,seeDePriscoetal. (2007). 4For some exotic options, the price at t will depend on the entire path of X on s = (0,t]. This is why we need s the filtration F and not just X . t t 2

The valuations are expressed in currency units, so there is no need for portfolio weights.5 We set V (t) = 0 whenever t ≥ T . k 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 X K (cid:18) 1 (cid:18) Z H (cid:19)(cid:19) Y = V (0)− V (H)+ B (H)dC (t) . k k t k B (H) 0 0 k=0 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. DrawasinglepathX fort = (0,H]underthephysicalmeasure. Letξ representthefiltration t F that is generated by this path. H 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 “inner step” trials for remaining positions k = 1,...,K. These k paths are simulated under the risk-neutral measure. 4. Discount back to time 0 to get our loss Y(ξ). 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 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 on (H,max{T }]. However, if we use the same t k samples of X across inner step simulations, pricing errors are no longer independent across the t 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 that influence that instrument. This may greatly reduce the memory footprint of t the simulation, in particular when the number of state variables (m) is large and when some of the maturities {T } are very long relative to the horizon H. k 5Values can be negative, so it isn’t always natural to decompose portfolio value into a vector of weights and a vector of “returns”. 3

We have assumed that initial prices V (0) are already known and can be taken as constants in k our algorithm. Of course, this can be relaxed. Inthe followingthree subsections, we discuss estimationoflargelossprobabilities (§2.1), Valueat-Risk (§2.2), and Expected Shortfall (§2.3). For simplicity, we impose a single value of N across all positions (i.e., N = N for k = 1,...,K). This restriction is relaxed in Section 2.4. In Section k 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 Wefirstconsidertheproblemofefficientestimationofα = P(Y(ξ) > u)viasimulationforagivenu. Ifforeachgeneratedξ, themark-to-marketvaluesofeachpositionwereknown, theassociatedY(ξ) wouldbeknownandsimulationwouldinvolvegeneratingi.i.d. samplesY(ξ ),Y(ξ ),...,Y(ξ )and 1 2 L taking the average L 1 X 1[Y(ξ ) > u] i L i=1 as an estimator of α. 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(ξ). Let ζ (ξ) denote the zero-mean pricing error associated with the ith such ki sample for position k, let Z (ξ) denote the portfolio pricing error for the ith inner step sample, and i finally let N N K 1 X 1 XX Z¯N(ξ) = Z (ξ) = ζ (ξ) i ki N N i=1 i=1 k=1 be the zero-mean average pricing error for the portfolio as a whole. In place of Y(ξ), we take as its surrogate Y˜(ξ) ≡ Y(ξ)+Z¯N(ξ) as an estimate of loss in the portfolio. By the law of large numbers, Z¯N(ξ) → 0 a.s. as N → ∞ for any fixed K (assuming that E[|ζ |] < ∞ for each k). The estimator for P(Y(ξ) > u) ki then involves generating i.i.d. samples (Y˜ (ξ ),...,Y˜ (ξ )) via outer and inner step simulation and 1 1 L L taking the average L 1 X αˆ = 1[Y˜ (ξ ) > u]. L,N L ‘ ‘ ‘=1 We now examine the mean square error of αˆ . Let α denote P(Y˜(ξ) > u). The mean square L,N N error of the estimator αˆ separates into L,N E[(αˆ −α)2] = E[(αˆ −α +α −α)2] = E[(αˆ −α )2]+(α −α)2. L,N L,N N N L,N N N Further, note that α (1−α ) E[(αˆ −α )2] = V[αˆ ] = N N L,N N L,N L whereforanyrandomvariableX, V[X]denotesitsvariance. Supposethatthecomputationaleffort to generate one sample of ξ (outer step simulation) on an average equals γ > 0, and to generate 0 4

an inner step simulation sample on the average equals γ . Then, average effort per iteration of 1 simulation equals Nγ +γ and average effort for L iterations equals L(Nγ +γ ). By law of large 1 0 1 0 numbersthisisclosetotheactualeffortwhenLislarge. Wesupposethattheoverallcomputational budget is fixed at Γ. We analyze the problem of choosing (N,L) to minimize the mean square error oftheestimatorαˆ whenL(Nγ +γ ) = ΓforlargeΓ. Thusweconsidertheoptimizationproblem L,N 1 0 α (1−α ) min N N +(α −α)2 subject to L(Nγ +γ ) = Γ, (1) N,L≥0 L N 1 0 as Γ → ∞. Proposition 1 states our essential result for solving the optimization problem. We need some notation and a technical assumption for this. √ √ Let Z˜ = Z¯N· N so that Z˜ has a non-trivial limit as N → ∞. Then α = P(Y +Z˜ / N > N N N N u). OurasymptoticanalysisreliesonTaylorseriesexpansionofthejointdensityfunctiong (y,z)of N (Y,Z˜ )anditspartialderivatives. Assumption1ensuresthathigherordertermsinsuchexpansions N can be ignored. Assumption 1 1. The joint pdf g (·,·) of Y and Z˜ and its partial derivatives ∂ g (y,z) and N N ∂y N ∂2 g (y,z) exist for each N and for all (y,z). ∂y2 N 2. For N ≥ 1, there exist non-negative functions p (·), p (·) and p (·) such that 0,N 1,N 2,N g (y,z) ≤ p (z), N 0,N (cid:12) (cid:12) (cid:12) ∂ (cid:12) (cid:12) (cid:12)∂y g N (y,z)(cid:12) (cid:12) ≤ p 1,N (z) (cid:12) (cid:12) ∂2 (cid:12) (cid:12) (cid:12) (cid:12)∂y2 g N (y,z)(cid:12) (cid:12) ≤ p 2,N (z) for all y,z. In addition, Z ∞ sup |z|rp (z)dz < ∞ i,N N −∞ for i = 0,1,2, and 0 ≤ r ≤ 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 Z˜ through adding to both of them mean zero, variance (cid:15) independent Gaussian N random variables, also independent of Y and Z˜ . For small (cid:15) this has a negligible impact on the N tail measures. Then, if Z ∞ sup z4dF (z) < ∞, (2) N −∞ Z˜ N where F (·) denotes the distribution function of Z˜ , Assumption 1 can be seen to hold. To see Z˜ N N this, let Y and Z˜ denote the random variables obtained by perturbing Y and Z˜ as described (cid:15) N,(cid:15) N above. Then, the joint pdf of (Y ,Z˜ ) equals (cid:15) N,(cid:15) Z g (y,z) = φ (y−a)φ (z−b)dG (a,b), N,(cid:15) (cid:15) (cid:15) N <2 5

where φ (·) denotes the pdf of Gaussian random variable with mean zero and variance (cid:15), and (cid:15) G (·,·) denotes the joint distribution function of (Y,Z˜ ). Since, φ (·) and its first two derivatives N N (cid:15) are bounded, it is easy to see that equation (2) implies that Assumption 1 holds for (Y ,Z˜ ). (cid:15) N,(cid:15) 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 f˜ and F˜ denote the N N density and cumulative distribution function for Y˜. Now let y be some sequence of real numbers N that converges to a real number y. In Appendix A.1, we prove the following lemma: Lemma 1 Under Assumption 1, if y → y, then f˜ (y ) → f(y) and f˜0 (y ) → f0(y) as N → ∞. N N N N N We now approximate α in orders of 1/N. We define the function N 1 Θ(u) = f(u)E[σ2|Y = u], (3) 2 ξ where σ2 denotes the conditional variance of Z (conditioned on ξ). Our approximation is given ξ i by:6 Proposition 1 Under Assumption 1, α = α+θ /N +O (1/N3/2) N u N where θ = −Θ0(u). 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., θ > 0). By construction, the distribution of Y(ξ)+Z¯N(ξ) differs from the u distribution of Y(ξ) 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 αˆ ≥ α for all u > u∗. L,N Applying Proposition 1, the objective function reduces to finding N that minimizes Nγ +γ θ2 1 0 (α(1−α)+O (1/N))+ u +O (1/N5/2). Γ N N2 N It is easy to see that an optimal N for this has the form (cid:18) 2θ2 (cid:19)1/3 N∗ = u Γ1/3+o (Γ1/3). (4) Γ α(1−α)γ 1 Therefore optimal L has the form (cid:18) α(1−α) (cid:19)1/3 L∗ = Γ2/3+o (Γ2/3), (5) 2γ2θ2 Γ 1 u and the mean square error at optimal N∗ equals (cid:18) θ α(1−α)γ (cid:19)2/3 3 u 1 +o (Γ−2/3). Γ 2Γ 6WesaythatafunctionisO (h(m))ifforallmsufficientlylargeitsabsolutevalueisupperboundedbyaconstant m times h(m). We say that it is o (h(m)) if for any (cid:15)>0 its absolute value is upper bounded by (cid:15) times h(m). m 6

For large computational budgets, we see that L∗ grows with the square of N∗. Thus, marginal increments to Γ 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 Γ, 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 α, VaR is the value y given by α y = VaR [Y] = inf{y : P(Y ≤ y) ≥ 1−α}. α α Under Assumption 1, Y is a continuous random variable so that P(Y ≥ y ) = α. As before, α our nested simulation generates samples (Y˜ (ξ ),...,Y˜ (ξ )) where Y˜(ξ) ≡ Y(ξ) + Z¯N(ξ). We 1 1 L L sort these draws as Y˜ ≥ ... ≥ Y˜ , so that Y˜ provides an estimate of y , where dae denotes [1] [L] dαLe α the integer ceiling of the real number a. Our interest is in characterizing the mean square error E[(Y˜ −y )2] and then minimizing it. As before, we decompose MSE into variance and squared dαLe α bias: E[(Y˜ −y )2] = V[Y˜ ]+E[Y˜ −y ]2. dαLe α dαLe dαLe α To approximate bias and variance, we use the following result: Proposition 2 Under Assumption 1, θ E[Y˜ ]−y = α +o (1/N)+O (1/L)+o (1)O (1/L), dαLe α Nf(y ) N L N L α where θ = −Θ0(y ) and α α α(1−α) V[Y˜ ] = +O (1/L2)+o (1)O (1/L). dαLe (L+2)f(y )2 L N 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 (Martin and Wilde, 2002; Gordy, 2004). 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 α(1−α) θ2 + α +o (1/N2)+o (1)O (1/L)+O (1/L2), (L+2)f(y )2 N2f(y )2 N N L L α α subject to L(Nγ +γ ) = Γ. It is easy to see that the optimal solution is 1 0 (cid:18) 2θ2 (cid:19)1/3 N∗ = α Γ1/3+o (Γ1/3). Γ α(1−α)γ 1 and (cid:18) α(1−α) (cid:19)1/3 L∗ = Γ2/3+o (Γ2/3). 2γ2θ2 Γ 1 α 7

Thesevaluesareidenticaluptotermsofsizeo (Γ2/3)totheoptimalvaluesforestimatingP(Y > u) Γ derived in the previous section when u = y . The mean square error at optimal N∗ equals α 3 (cid:18) θα(1−α)γ (cid:19)2/3 1 +o (Γ−2/3). f(y )2 2Γ Γ α 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. Moreformally, Value-at-Riskfailstosatisfythesub-addivityproperty, soamergeroftwoportfolios 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). AsanalternativetoVaR,AcerbiandTasche(2002)proposeusinggeneralizedExpected Shortfall (“ES”), defined by ES [Y] = α−1(E[Y ·1[Y ≥ y ]]+y (1−α−Pr(Y < y ))), (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 . In our setting, Y and Y˜ are continuous in distribution, so α 1 ES [Y] = E[Y ·1[Y > y ]] α α α 1 ES [Y˜] = E[Y˜ ·1[Y˜ > y˜ ]] α α α AcerbiandTasche(2002)showthatESiscoherentandequivalenttothe“conditionalVaR”(CVaR) measure of Rockafellar and Uryasev (2002). We begin with the more general problem of optimally allocating a computational budget to efficiently estimate Υ(u) = E[Y · 1[Y > u]] for arbitrary u. This is easier than the problem of estimating E[Y ·1[Y > y ]] since here u is specified while in the latter case y is estimated. We α α return later to analyze the bias associated with the estimate of E[Y ·1[Y > y ]]. α Again, our sample output from the simulation to estimate Y ·1[Y > u] equals Y˜ ·1[Y˜ > u]. Let Υ (u) denote E[Y˜ ·1[Y˜ > u]]. The following proposition evaluates the bias associated with this N term. Proposition 3 Under Assumption 1, Υ (u) = Υ(u)+(Θ(u)−uΘ0(u))/N +O (1/N3/2) N N The proof is given in Appendix A.4. Using similar analysis to the proof of Proposition 3, we can establish that E[Y˜2·1[Y˜ > u]]−E[Y2·1[Y > u]] = O (1/N), N and therefore that V[Y˜ ·1[Y˜ > u]] = V[Y ·1[Y > u]]+O (1/N). (7) N 8

Thus, analogous to Section 2.1, we consider the optimization problem V[Y˜ ·1[Y˜ > u]] min +(Υ (u)−Υ(u))2 subject to L(Nγ +γ ) = Γ, N,L≥0 L N 1 0 as Γ → ∞. Applying Proposition 3 and (7), the objective function reduces to finding N that minimizes V[Y ·1[Y > u]]+O (1/N) (Θ(u)−uΘ0(u))2 N + +O (1/N5/2). L N2 N It is easy to see that an optimal N for this has the form !1/3 2(Θ(u)−uΘ0(u))2 N∗ = Γ1/3+o (Γ1/3), (8) Γ V[Y ·1[Y > u]]γ 1 and optimal L has the form !1/3 V[Y ·1[Y > u]] L∗ = Γ2/3+o (Γ2/3). (9) 2γ2(Θ(u)−uΘ0(u))2 Γ 1 The mean square error at optimal N∗ equals (cid:18) (Θ(u)−uΘ0(u))V[Y ·1[Y > u]]γ (cid:19)2/3 3 1 +o (Γ−2/3). Γ 2Γ We return now to the problem of the bias of ES [Y˜]. We can write the difference between the α Expected Shortfall of random variables Y˜ and Y as 1 E[Y˜|Y˜ > y˜ ]−E[Y|Y > y ] = (Υ (y˜ )−Υ(y )) α α α N α α 1 = ((Υ (y˜ )−Υ (y ))+(Υ (y )−Υ(y ))) (10) α N α N α N α α From Proposition 3, we have Υ (y )−Υ(y ) = 1 (cid:0) Θ(y )−y Θ0(y ) (cid:1) +O (1/N3/2). (11) N α α N α α α N Now from the mean value theorem Υ (y˜ )−Υ (y ) = (y˜ −y )Υ 0(y˘) N α N α α α N where y˘ lies between y and y˜ . Note that Υ 0(y) = −yf˜ (y) and that y˘→ y as N → ∞. From α α N N α Lemma 1 it follows that f˜ (y˘) → f(y). Therefore, N y θ α α Υ (y˜ )−Υ (y ) = −(y˜ −y )(y f(y )+o (1)) = − +o (1/N). (12) N α N α α α α α N N N where the last equality follows from Proposition 2. By substituting equations (11) and (12) into (10), we arrive at ES [Y˜]−ES [Y] = Θ(y )/N +o (1/N). (13) α α α N A similar result is noted by Martin and Tasche (2007) and Gordy (2004). 9

2.4 Optimal allocation within the inner step In this subsection we relax the restriction that N is equal across k. We focus on estimation of k large loss probabilities. Similar analysis would allow us to vary N across positions in estimating k 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,...,K by allocating p N simulations for position k where each k P p ≥ 0 and p = 1. Suppose that the average effort to generate a single such inner step k k≤K k simulation for k is γ . Then, total inner step simulation effort equals Nγ where 1,k 1 K X γ = p γ . 1 k 1,k k=1 For a single trial of the outer step simulation, the inner step simulation generates the estimate X K 1 p Xk N Y˜(ξ) = Y(ξ)+ ζ (ξ). ki p N k k=1 i=1 Here we ignore minor technicalities associated with p N not being an integer. k TheanalysistocomputethemeansquareerrorproceedsexactlyasinSection2.1. Theresultant θ in this setting is u K K θ = − 1 d f(u) X 1 E[σ2 |Y = u] = X θ u,k (14) u 2du p k,ξ p k k k=1 k=1 where σ2 denotes the variance of ζ (ξ) conditioned on ξ, and where we define k,ξ k,· 1 d θ = − f(u)E[σ2 |Y = u]. u,k 2du k,ξ Recall from Section 2.1 that the mean square error at optimal N∗ equals (cid:18) θ α(1−α)γ (cid:19)2/3 3 u 1 +o (Γ−2/3). Γ 2Γ Holding N∗ fixed, we now consider the problem of determining an approximation to optimal {p ≥ k 0 : k ≤ K}. As u and hence α are fixed, it is reasonable to ignore the o (Γ−2/3) residual and simply Γ minimize the product θ γ , which we can write as u 1   K ! K X θ u,k X θ u γ 1 = × p j γ 1,j. p k k=1 j=1 P Since the terms p and p appear as ratios in the objective, the constraint p = 1 simply k j k≤K k involves normalizing any solution of the unconstrained problem. From the first order conditions, we can easily verify that the solution (p∗ : k ≤ K) is k p θ /γ p∗ = u,k 1,k . k PK p θ /γ j=1 u,j 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 θ in the u,k numerator and γ in the denominator. 1,k 10

2.5 Large portfolio asymptotics Intuition suggests that as the portfolio size increases, the optimal number of inner loops needed becomessmall,evenfallingtoN∗ = 1forasufficientlylargeportfolio. Weformalizethisintuitionby 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 be the loss on k position k. Let Y¯K be the average loss per position on a portfolio consisting of the first K positions in the sequence, i.e., K 1 X Y¯K = W . k K k=1 We assume that Assumption 1 holds for the individual W and their respective pricing errors ζ . k ki As before, instead of observing W (ξ), we generate N inner step samples W (ξ) + ζ (ξ) for k k ki i = 1,...,N, so that our simulation provides an unbiased estimator for the probability P(Y¯K +Z¯K,N > u), where N K 1 XX Z¯K,N = ζ . ki KN i=1 k=1 We approximate the bias P(Y¯K +Z¯K,N > u)−P(Y¯K > u) as K and N grow large. Applying the same arguments as in proof of Proposition 1, we can show P(Y¯K +Z¯K,N > u) = P(Y¯K > u)+θK/(KN)+o (1/(KN)), u KN where θK = −Θ0 (u) and u K 1 d Θ (u) = f¯ (u)E[σ˜2|Y¯K = u] K 2du K ξ where f¯ (·) denotes the pdf of Y¯K, and σ˜2 = V[ζ |ξ]. K ξ ·,· Bythelawoflargenumbers,|Y¯K−E[W |ξ]| → 0,almostsurely,sothecdfF¯ convergestoanon- 1 K degenerate limiting distribution F¯, which is the distribution of E[W |ξ]. Similarly, E[σ˜2|Y¯K = u] 1 ξ converges to E[σ˜2|E[W |ξ] = u]. Under suitable regularity conditions, θK → θ¯ = −Θ¯0(u), where ξ 1 u u 1 d Θ¯(u) = f¯(u)E[σ˜2|E[W |ξ] = u]. 2du ξ 1 Therefore, the bias has the form θ¯ /(KN)+o (1/KN). u KN We assume that the computational budget Γ = χKβ for χ > 0 and β ≥ 1. The value of β captures the size of computational budget available relative to the time taken to generate a single innerloopsample. Notethatifβ < 1thenasymptoticallyevenasinglesamplecannotbegenerated. 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 ψ(m,K) for some function ψ, and to generate an inner step simulation sample 11

on the average equals Kγ for a constant γ > 0. Average effort per outer step trial then equals 1 1 KNγ +ψ(m,K) and average effort for L such trials equals L(KNγ +ψ(m,K)). We analyze the 1 1 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 α(1−α) θ¯2 + u L K2N2 plus terms that are relatively negligible for large values of K. Noting that L = χKβ−1 , it is Nγ1+ψ(m,K)/K easily seen that the value of N ≥ 1 that minimizes the dominant terms in the mean square error equals ! (cid:18) 2θ¯2χ (cid:19)1/3 N∗ = max 1, u Kβ/3−1 α(1−α)γ 1 In particular, if β < 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 ψ(m,K). In this case L∗ is of order min(Kβ−1,Kβ/ψ(m,K)). Only when β > 3 does N grow with K. 2.6 Simultaneous repricing Uptothispoint,wehavestipulatedthattheinnerstepsamplesforeachpositionintheportfolioare generated independently (conditional on ξ) across positions. In application to derivative portfolios, theremaybefactorscommontomanypositions(e.g.,thepricesofunderlyingsecurities)anditmay 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 PK X via simulation. Suppose that average computational effort needed i=1 i to generate a sample of PK X by generating independent samples of (X ,...,X ) equals Kγ i=1 i 1 K for some constant γ > 0. Let V[X] denote the variance of these X ’s (to keep the discussion simple i 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·V[X] times the expected effort required to generate a single sample Kγ (see Glynn and Whitt, 1992), i.e., K2γV[X]. We refer to this measure as the simulation efficiency. Nowconsiderthecasewherewegenerate PK X bygeneratingdependentsamplesof(X ,...,X ). i=1 i 1 K Suppose that the computational effort to generate these samples on average equals Kγ(1−β) for some β < 1. Further suppose that the correlation between any two random variables X and X for i j i 6= j is ρ ∈ (0,1). Then the variance of PK X equals K(1+ρ(K −1))V[X]. So the simulation i=1 i efficiency equals K2γ(1−β)(1+ρ(K −1))V[X]. We therefore prefer to draw dependent samples whenever 1 β > 1− . (1+ρ(K −1)) 12

Unless 1−β = o (1/K), we will prefer to draw independent samples for any K sufficiently large. K 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 , represent a singlet dimensional market risk factor, and assume X ∼ N(0,1). Let (cid:15) be the idiosyncratic component H k to the return on position k at the horizon, so that the loss on position k is (X +(cid:15) ) per unit of H k exposure. We assume that the (cid:15) are i.i.d. N(0,ν2). To facilitate comparative statics on K, we k scale exposure sizes by 1/K. The exact distribution for portfolio loss Y is N(0,1+ν2/K). We assume that the position-level inner step pricing errors ζ are i.i.d. N(0,η2) per unit of exposure, so that the portfolio pricing k,· error has variance σ2 = η2/K across inner step trials. This implies that the simulated loss variable Y˜ is distributed N(0,1 +ν2/K +σ2/N). Figure 1 shows how the density of the simulated loss distribution varies with the choice of N. For the baseline parameter values ν = 3, η = 10 and K = 100, we observe that the density of Y˜ for N = 64 is a close approximation to the “true” density for Y. Even for N = 16, the error due to inner step simulation appears modest. p Weconsiderourestimatorαˆ ofthelargelossprobabilityα = P(Y > u) = Φ(−u/ 1+ν2/K) L,N forafixedlosslevelu. Theexpectedvalueofαˆ isα = P(Y˜ > u) = Φ (cid:16) −u/ p 1+ν2/K +σ2/N (cid:17) . L,N N Applying Proposition 1, bias in αˆ expands as α −α ≈ θ/N where L,N N σ2 p θ = u(1+ν2/K)−3/2φ(−u/ 1+ν2/K). 2 and where φ is the standard normal density. Figure 2 plots the exact bias α −α and the approx- N imation θ/N as a function of 1/N. The first-order approximation to the bias is quite accurate at modest values of N. At N = 8, the relative error of the approximation is roughly 6%, and even at N = 2, the relative error is under 17%. In Figure 3, we plot the exact root mean square error of αˆ as a function of N and with the L,N computational budget held fixed. Here we make the assumption that the fixed cost of the outer step is negligible (i.e., γ ≈ 0), so that Γ ∝ N ·L. We observe that the optimal N is increasing 0 with the budget, but remains quite modest even for the largest budget depicted (N · L = 216). The relationship between the computational budget and optimal N is explored further in Figure 4. We solve for the N that minimizes the (exact) mean square error, and plot N∗ as a function of the budget Γ = N ·L. Again, we see that N∗ is much smaller than L∗ and grows at a slower rate with Γ. For example, when N ·L = 210, we find N∗ is under 6 and L∗ is over 165. Increasing the budget by a factor of 64, we find N∗ roughly quadruples while L∗ increases by a factor of roughly 16. The figure demonstrates the accuracy of the approximation to N∗ given by equation (4). When 13

Figure 1: Density of the loss distribution Gaussian example with parameters: ν =3, η=10, K =100. Figure 2: Exact and approximated bias SolidlineplotsexactbiasforGaussianexample,dashedlineplotsfirst-orderapproximation of Proposition 1. Parameters: ν =3, η=10, K =100 and u=VaR [Y]. 0.01 14

Figure 3: Root Mean Square Error RMSEinGaussianexample. EachlinedepictstherelationshipbetweenN andRMSEfor a fixed computational budget. Parameters: ν =3, η=10, K =100 and u=VaR [Y]. 0.01 Γ = 210, the relative error of the approximation is 8.3%. Increasing the budget to Γ = 216 shrinks the relative error to under 2.5%. TheoptimalN∗ andtheaccuracyofitsapproximationmaydependontheexceedancethreshold u, and not necessarily in a monotonic fashion. This is demonstrated in the top panel of Figure 5. The budget here is large, and we see that the approximation to N∗ is accurate over the entire range of interest (say, for u < VaR [Y]). When the budget is smaller (bottom panel), the accuracy 0.0001 of the approximation is most severely degraded in the tail. In Figure 6, we explore the relationshipbetweenportfolio size K andoptimalN. For simplicity, weassumeherethatthebudgetgrowslinearlywithK. InthebaselinecaseofK = 100,wefindthat N∗ is roughly 23. When we triple the portfolio size (and budget), we find that N∗ falls to under six. If we increase the portfolio size by a factor of 10, we find that N∗ is under two. These results suggest that the large portfolio asymptotics of section 2.5 may pertain to portfolios of realistic size. As Value-at-Risk is the risk-measure most commonly used in practice, we turn briefly to the estimation of VaR. The results of Section 2.2 show that the bias in yˆ vanishes with 1/N. This is α demonstrated in Figure 7 for three values of α. For each line, the y-axis intercept at 1/N = 0 is the unbiasedbenchmark, i.e., VaR [Y]. Thedistanceonthey-axisbetweenVaR [Y˜](onthesolidline) α α for any finite N and the corresponding unbiased benchmark VaR [Y] is the exact bias attributable α to errors in pricing at the horizon. The dashed lines show the approximated VaR based on 1 θ E[Y˜ ] ≈ y + α . dαLe α N f(y ) α 15

Figure 4: Optimal N grows with computational budget Budget is Γ=N ·L. Parameters: ν =3, η=10, K =100 and u=VaR [Y]. 0.01 Figure 5: Optimal N depends on exceedance threshold QuantilesofthedistributionofY markedinbasispoints. BudgetisΓ=N·L. Parameters: ν =3, η=10 and K =100. 16

Figure 6: Optimal N as portfolio size varies Budget is Γ=N ·L=216 for K =100, and grows linearly with K. Parameters: ν =3, η=10, and u=VaR [Y]. 0.01 For the distributional assumptions of our example, we have θ σ2 Φ−1(1−α) α = . p f(y α ) 2 1+ν2/K Clearly, the linearized approximation to the bias is highly accurate in this example. 4 Jackknife estimation Jackknife estimators, introduced over 50 years ago by Quenouille (1956), are commonly applied when bias reduction is desired. We show how jackknife methods can be applied in our setting to the estimation of large loss probabilities. Parallel methods would apply to VaR and ES. We divide the inner loop sample of N draws into I non-overlapping sections (N and I are selected so that N/I is an integer). Section i covers the N/I draws N (i−1) + 1,...,N i. For I I a single outer step trial, represented by ξ, let a denote the sample output from the inner step N 1[Y˜ > u] as proposed in earlier sections. Let Y˜ (−i) be the estimate of Y that is obtained when ξ ξ ξ section i is omitted, and define a (−i) similarly, e.g., a (−I) = 1[Y˜ (−I) > u]. Observe that the N N ξ bias in a (−i) is θ/(N(I −1)/I) plus O (1/N3/2) terms, and furthermore that we can construct N N I different sample outputs this way. 17

Figure 7: Bias of estimated Value-at-Risk Solid lines are exact VaR [Y˜], dashed lines are approximations based on Proposition 2. α Intercept on y-axis is VaR [Y]. Parameters: ν =3, η=10 and K =100. α We now propose the jackknife inner step simulation output I I 1 X I −1 X a† = Ia −(I −1) a (−i) = a + (a −a (−i)). (15) N I N N I N N i=1 i=1 The bias in a† is E[a†]−α = Iα −(I −1)α −α N N(I−1)/I = I(α+θ/N +O (1/N3/2))−(I −1)(α+θ/(N(I −1)/I)+O (1/N3/2))−α N N (cid:18) (cid:19) I I −1 = θ − +O (1/N3/2) = O (1/N3/2). N N N (N(I −1)/I) Thus, the first-order term in the bias is eliminated. Define b (−i) = a −a (−i) for i ≤ I. From the second representation of a† in (15), we see N N N that the variance of the estimator a† is (I −1)2 (I −1)3 V[a†] = V[a ]+ V[b (−I)]+2(I−1)Cov[a ,b (−I)]+ Cov[b (−1),b (−I)], (16) N I N N N I N N where Cov[X ,X ] denotes the covariance between X and X . We now evaluate the dominant 1 2 1 2 term in this as N → ∞. To keep the analysis simple, we assume that I is fixed even as N increases. As we note later, keeping I small has computational benefits. In Appendix A.5, we demonstrate that 18

Proposition 4 V[b (−i)] = O (1/N1/2). N N From this result and the Cauchy-Schwartz inequality, we have p |Cov[a ,b (−I)]| ≤ V[a ]·V[b (−I)] = O (1/N1/4). N N N N N Similarly, Cov[b (−1),b (−I)] is O (1/N1/2). Therefore, from equation (16) for fixed I indepen- N N N dent of N, V[a†] = α(1−α)+O (1/N1/4). N Our jackknife estimator is the average of L samples of a†, i.e., αˆ† = (1/L) PL a† (ξ ). The ‘=1 ‘ ‘ contribution of variance to its mean square error is α(1−α)+O (1/N1/4) V[αˆ†] = N . L The contribution from the squared bias term is O (1/N3). Suppose that the dominant term in N the squared bias has the form π (an expression for π could be obtained by taking one more term N3 in the Taylor series expansion of α in the proof of Proposition 1). Proceeding as in section 2.1, N we assume that the computation effort takes the form Γ = L(Nγ +γ ), and choose N∗and L∗ to 1 0 minimize the mean square error subject to the budget constraint. We find (cid:18) 3π (cid:19)1/4 N∗ = Γ +o (Γ1/4) Γ α(1−α)γ 1 (cid:18) α(1−α) (cid:19)1/4 L∗ = Γ3/4+o (Γ3/4), 3πγ3 Γ 1 so even fewer inner step samples are needed for the jackknife estimator. The mean square error at optimal N∗ equals (cid:18) α(1−α)γ (cid:19)3/4 4π1/4 1 +o (Γ−3/4). Γ 3Γ Thus, the rate of convergence of the mean square error to zero reduces from order Γ−2/3 to a faster Γ−3/4. In most applications of jackknife methods, bias reduction comes at the price of increased variance. This trade-off applies here as well. To illustrate, we return to the Gaussian example of Section 3. For any parametric example, it is useful to re-write the terms in equation (16) for the variance of a† as: V[b (−I)] = V[a ]+V[a (−I)]−2Cov[a ,a (−I)] N N N N N Cov[a ,b (−I)] = V[a ]−Cov[a ,a (−I)] N N N N N Cov[b (−1),b (−I)] = V[a ]−Cov[a ,a (−I)]−Cov[a ,a (−1)]+Cov[a (−1),a (−I)]. N N N N N N N N N By the symmetry of the a (−i), we have N Cov[a ,a (−1)] = Cov[a ,a (−I)] = P(Y˜ > u,Y˜(−I) > u)−α α . (17) N N N N N N(I−1)/I Similarly, we have Cov[a (−1),a (−I)] = P(Y˜(−1) > u,Y˜(−I) > u)−α2 . (18) N N N(I−1)/I 19

For the Gaussian example, let −u U = n p 1+ν2/K +σ2/n so that α = Φ(U ) and α = Φ(U ). The bivariate probabilities in Equations (17) and (18) are ∞ N N s ! 1+ν2/K +σ2/N P(Y˜ > u,Y˜(−I) > u) = Φ U ,U , 2 N N(I−1)/I 1+ν2/K +σ2/(N(I −1)/I)   1+ν2/K + σ2 I−2 P(Y˜ L (−1) > u,Y˜ L (−I) > u) = Φ 2U N(I−1)/I ,U N(I−1)/I , 1+ν2/K + N(I−1 σ ) 2 /I I−1  N(I−1)/I where Φ (z ,z ,ρ) is the bivariate normal cdf for standard normal marginals and correlation ρ. 2 1 2 Holding fixed N = 32, we plot the bias of the jackknife estimator as a function of the number of sections (I) in the top panel of Figure 8. The bias is decreasing in I, but the sensitivity to I is modest in both absolute terms and relative to the bias of the uncorrected estimator. The bias of the uncorrected a is 9.04 basis points, whereas the bias of a† is -0.29 basis points when I = 2 and N -0.15 basis points when I = N = 32. The standard deviation of a is plotted as a function of I in L,N the bottom panel of the figure. We find the standard deviation increases in roughly linear fashion p p with I. For the uncorrected estimator, we find V[a ] = 0.104, whereas V[a†] = 0.128 for I = 2 N p and V[a†] = 0.606 for I = N = 32. Figure 8: Bias and variance of jackknife estimators Parameters: N =32, ν =3, η=10, K =100 and u=VaR [Y]. 0.01 The optimal choice of I for minimizing mean square error will depend on L. The larger is L, the smaller the contribution of variance to MSE, so the larger the optimal I. As a practical matter, 20

we advocate setting I = 2, which eliminates nearly all the bias at very little cost to variance. Setting I = 2 has the further advantage of minimal memory overhead, in that one can estimate 1[Y˜ (−1) > u], 1[Y˜ (−2) > u], and 1[Y˜ > u] for each outer step ξ in a single pass through the inner ξ ξ ξ loop. By contrast, to implement the jackknife with I = N sections, we need to save each of the N inner loop draws in order to calculate each of the N estimates 1[Y˜ (−1) > u],...,1[Y˜ (−N) > u]. ξ ξ 5 Dynamic allocation Through dynamic allocation of workload in the inner step we can further reduce the computational effort in the inner step while increasing the bias by a negligible controlled amount or even (in many cases) decreasing the bias. Consider the estimation of large loss probabilities P(Y > u) for a given u. Bias in our estimated αˆ is increasing with the likelihood that the indicator 1[Y˜(ξ) > u] obtained from the inner step simulation is not equal to the true 1[Y(ξ) > u] for a given realization ξ of the outer step. Say we form a preliminary estimate Y˜n based on the average of a small number n of inner step trials. If this estimate is much smaller or much larger than u, it may be a waste of effort to generate many more samples in the inner simulation step. However, if this average is close to u, it makes sense to generate many more inner step samples in order to increase the probability that the estimated 1[Y˜(ξ) > u] is equal to the true value 1[Y(ξ) > u]. The variance of αˆ is dominated by α(1−α)/L, so dynamic allocation in the inner step will have little effect on the variance. As in Section 4, we develop our dynamic allocation methods for the estimation of large loss probabilities. Similar methods would apply to VaR and ES. Our proposed dynamic allocation (DA) scheme is very simple. For each trial ξ of the outer step, we generate δN inner step trials (selecting δ so that δN is an integer), and let the resulting preliminary loss estimate be denoted Y˜δ. If Y˜δ < u − (cid:15) for some well-chosen (cid:15) > 0 then we terminate the inner step and our sample output is zero.7 Otherwise, we generate a second estimate Y˜1−δ based on an additional (1−δ)N conditionally independent samples. Defining Y˜∗ = δY˜δ +(1−δ)Y˜1−δ, our final estimate Y˜DA is taken as ξ ( Y˜δ, if Y˜δ < u−(cid:15) Y˜DA = ξ Y˜∗, otherwise. The cost of computation per inner loop is proportional to the average number of inner step trials executed, which is N¯DA ≡ N ·(δ+(1−δ)P(Y˜δ > u−(cid:15))) < N, so the DA estimator always has lower computational effort than the static estimator (for any N held fixed). With dynamic allocation, the bias in the estimated exceedance probability is P(Y˜DA > u,Y˜δ > u−(cid:15))−P(Y > u) = P(Y˜∗ > u)−P(Y > u)−P(Y˜∗ > u,Y˜δ ≤ u−(cid:15)) = (α −α)−P(Y˜∗ > u,Y˜δ ≤ u−(cid:15)). (19) N 7We could improve the scheme slightly by incorporating an upper bound stopping rule as well, setting sample output to 1 if, say, tYδ >u+(cid:15). As this event occurs only rarely, the computational savings are generally modest so we do not develop this modification here. 21

Thus, dynamic allocation introduces a negative increment to bias, relative to the static estimator. As we have seen earlier, in typical application the static estimator has a positive bias. The joint probability P(Y˜∗ > u,Y˜δ ≤ u−(cid:15)) can be made arbitrarily small by choosing δ large (near one) and/or choosing a large (cid:15). Therefore, when α −α > 0, there must exist pairs (δ,(cid:15)) such that the N DA estimator has a lower bias than the static estimator. In this case, dynamic allocation reduces both computational effort and mean square error. If α ≤ α then dynamic allocation will increase absolute bias. Even in this general case, we N can derive an upper bound on the increase in absolute bias. As shown in Appendix B, for every δ one can find an (cid:15) so that the resultant increase in bias is within a specified tolerance. The issue then is to select a δ that most reduces the computational effort. For numerical illustration, we return to our Gaussian example of Section 3. The exact expected value of αˆDA in this example is E[αˆDA] = P(Y˜DA > u) s ! −u (cid:15)−u 1+ν2/K +σ2/N = Φ , , . (20) 2 p 1+ν2/K +σ2/N p 1+ν2/K +σ2/(δN) 1+ν2/K +σ2/(δN) As in our baseline examples, we fix ν = 3, η = 10, K = 100 and u = VaR [Y]. For the static 0.01 scheme with N = 32, we find a bias of 9.0 basis points. For the dynamic allocation scheme, say p we choose δ = 1/32 and (cid:15) = V[Y] = 1.044. Observe that the preliminary sample in the inner step will consist of a single trial. We find the average number of inner step trials executed is N¯DA = 6.24, so computational effort is reduced by over 80%. Bias is -0.4 basis points, so also a large improvement over the static scheme. More complex dynamic schemes are possible. For example, instead of a single stopping test at δN, one can have a series of j stopping tests after δ N,δ N,...,δ N inner step trials. The 1 2 j bandwidths would narrow successively, i.e., (cid:15) > (cid:15) > ...(cid:15) . It may also be possible to design 1 2 j adaptive strategies for choosing δ and (cid:15) that will converge to the optimal values as the outer step simulation progresses. 6 Conclusion We have shown that nested simulation of loss distributions poses a much less formidable computational obstacle than it might initially appear. The essential intuition is similar to the intuition for diversification in the long-established literature on portfolio choice. In the context of a large, well-diversified portfolio, risk-averse investors need not avoid high-variance investments so long as most of the variance is idiosyncratic to the position. In our risk-measurement problem, we see that large errors in pricing at the model horizon can be tolerated so long as the errors are zero mean and mainly idiosyncratic. In the aggregate, such errors have modest impact on our estimated loss distribution. More formally, we are able to quantify that impact in terms of bias and variance of theresultingestimator, andallocateworkloadinthesimulationalgorithmtominimizemeansquare error. Simple extensions of the basic nested algorithm can eliminate much of the bias at modest cost. Our results suggest that current practice is misguided. In order to avoid the “inner step” simulation for repricing at the model horizon, practitioners typically rely on overly stylized pricing models with closed-form solution. Unlike the simulation pricing error in our nested approach, 22

the pricing errors that arise due to model misspecification are difficult to quantify, need not be zero mean, and are likely to be correlated across positions in the portfolio. At the portfolio level, therefore, the error in estimates of Value-at-Risk (or other quantities of interest) cannot readily be bounded and does not vanish asymptotically. Our results imply that practitioners should retain the best pricing models that are available, regardless of their computational tractability. A single trial of a simulation algorithm for the preferred model will often be less costly than a single call to the stylized pricing function,8 so running a nested simulation with a small number of trials in the inner step may be comparable in computational burden to the traditional approach. Despite the high likelihood of grotesque pricing errors at the instrument level, the impact on estimated VaR is small. In the limit of an asymptotically large, fine-grained portfolio, even a single inner step trial per instrument is sufficient to obtain exact pricing at the portfolio level. Our methods have application to other problems in finance. Nested simulation may arise in pricing options on complex derivatives (e.g., a European call option on a CDO tranche). When parameters in an option pricing model are estimated with uncertainty, nested simulation may also be needed to determine confidence intervals on model prices for thinly-traded complex derivatives. Similar problems arise in the rating of CDOs and other structured debt instruments when model parameters are subject to uncertainty. These applications will be developed in future work. A Proofs In this section we prove Lemma 1 and Proposition 1. We then derive the expression for the mean square error of Value-at-Risk discussed in Section 2.2. We then provide a proof of Proposition 3. In our notation, we suppress the dependence on ξ and N wherever it is visually convenient. A.1 Proof of Lemma 1 Note that √ Z Z y−z/ N P(Y˜ ≤ y) = g (v,z)dvdz. N < −∞ Differentiating both sides by y, noting that in the RHS the interchange of integral and the derivative follows from Assumption 1 (Lang, 1968, p. 249), we get Z √ f˜ (y) = g (y−z/ N,z)dz. N N < Therefore, Z (cid:16) √ (cid:17) f˜ (y )−f(y) = g (y −z/ N,z)−g (y,z) dz. N N N N N < This equals Z Z 1 (y −y) g (y˘ ,z)dz− √ zg (y˘ ,z)dz, N N N N N N < < √ where y˘ lies between y −z/ N and y. Clearly, due to Assumption 1, both these terms converge N N to zero. Similarly we see that f˜0 (y ) → f0(y) as N → ∞. N N 8Closed-form solutions to even the simplest pricing models may require special functions, e.g., the confluent hypergeometric function. 23

A.2 Proof of Proposition 1 Note that Z Z ∞ α = g (y,z)dydz, N √ N < u−z/ N so that Z Z u α −α = g (y,z)dydz. (21) N √ N < u−z/ N Consider the Taylor’s series expansion of the pdf in (21), ∂ (y−u)2 ∂2 g (y,z) = g (u,z)+(y−u) g (u,z)+ g (u ,z), N N ∂y N 2 ∂y2 N y where u denotes an appropriate number between u and y. From this and Assumption 1, it follows y that Z z Z z2 ∂ α −α = √ g (u,z)dz− g (u,z)dz+O (1/N3/2). N N N 2N ∂y N N < < Note that the first term above equals f(u) √ E[Z˜ |Y = u]. N N This term equals zero since E[Z˜ |Y = u] = E[E[Z˜ (ξ)|Y(ξ) = u,ξ]] = 0. N N The second term can be re-written as: Z Z ∂ d z2 g (u,z)dz = z2g (u,z)dz ∂y N du N < < d = f(u)E[Z˜2 |Y = u] du N d = f(u)E[E[Z˜2 |ξ]|Y = u] du N d = f(u)E[σ2|Y = u], du ξ which completes the proof. A.3 Mean square error for Value-at-Risk We derive the bias and variance approximations of Proposition 2. Our treatment in places is heuristical to avoid lengthy technical issues. We first develop an asymptotic expression for the bias E[Y˜ ]−y . To evaluate E[Y˜ ], it is dαLe α dαLe useful to consider the order statistics U ≥ ... ≥ U of L random variables uniformly distributed [1] [L] overtheunitinterval. Withinthisappendixonly,letG˜ denotetheinverseofthecdfofY˜. Observe N that Y˜ has the same distribution as G˜ (U ). dαLe N dαLe 24

Consider the Taylor’s series expansion G˜ (U ) = G˜ (E[U ])+(U −E[U ])G˜0 (E[U ]) N dαLe N dαLe dαLe dαLe N dαLe 1 + (U −E[U ])2G˜00(E[U ])+ remainder terms. (22) 2 dαLe dαLe N dαLe Taking expectations E[Y˜ ] = E[G˜ (U )] dαLe N dαLe 1 = G˜ (E[U ])+ V[U ]G˜00(E[U ])+E[remainder terms]. (23) N dαLe 2 dαLe N dαLe It is well known (see, e.g., David, 1981) that E[U ] = 1−dαLe/(L+1) = 1−α+O (1/L) (24) dαLe L α(1−α) V[U ] = +O (1/L2). (25) dαLe L+2 L Through differentiation of the relation F˜ (G˜ (u)) = u, it can be seen that N N f˜0 (G˜ (u)) G˜00(u) = − N N . N f˜ (G˜ (u))3 N N Letting yˆ = G˜ (E[U ]), it can be seen that α N dαLe α(1−α) f˜0 (yˆ ) E[Y˜ ] = yˆ − N α +O (1/L2)(constant+o (1)). (26) dαLe α 2(L+2)f˜ (yˆ )3 L N N α Here we omit the lengthy discussion on the technical assumptions needed to ensure that the expectation of the remainder terms in (23) has the form O (1/L2)(constant+o (1)). L N Let y˜ = G˜ (1−α) denote the α quantile corresponding to the random variable Y˜. If we can α N show that yˆ = y˜ +O (1/L) (27) α α L and that θ α y˜ −y = +o (1/N), (28) α α N Nf(y ) α where θ = −Θ(y ) = −1 d f(u)E[σ2|Y = u]| , then α α 2 du ξ u=yα θ α yˆ = y + +O (1/L)+o (1/N). (29) α α L N Nf(y ) α To see that α(1−α) f˜ N 0 (yˆα) is O (1/L)(1 +o (1)), note that due to Assumption 1, f˜ (·) and 2(L+2)f˜ N(yˆα)3 L N N f˜0 (·) have bounded derivatives, which implies N f˜ (yˆ ) = f˜ (y )+O (1/N)+O (1/L) = f(y )+o (1)+O (1/L). N α N α N L α N L Similarly, f˜0 (yˆ ) = f0(y )+o (1)+O (1/L). N α α N L 25

In particular, from this, (26) and (29), it follows that the bias θ E[Y˜ ] = y + α +O (1/L)+O (1/L)o (1)+o (1/N). dαLe α Nf(y ) L L N N α Equation (27) follows since equation (24) implies that yˆ = G˜ (1−α)+O (1/L) = y˜ +O (1/L) α N L α L for N sufficiently large. We now show (28). Using the Taylor’s series expansion, α = P(Y˜ > y˜ ) = P(Y˜ > y )−(y˜ −y )f˜ (y )+(y˜ −y )2f˜0 (y˘ ) α α α α N α α α N α where y˘ lies between y˜ and y . From Assumption 1 it follows that f˜0 (y) is uniformly bounded α α α N for all y. By Proposition 1, we have P(Y˜ > y ) = P(Y > y )+θ /N +O (1/N3/2). α α α N Since P(Y > y ) = α, it follows that α θ α y˜ −y = +o (1/N). α α Nf˜ (y ) N N α Since, f˜ (y ) → f(y ), (28) follows. N α α The expression for variance can be determined by subtracting (26) from (22) and squaring the difference to get α(1−α) V[Y˜ ] = +O (1/L2)+O (1/L)·o (1). dαLe (L+2)f˜ (y˜ )2 L L N N α Again, from Lemma 1, f˜ (y˜ ) → f(y ) as N → ∞. Therefore, N α α α(1−α) V[Y˜ ] = +O (1/L2)+O (1/L)·o (1). dαLe (L+2)f(y )2 L L N α A.4 Proof of Proposition 3 The Υ(·) and Υ (·) can be written as N Z Z ∞ √ Υ (u) = E[Y˜ ·1[Y˜ > u]] = (y+z/ N)g (y,z)dydz (30) N √ N < u−z/ N Z Z ∞ Υ(u) = E[Y ·1[Y > u]] = yg (y,z)dydz. (31) N < u As E[Z˜ |Y = y] = 0 for each y, we have N Z Z ∞ Z ∞ zg (y,z)dydz = E[Z˜ |Y = y]f(y)dy = 0. N N < u u 26

Substituting into (30), we have Z Z u Z √ Z u Υ (u)−Υ(u) = yg (y,z)dydz+ z/ N g (y,z)dydz. (32) N √ N √ N < u−z/ N < u−z/ N We apply a Taylor series expansion to the first term in (32): Z Z u yg (y,z)dydz √ N < u−z/ N Z Z u (cid:18) ∂ (y−u)2 ∂2 (cid:19) = y g (u,z)+(y−u) g (u,z)+ g (u ,z) dydz √ N ∂y N 2 ∂y2 N y < u−z/ N Z (cid:18) uz z2 (cid:19) Z (cid:18) uz2 z3 (cid:19) ∂ = − g (u,z)dz− − g (u,z)dz+o (1/N3/2) N1/2 2N N 2N 3N3/2 ∂y N N < < 1 u d (cid:16) (cid:17) = − f(u)E[Z˜2|Y = u]− f(u)E[Z˜2|Y = u] +O (1/N3/2). 2N N 2N du N N Here, in the first equality u denotes an appropriate number between u and y, in the third equality y R we use the fact that zg (u,z)dz = 0. < N For the second term in (32), Z √ Z u z/ N g (y,z)dydz √ N < u−z/ N Z √ Z u (cid:18) ∂ (y−u)2 ∂2 (cid:19) = z/ N g (u,z)+(y−u) g (u,z)+ g (u ,z) dydz √ N ∂y N 2 ∂y2 N y < u−z/ N Z z2 Z z3 ∂ = g (u,z)dz− g (u,z)dz+o (1/N3/2) N N 2N3/2∂y N N < < 1 = f(u)E[Z˜2|Y = u]+O (1/N3/2). N N N Proposition 3 follows by noting that E[Z˜2|Y = u] = E[E[Z˜2|Y,ξ]|Y = u] = E[σ2|Y = u]. N N ξ A.5 Proof of Proposition 4 We decompose the variance of b (−i) as E[b (−i)2]−E[b (−i)]2. The expectation of b (−i) is N N N N E[b (−i)] = E[a ]−E[a (−i)] N N N θ θ = α+ u +O (1/N3/2)−(α+ u +O (1/N3/2)) N N N N(I −1)/I −θ = u +O (1/N3/2), N N(I −1) which implies that E[b (−i)]2 = O (1/N2). N N We now argue that E[b (−I)2] is O (1/N1/2), ignoring some minor technical issues. From N N expanding the square and taking expectations, observe that E[b (−I)2] = P(Y +Z¯N > u,Y +Z¯N(I−1)/I ≤ u)+P(Y +Z¯N ≤ u,Y +Z¯N(I−1)/I > u). (33) N 27

We now show that these two probabilities are equal in their dominant term and are O(1/N1/2). Some notation is needed for this purpose. Let (cid:18) I (cid:19)1/2 X N A ≡ Z N i N i=N(I−1)/I (cid:18) I (cid:19)1/2N(I X −1)/I B ≡ Z . N i N(I −1) i=1 We assume that both A and B have limiting distributions with finite expectations. Then, N N 1 (cid:16) (cid:17) Z¯N = (1/I)1/2A +(I −1/I)1/2B . N N N1/2 Let g† denote the joint pdf of (Y,A ,B ). Then the first term in (33) can be written as N N N P(Y +Z¯N > u,Y +Z¯N(I−1)/I ≤ u) Z Z u−b/(N(I−1)/I)1/2 = g†(y,a,b)dy da db. N a≥b/(I−1)1/2 u−a/(NI)1/2−b/(NI/(I−1))1/2 Taking a Taylor series expansion of g†(y,a,b) at the first argument set to u and under assumptions N similar to Assumption 1, this can be seen to equal 1 f(u)E[(A −B /(I −1)1/2)·1[A −B /(I −1)1/2 > 0]|Y = u]+O (1/N). N N N N N (NI)1/2 Since A and B have limiting distributions with finite expectations, the expectation term con- N N verges to a constant as N → ∞, which implies that 1 f(u)E[(A −B /(I −1)1/2)·1[A −B /(I −1)1/2 > 0]|Y = u] = O (1/N1/2). N N N N N (NI)1/2 We proceed in exactly the same fashion for the second term in (33), and find that this term is also O (1/N1/2). N B Bounding the bias under dynamic allocation We can bound the absolute bias under dynamic allocation via the triangle inequality: |P(Y˜∗ > u,Y˜δ > u−(cid:15))−P(Y > u)| = |P(Y˜∗ > u)−P(Y > u)−P(Y˜∗ > u,Y˜δ ≤ u−(cid:15))| ≤ |P(Y˜∗ > u)−P(Y > u)|+P(Y˜∗ > u,Y˜δ ≤ u−(cid:15)) so the increase in the absolute bias is no greater than P(Y˜∗ > u,Y˜δ ≤ u−(cid:15)). Decomposing this as P(Y˜∗ > u,Y˜δ ≤ u−(cid:15)) = P(Y˜∗ > u,Y˜δ ≤ u−(cid:15),Y > u)+P(Y˜∗ > u,Y˜δ ≤ u−(cid:15),Y ≤ u), (34) 28

we can upper bound each of these terms. Let Z¯δ and Z¯1−δ be the average portfolio pricing errors for the preliminary sample of δN and the continuation sample of (1−δ)N inner step trials. For the first term in equation (34), we have the inequality P(Y˜∗ > u,Y˜δ ≤ u−(cid:15),Y > u) ≤ P(Z¯δ ≤ −(cid:15)). For the second term, we have (cid:18) (cid:19) (cid:18) (cid:19) δ δ P(Y˜∗ > u,Y˜δ ≤ u−(cid:15),Y ≤ u) ≤ P Y˜1−δ > u+ (cid:15),Y ≤ u ≤ P Z¯1−δ > (cid:15) . 1−δ 1−δ Hence, an upper bound on the additional bias is (cid:18) (cid:19) δ P(Z¯δ ≤ −(cid:15))+P Z¯1−δ > (cid:15) . (35) 1−δ From here we can take two different approaches: 1. Hoeffding’sinequality(Hoeffding,1963)canbeusedtodevelopexactboundsifthezeromean noise associated with each position is upper as well as lower bounded. This is generally true of put options, credit derivatives, and many other types of derivatives. Let a < ζ < b k k,· k denote the bounds on position k. Then, ! 2δN(cid:15)2 P(Z¯δ ≤ −(cid:15)) ≤ exp − PK (b −a )2 i=1 i i ! (cid:18) δ (cid:19) 2δ2N(cid:15)2 P Z¯1−δ > (cid:15) ≤ exp − . 1−δ (1−δ) PK (b −a )2 i=1 i i 2. As Z is the sum of zero mean errors from K positions, we can appeal to the central limit i theorem and take each Z (ξ) as approximately normal in distribution with mean zero and i variance σ2.9 Then, ξ ! −(cid:15) P(Z¯δ ≤ −(cid:15)) = Φ √ σ / δN ξ ! (cid:18) (cid:19) δ −(cid:15)δ/(1−δ) P Z¯1−δ > (cid:15) = Φ . p 1−δ σ / (1−δ)N ξ Thus, given an acceptable threshold on the increase in bias, for every δ one can find an (cid:15) so that the resultant increase in bias is controlled. The issue then, is to select a δ that most reduces the computational effort. We illustrate this computation using our simple Gaussian example with baseline parameter √ values K = 100, ν = 3, and η = 10. This gives us σ = η/ K = 1 as the standard deviation of the p Z . Let u = VaR [Y] = 1+ν2/KΦ−1(0.99), so that α = 0.01. Say we fix N = 30, δ = 1/3 and i 0.01 9One can, if necessary, batch the sum of a few Z ’s to improve the approximation to the normal distribution. If i σ2 isunknown,onecanestimateitfromthesamplevariancewithintheinnerstep. Inthiscase,thestandardizedZ ξ i have Student t distribution. 29

(cid:15) = 2. The upper bound in equation (35) for the increase in bias for αˆDA over the static estimator αˆ is under 4×10−6. The probability of stopping with the preliminary sample is L,N ! u−(cid:15) P(Y˜δ ≤ u−(cid:15)) = Φ ≈ 0.653. p 1+ν2/K +σ2/(δN) The average number of inner step trials executed works out to N¯DA ≈ 0.565N. For a negligible upper bound on the increase in bias, we cut computational effort nearly in half. References Carlo Acerbi and Dirk Tasche. On the coherence of expected shortfall. Journal of Banking and Finance, 26(7):1487–1503, 2002. Philippe Artzner, Freddy Delbaen, Jean-Marc Eber, and David Heath. Coherent measures of risk. Mathematical Finance, 9(3):203–228, 1999. Herbert Aron David. Order Statistics. John Wiley & Sons, second edition, 1981. Ben De Prisco, Ian Iscoe, Yijun Jiang, and Helmut Mausser. Compound scenarios: An efficient framework for integrated market-credit risk. Research paper series, Algorithmics, May 2007. Paul Glasserman. Monte Carlo Methods in Financial Engineering. Springer-Verlag, New York, 2004. Peter W. Glynn and Ward Whitt. The asymptotic efficiency of simulation estimators. Operations Research, 40:505–520, 1992. Michael Gordy and Sandeep Juneja. Efficient simulation for risk measurement in a portfolio of cdos. In L. Felipe Perrone, Barry G. Lawson, Jason Liu, and Frederick P. Wieland, editors, Proceedings of the 2006 Winter Simulation Conference, Piscataway, NJ, 2006. IEEE Press. Michael B. Gordy. Granularity adjustment in portfolio credit risk measurement. In Giorgio P. Szeg¨o, editor, Risk Measures for the 21st Century. John Wiley & Sons, 2004. C. Gouri´eroux, J.P. Laurent, and O. Scaillet. Sensitivity analysis of values at risk. Journal of Empirical Finance, 7:225–245, 2000. Wassily Hoeffding. Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association, 58:13–30, March 1963. Serge Lang. Analysis I. Addison-Wesley, 1968. Shing-Hoi Lee. Monte Carlo Computation of Conditional Expectation Quantiles. PhD thesis, Stanford University, October 1998. 30

Francis A. Longstaff and Eduardo S. Schwartz. Valuing American options using simulation: A simple least-squares approach. Review of Financial Studies, 14:113–147, 2001. Richard Martin and Dirk Tasche. Shortfall: a tail of two parts. Risk, 20(2):84–89, February 2007. Richard Martin and Tom Wilde. Unsystematic credit risk. Risk, 15(11):123–128, November 2002. Maurice Quenouille. Notes on bias and estimation. Biometrika, 43:353–360, 1956. R.TyrrellRockafellarandStanislavUryasev.Conditionalvalue-at-riskforgenerallossdistributions. Journal of Banking and Finance, 26(7):1443–1471, 2002. Michael Rothschild and Joseph E. Stiglitz. Increasing risk I: A definition. Journal of Economic Theory, 2(3):225–243, 1970. 31