Efficient Monte Carlo Counterparty Credit Risk Pricing and Measurement

Finance and Economics Discussion Series Divisions of Research & Statistics and Monetary Affairs Federal Reserve Board, Washington, D.C. Efficient Monte Carlo Counterparty Credit Risk Pricing and Measurement Samim Ghamami and Bo Zhang 2014-114 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.

Efficient Monte Carlo Counterparty Credit Risk Pricing and Measurement∗ Samim Ghamami† and Bo Zhang‡ December 17, 2014 Abstract Counterpartycreditrisk(CCR),akeydriverofthe2007-08creditcrisis, hasbecomeone ofthemainfocusesofthemajorglobalandU.S.regulatorystandards. Financialinstitutions invest large amounts of resources employing Monte Carlo simulation to measure and price their counterparty credit risk. We develop efficient Monte Carlo CCR estimation frameworks by focusing on the most widely used and regulatory-driven CCR measures: expected positiveexposure(EPE),creditvalueadjustment(CVA),andeffectiveexpectedpositiveexposure (EEPE). Our numerical examples illustrate that our proposed efficient Monte Carlo estimators outperform the existing crude estimators of these CCR measures substantially in terms of mean square error (MSE). We also demonstrate that the two widely used sampling methods, the so-called Path Dependent Simulation (PDS) and Direct Jump to Simulation date (DJS), are not equivalent in that they lead to Monte Carlo CCR estimators which are drastically different in terms of their MSE. Keywords: RiskManagement,CounterpartyCreditRisk,OTCDerivativesMarket,Credit Value Adjustment, Efficient Monte Carlo Simulation, Basel II-III 1 Introduction and a Summary of Important CCR Measures Counterpartycreditrisk(CCR)istheriskthatapartytoaderivativecontractmaydefaultprior to the expiration of the contract and fail to make the required contractual payments, (see [5] for the basic CCR definitions). Counterparty credit risk has been widely considered as one of the keydriversofthe2007-08creditcrisis, andithasbecomeoneofmainfocusesofthemajorglobal ∗The opinions expressed in this paper are those of the authors and do not necessarily reflect the views of the Board of Governors or its staff. Samim Ghamami is grateful to Robert Anderson, Lisa Goldberg, and Travis Nesmith for their comments. †FederalReserveBoard,Washington,D.C.,email: samim.ghamami@frb.gov,andCenterforRiskManagement Research, University of California, Berkeley, CA, email: samim ghamami@berkeley.edu. ‡BusinessSolutionsandMathematicalSciencesDepartment, IBMThomasJ.WatsonResearchCenter, Yorktown Heights, NY, 10598, email:bozhang@gatech.edu. 1

andU.S.regulatoryframeworks(BaselIII1andtheDodd-FrankActof2009-10; see, forinstance, [3]). It is well known that pricing and measuring counterparty credit risk is computationally extremely intensive; financial institutions (derivative dealers) invest large amounts of resources developing and maintaining Monte Carlo simulation “engines” to manage their counterparty risk, (see [24], [18], and [5]). While various aspects of counterparty credit risk have been subject of extensive research post 2007-08 financial crisis, statistical efficiency of the CCR estimators has received no attention in the literature. In this paper we develop efficient Monte Carlo frameworks for pricing and measuring counterparty risk. More specifically, we focus on efficient Monte Carlo estimation of the most widely used and regulatory-driven CCR measures, expected positive exposure (EPE), credit value adjustment (CVA), and effective EPE (EEPE), as defined below. Efficiency criteria under consideration are variance, bias, and computing time of the Monte Carlo estimators. Our proposed Monte Carlo estimators of EPE, CVA, and EEPE outperform the existing crude estimators of these CCR measures substantially in terms of mean square error (MSE). To the best of our knowledge, this paper is the first to develop efficient Monte Carlo counterparty credit risk estimators. Currently, CVA is a CCR measure that is only applied to bilateral derivatives transactions. However, EPE and EEPE are CCR measures applicable to both bilateral derivatives transactions and centrally-cleared derivatives transactions. Specifically, Basel Committee on Banking Supervision (BCBS) has devised regulatory capital charges on clearing member banks against their central counterparty credit risk; EPE and EEPE are components of these capital charges to be estimated via Monte Carlo simulation by large dealer banks. Counterparty credit exposure [5], denoted by V, of a financial institution against one of its counterparties, is the larger of zero and the market value of the portfolio of derivatives contracts the financial institution holds with this counterparty. To effectively introduce our efficient Monte Carlo procedures, we consider credit exposures in the absence of the commonly used risk mitigants, collateral and netting agreements. This simple setting facilitates the communication of our main results. EPE is a widely used counterparty credit risk measure for regulatory and economic capital calculations, (see Chapters 2 and 11 of [18]). It is defined as follows, (cid:90) T EPE ≡ E[V ]dt, (1) t 0 where E[V ] is the expected value of the (credit) exposure at time t ≥ 0, and T > 0 denotes the t time to maturity of the longest transaction in the derivatives portfolio. Effective EPE (EEPE), another widely used regulatory and economic capital-related counterparty risk measure [18] is defined as follows in the CCR literature: n (cid:88) EEPE ≡ max E[V ]∆ . (2) dst j i 1≤j≤i i=1 1Basel III is a global regulatory standard on bank capital adequacy, stress testing and market liquidity risk agreed upon by the members of the Basel Committee on Banking Supervision in 2010-11, and scheduled to be introduced from 2013 until 2018. 2

This definition is based on a discrete time grid, 0 ≡ t < t < ... < t ≡ T with ∆ = t −t , 0 1 n i i i−1 i = 1,...,n. We prefer and propose the following continuous version of EEPE: (cid:90) T EEPE ≡ max E[V ]dt, (3) u 0 0≤u≤t which is consistent with the definition of EPE and has the advantage of not requiring an a priori specification of a discrete time grid. Our results in Section 5 apply to EEPE as well as EEPE . dst EEPE is the “conservative” version of EPE that accounts for roll-over risk. Roll-over risk refers to the following scenario. Expiration of some of the short-term trades in the derivatives portfolio before T would decrease some of the E[V ] and so EPE. However, it is likely that these i short-term trades are replaced by new ones. When these replacements are not captured by the Monte Carlo CCR “engine”, EPE is underestimated, (see [24]). CVA, which is the difference between the risk free portfolio value and the true counterparty defaultriskyportfoliovalue, (see[23]), hasbecomeoneofthemostimportantCCRmeasuresfor derivativesdealersandotherlargefinancialinstitutions. CVAisalsooneofthemaincomponents of the Basel III’s counterparty credit risk capital framework. Let τ, a positive random variable, denote the default time of the counterparty. It can be shown that CVA, the price of the counterparty credit risk, is equal to the risk neutral expected discounted loss, i.e., CVA ≡ E[(1−R)D V 1{τ ≤ T}], (4) τ τ where 1{A} is the indicator of the event A, D = B /B is the stochastic discount factor at time t 0 t t, B is the value of the money market account at time t, and R is the financial institution’s t recovery rate, (see, for instance, Chapter 7 of [18] for a derivation of this formula). Hereafter we suppress the dependence of the CVA on the recovery rate, R. When V and τ are assumed to be independent, we refer to CVA as independent CVA. Let F denote the cumulative distribution function of τ. Independent CVA can be written as follows, (cid:90) T CVA ≡ E[D V 1{τ ≤ T}] = E[D V ]dF , (5) I τ τ t t t 0 where the last equality follows from conditioning on τ, the independence of V and τ, and the independence of D and τ. We focus on efficient Monte Carlo estimation of independent CVA in this paper. 2 EPE,EEPE,andindependentCVAareestimatedbasedontheReimannsumapproximation of the integrals in (1), (3), and (5) and Monte Carlo estimation of expected exposures, E[V ], t and expected discounted exposures, E[D V ]. t t Section 2 summarizes the common features of the Monte Carlo CCR framework widely used by financial institutions and introduces the notion of Marginal Matching, which enables us to defineanddifferentiatethetwowidelyusedCCRsamplingmethods,Path Dependent Simulation 2Wrong (right) way risk are referred to as cases where credit exposures are negatively (positively) correlated with the credit quality of the counterparty, (see [12], [5], and [19]). 3

(PDS) and Direct Jump to Simulation date (DJS). These two terms were first introduced by Pykhtin and Zhu in 2006 [24]. Practitioners do not usually evaluate and compare the statistical efficiency and computing time associated with the PDS and DJS-based Monte Carlo estimators of CCR measures. A recurring theme of Sections 3 through 5 of this paper is to illustrate that PDS and DJS-based CCR estimators have drastically different MSE.3 Section 3 introduces an efficient Monte Carlo framework for estimating EPE, which also directly applies to efficient CVA estimation. Using our results in Section 3, we summarize our proposed Monte Carlo I framework for efficient estimation of CVA in Section 4. Using our results in Section 3, Section I 5 considers efficient Monte Carlo estimation of EEPE. Our numerical examples indicate that employing our Monte Carlo CCR schemes leads to substantial MSE reduction. We would like to emphasize that Sections 4 and 5 should not be read independently. The main components of the proposed efficient CCR framework are developed in Section 3. As will be seen in the sequel, this is because EPE and CVA are both weighted sums of expected exposures, and the I ideas developed for efficient EPE and CVA estimation have implications for efficient EEPE I estimation. Appendix H discusses the scope of our study by considering cases that have not been explicitly formulated in the main body of the paper. 2 Monte Carlo Counterparty Credit Risk Estimation Contract level credit exposure at time t > 0 is the maximum of the contract’s market value and zero, max{C ,0}, where C denotes the time-t value of the derivative contract. Consider t t a financial institution that holds a portfolio of k derivative contracts with its counterparty. Counterparty level credit exposure is k (cid:88) V = max{Ci,0}, (6) t t i=1 where Ci denotes the time-t value of the i’th derivative contract in the derivatives portfolio. t When risk mitigants are employed, V is defined differently. For instance, in the presence of t netting agreements, credit exposure becomes, (see [23]), k (cid:88) V = max{ Ci,0}. (7) t t i=1 A typical Monte Carlo counterparty risk engine of a derivatives dealer estimates various types of CCR measures based on sampling from the credit exposure process on a time grid, 0 < t < ... < t = T, where T denotes the maturity of the longest transaction in a portfolio of 1 n derivatives and t ,...,t are sometimes referred to as valuation points. Set V ≡ V . 1 n i ti 3The first version of the paper had used the first edition of the excellent book by Gregory [17] as one of its referencesasatthetimeofthecompletionofourmainresultsthesecondeditionofthebookbyGregory[18]had not been published. While [17] does not discuss statistical efficiency of CCR measure, [18] does. The statistical efficiency discussions of [18] are summarized at the Section 3.1.3’s footnote. 4

SomeoftheCCRmeasuresarestatic inthesensethattheyaredefinedbasedonagivenfixed time point. Expected exposure (EE) at time t , is simply E[V ]. Also, Value at Risk (VaR) type i i of measures for a given valuation point t is referred to as potential future exposure. Derivatives i dealers use Monte Carlo simulation to estimate EE and PFE for all the given valuation points t ,...,t on a frequent basis, (see [18] and [23] for more details). Note that EPE, CVA, and 1 n EEPE, as defined in Section 1, are dynamic in the sense that they depend on the time evolution of the credit exposure process. This paper focuses on developing a Monte Carlo CCR framework that provably improves the efficiency of estimating these three risk measures. PFE and its corresponding dynamic measure Maximum PFE are other important CCR measures (see [18]), whose efficient Monte Carlo estimation is not addressed in this paper. Efficient Monte Carlo quantile estimation is a well-studied topic, (see Chapter 9 of [13] and the references therein). As will be illustrated in the sequel, the mathematical analysis of this paper has benefited from some similarities in the functional forms of EPE, CVA, and EEPE. The functional form of the maximum PFE, i.e., the maximum of the quantiles of the exposure process on a discrete time grid, and its Monte Carlo estimator has a distinct nature. We leave efficient Monte Carlo Maximum PFE estimation for future study. In what follows we first summarize the simulation of the credit exposure process. Then, we introduce the notion of Marginal Matching in sampling from the time evolution of the credit exposure process. 2.1 Simulating the Credit Exposure Process Suppose that credit exposure is a stochastic process {V ;t ≥ 0} defined on a given filtered t probability space (Ω,F,(F ) ,P). Given (6) and (7), V can be viewed as a function of t 0≤t≤∞ t the stochastic processes that drive the values of the derivative contracts, C1,...,Ck. In risk t t management, these underlying stochastic processes are usually referred to as risk factors, e.g., interest rates, commodity prices, and equity prices. To generate a Monte Carlo realization of V , for a fixed t > 0, first, the underlying risk factors should be sampled from up to time t > 0. t Next, given the Monte Carlo realization of the risk factors up to time t > 0, the derivative contracts should be valued. This two-step procedure leads to a single Monte Carlo realization of V . It is a risk management common practice to use the physical probability measure in the t first step and the risk-neutral measure in the second. This applies to Monte Carlo estimation of EPE and EEPE. However, since CVA is usually viewed as the market price of counterparty credit risk, risk-neutral measure is usually used in both steps. Depending on the complexity of the payoff function of the derivative contracts, the valuation step could take straightforward Black-Scholes-type analytical calculations, or it could demand approximations that depending on the desired level of accuracy might be computationally intensive. These approximations could also involve Monte Carlo simulation: Nested Monte Carlo refers to the use of a second layer of Monte Carlo simulation in the valuation step of the above procedure, (see [16]), and regression-based Monte Carlo (see [4]) uses ideas from regression-based Monte Carlo American option pricing, (see Chapter 8 of [13]). 5

2.2 Marginal Matching Let X = (X ,...,X ) denote a random vector with distribution function F . Let ω ≡ 1 n X X (E[h (X )],...,E[h (X )]) for some functions h ,...,h . And let θ ≡ g(ω ) for a function 1 1 n n 1 n X X g that maps ω from Rn to R. Two simple examples of θ are as follows, X X n (cid:88) E[h(X )] and max{E[h(X )],...,E[h(X )]}, i 1 n i=1 that is θ is defined based on the marginal distribution of (functions of) X ,...,X . Let Y = X 1 n (Y ,...,Y ) denote another random vector with distribution function F such that, 1 n Y X (cid:54)=d Y , X =d Y for all i = 1,...,n, (8) i i where =d denotes “being equal in distribution”. Simply note that since the marginal distributions of X and Y match, θ = θ . Now, suppose that θ is to be estimated with Monte Carlo X Y X simulation. Given (8), samples can be drawn from F or F . Let θˆ and θˆ denote Monte X Y X,m Y,m Carlo estimators of θ based on m simulation runs when samples are drawn from F and F , X X Y respectively. Obviously, θˆ (cid:54)=d θˆ , X,m Y,m and so between θˆ and θˆ , i.e., when deciding on whether to sample from F or F , the X,m Y,m X Y estimator with a lower mean square error (MSE) should be chosen. Example: Finite-Dimensional Distributions of Brownian Motion Let {X ;t ≥ 0} t denote a Brownian motion with drift µ and volatility parameter σ. Consider the random vector X = (X ,...,X ) ≡ (X ,...,X ) on the time grid, 0 < t < t < ... < t . That is, following the 1 n t1 tn 1 2 n basic definition of a Brownian motion, X is a multivariate normal random vector with E[X ] = ti µt and Var(X ) = σ2t , and cov(X ,X ) = t > 0 for t < t . Now, let Y = (Y ,...,Y ) denote i ti i ti tj i i j 1 n a multivariate normal random vector whose marginal distributions match that of X but with cov(Y ,Y ) = 0, i.e., components of Y are independent normals. i j Stochastic Models of the Risk Factors Let {R ;t ≥ 0}, representing the dynamics t of a risk factor, denote a stochastic process defined on a given filtered probability space, (Ω,F,(F ) ,P). In this paper we assume that {R ;t ≥ 0} is in the following class:4 t 0≤t≤∞ t a Gauss-Markov process (see Chapter 5 of [20] or Chapter 3 of [13]) specified by dR = (g +h R )dt+σ dB , t t t t t t with g,h, and σ all deterministic functions of time and B a standard one-dimensional Brownian motion, a geometric Brownian motion (GBM), 4This assumption is only used in the proof of Proposition 1 and Proposition 2. 6

dR = µR dt+σdB , t t t with given constants µ and σ, or a square-root diffusion specified as (cid:112) dR = α(b−R )dt+σ R dB , t t t t in which α and b are positive. Many of the widely used continuous time stochastic processes in finance and economics are in this class. Consider the finite dimensional distribution of R on a time grid, t ,...,t and set R ≡ R . Suppose that R = (R ,...,R ) can be sampled from 1 n i ti 1 n exactly in the sense that the distribution of the simulated R is precisely that of the R process at times t ,...,t ; examples are Brownian motion, Ornstein-Uhlenbeck processes, GBM, and the 1 n square-root diffusion specified above whose simulations involve generating positively correlated normal random variables. Let R˜ = (R˜ ,...,R˜ ) denote a random vector for which R˜ (cid:54)=d R but 1 n R˜ =d R for all i = 1,...,n and cov(R˜ ,R˜ ) = 0 for all i (cid:54)= j. That is, simulation of R˜ ,...,R˜ i i i j 1 n can be done by generating n uncorrelated or simply independent normal random variables. PDS Sampling versus DJS Sampling In the CCR literature when counterparty risk measuresareestimatedbasedonsamplingfromthefinite-dimensionaldistributionsoftheunderlying risk factors, the sampling is referred to as Path Dependent Simulation (PDS sampling). When the notion of marginal matching is used, the sampling is referred to as Direct Jump to Simulation date (DJS). For instance, in the Brownian motion example above, sampling from X and Y when estimating θ -type estimands are referred to as PDS and DJS sampling, respectively. In X Monte Carlo estimation of CCR measures, PDS and DJS sampling have been widely considered equivalent (see [24]). One of the main contributions of this paper is to differentiate DJS and PDS in terms of the mean square error of the estimators of EPE, CVA, and EEPE. 3 Efficient Monte Carlo Estimation of EPE In this section we consider efficient Monte Carlo estimation of EPE, (cid:90) T EPE = E[V ]dt, t 0 where V denotes the credit exposure process, and T > 0 represents the expiration time of the longest maturity derivative contract in an OTC derivatives portfolio. Consider a time grid, 0 ≡ t < t < ... < t ≡ T, with a fixed n. Set ∆ ≡ t −t and V ≡ V , i = 1,...,n. Let 0 1 n i i i−1 i ti θˆ denote a class of Monte Carlo estimators of EPE defined as follows, b,m,n,k n (cid:88) θˆ ≡ V¯∆ , b,m,n,k i i i=1 7

where V¯ ≡ (cid:80)m V /m and V ,...,V represent the m simulation samples at valuation point i j=1 ij i1 im t . The subscript b refers to the biased nature of the estimators, and the subscript k could take i p and d, referring to PDS and DJS based simulation of the credit exposure process, respectively. As mentioned in Section 2.1, simulating the credit exposure process involves sampling from the underlying risk factors. Hereafter, PDS and DJS-based simulations of the credit exposure process refer to the cases where the underlying risk factors are sampled from based on their finite dimensional distributions (PDS sampling) and based on the notion of marginal matching (DJS sampling), respectively. Note that, n (cid:32) n (cid:90) T (cid:33)2 (cid:88) (cid:88) MSE(θˆ ) = Var( V¯∆ )+ E[V¯]∆ − E[V ]dt . b,m,n,k i i i i t 0 i=1 i=1 We assume that Monte Carlo realizations of V are unbiased estimates of E[V ], i = 1,...,n. This i i implies that the bias part of the MSE of θˆ is not affected by the choice of the sampling b,m,n,k method (PDS or DJS). In Section 3.1, we assume that n, the number of valuation points, is fixed, and we compare the efficiency of θˆ and θˆ in terms of variance and computing b,m,n,p b,m,n,d time both for path independent and path dependent derivatives. Next, we introduce our efficient biased, yet consistent Monte Carlo estimators of EPE. In Section 3.3 we introduce efficient unbiased estimators of EPE. Numerical examples in Section 6.1 indicate that our proposed estimators substantivally outperform the crude estimators of EPE in terms of the mean square error. 3.1 Comparing PDS and DJS-based Estimation of EPE Suppose that the credit exposure process, V, is defined on a given filtered probability space (Ω,F,(F ) ,P), where (F ) denote the filtration generated by the underlying risk t 0≤t≤∞ t 0≤t≤∞ factors. Consider the setting where V denotes the contract level exposure and a financial institution takes a position in a maturity-T derivative contract with its counterparty. Let Π denote T the payoff function of the derivative contract. It is well known from martingale pricing that (cid:20) (cid:21) Π T C = n E |F , (9) t t t n T where n is a numeraire. Transactions between the financial institution and its counterparty for which V = max{C ,0} = C for all 0 < t ≤ T are referred to as unilateral transactions, e.g. the t t t financial institution takes a long position in a call option with its counterparty. Transactions for which V = max{C ,0} =(cid:54) C for some 0 < t ≤ T are referred to as bilateral transactions, e.g. t t t an interest rate swap between the financial institution and its counterparty. ThefollowingsimpleexamplereviewssimulationoftheexposureprocessunderPDSandDJS. Supposethat{S t ;t ≥ 0}isaGBM,S t = S 0 eXt,and{X t ;t ≥ 0}isaBrownianmotionwithdrift (cid:104) (cid:105) µ and volatility σ. Consider a unilateral transaction. Note that V = C = n E ΠT|S ≡ f(S ). t t t nT t t 8

That is, credit exposure is considered as a function of the risk factor.5 Consider the time grid, 0 ≡ t < t < ... < t ≡ T and let V ≡ V . Set 0 1 n i ti n (cid:88) θ ≡ E[V ]∆ . i i i=1 Recall that, n (cid:88) θˆ = V¯∆ , b,m,n,k i i i=1 where V¯ i is the m-simulation-run average of V i1 ,...,V im . With V i = f(S i ) and S i = S 0 eXi, Monte Carlo estimation of θ requires sampling from the multivariate normal random vector, X = (X ,...,X ). Thisistheso-calledPDSsamplingmethod. Analternativesamplingmethod, 1 n usingthenotionofmarginalmatching,istosamplefromthemultivariatenormalrandomvector, Y = (Y ,...,Y ), whose components are uncorrelated but marginal distributions match those of 1 n X. This is the so-called DJS method. To be more specific, in DJS sampling, S is generated i fromtimezero. Thatis, generateY , anormalrandomvariablewithmeanµt andvarianceσ2t , i i i and set S i = S 0 eYi. In PDS sampling, V i ’s are sampled based on generating the sample path of the GBM sequentially at i = 1,...,n. That is, to generate a realization of V , S is generated i i given the previously generated value of S . More specifically, to sample from S generate X˜ i−1 i i and set S i = S i−1 eX˜ i, where X˜ i is a normal random variable with mean µ∆ i and variance σ2∆ i . Note that since for any given t > 0, V t is a function of S t = S 0 eXt, DJS-based simulation of the exposure process implies that cov(V ,V ) = 0 for any i (cid:54)= j, i,j = 1,...,n. i j In what follows we compare the efficiency of θˆ and θˆ in terms of variance and b,m,n,p b,m,n,d computing time for path independent and path dependent derivatives. We consider unilateral and bilateral transactions in both single risk-factor and multi-risk factor settings. That is, we consider two cases: a stylized setting where (F ) is the filtration generated by a single risk t 0≤t≤∞ factor; we also consider the more general multi-risk factor settings. 3.1.1 Path Independent Case The above mentioned example shows that under DJS, cov(V ,V ) = 0 for any 0 < u < t < T. u t Proposition1andProposition2,whoseproofsaregivenintheAppendix,considerthiscovariance function of the contract level credit exposure process under the PDS method for unilateral and bilateral transactions, respectively, and identify conditions under which cov(V ,V ) > 0 for any u t 0 < u < t < T. Condition 2 of Proposition 1 below uses the well known changes of numeraire 5Consider, for instance, the payoff function Π =(S −K)+ of a maturity-T GBM-driven vanilla call option T T withstrikeK. Assumingzeroshort rate, C =E[(S −K)+|S ]=E[(S S −K)+|S ]. Notethatthe function t T t t T−t t f in f(S ) ≡ E[(S S −K)+|S ], which is well-defined for all values of t ≥ 0 given the payoff function Π t t T−t t T with a fix maturity T, is in fact a function of t and S . In Section 3, for notational simplicity, we suppress the t dependence of f on t in the definition C =n E[ΠT|S ]≡f(S ). t t nT t t 9

techniquesofGeman-ElKaroui-Rochet[10]foroptiontypecontractswithatmostthreedistinct sources of randomness: stochastic short rate and a maximum of two risky assets. Well known examples of these contracts are options written on stocks or bonds, e.g. European options and exchange options. Proposition 1. Consider the credit exposure process, {V ;t ≥ 0}, defined on a given filtered t probability space (Ω,F,(F ) ,P), and a T-maturity transaction between the financial instit 0≤t≤∞ tution and its counterparty that is unilateral, i.e. the credit exposure process is the price process, V = C > 0 for all 0 ≤ t ≤ T, where C denotes the time-t value of the derivative contract with t t t payoff Π . Then, T cov(V ,V ) > 0, u t for any 0 < u < t < T under any of the following conditions: Condition 1: Numeraire is the money market account, B, with deterministic short rate, r, and Π is a function of N ≥ 1 exogenously given risky assets. T Condition 2: Short rate is stochastic and the T-payoff function is a function of at most two risky assets as follows Π = (α S (T)+α S (T))+, where α and α are any real numbers, and T 1 1 2 2 1 2 S and/or S are risky assets. 1 2 In the case of bilateral transactions, e.g. interest rate swaps, for which the exposure process satisfies V = max{C ,0} (cid:54)= C for some 0 < t ≤ T, where C denotes the time-t value of the t t t t derivative contract with payoff function Π , stronger assumptions are required to analytically T show that cov(V ,V ) > 0 for any 0 < u < t < T. This is shown in Proposition 2 below. u t Proposition 2. Consider the credit exposure process, {V ;t ≥ 0}, defined on a given filtered t probability space (Ω,F,(F ) ,P), and a T-maturity transaction between the financial instit 0≤t≤∞ tution and its counterparty that is bilateral, i.e. the credit exposure process is the price process, V = max{C ,0} =(cid:54) C for some 0 < t ≤ T, where C denotes the time-t value of the derivative t t t t contract with payoff function Π . Then, T cov(V ,V ) > 0 u t for any 0 < u < t < T under the following condition: Numeraire is the money market account, B, with deterministic short rate, r, and Π is a T monotone function of a single risky asset whose dynamics is modeled by a GBM, a Gauss- Markov process or a square-root diffusion. Themonotonicityassumptionofthepayofffunctionissatisfiedformostoftheactivelytraded OTC derivative contracts; well-known exceptions are Barrier6 and Lookback options, (see, for instance, [22]). 6Morespecifically,thepayofffunctionofup-and-inanddown-and-outEuropeanbarriercalloptionsaremonotonefunctionsoftheunderlyingsecurityprices. Thismonotonicityassumptiondoesnotholdforup-and-outand down-and-in European barrier call options, (see Chapter 6 of [22] and the references there). 10

Propositions1 and 2 identifyconditions for unilateraland bilateral transactions under which the credit exposure process satisfies cov(V ,V ) > 0 for any 0 < u < t < T. This, then, implies u t that Var(θˆ ) ≤ Var(θˆ ). (10) b,m,n,d b,m,n,p Note that the above inequality holds since Var(θˆ ) = (cid:88) n Var(V i )∆2 i ≤ (cid:88) n Var(V i )∆2 i + 2 (cid:88)(cid:88) cov(V ,V )∆ ∆ = Var(θˆ ).(11) b,m,n,d i j i j b,m,n,p m m m i=1 i=1 i<j 3.1.2 Path Dependent Case Suppose that V is time t value of a maturity-T contract, where the payoff at the time T is t a function of S ,...S , (for instance, an arithmetic Asian option). That is, V = g(S ,...,S ), 1 n i 1 i where g is a function from Ri to R. The DJS sampling method is to make V = g(S ,...,S ) i 1 i and V = g(S ,...,S ), i < j, uncorrelated random variables. That is, sample from S ,...,S j 1 j 1 i to generate a single realization of V . To generate V , start again from time zero, and sample i j from S ,...,S ,...S . Under this DJS-type sampling method, V and V become uncorrelated, 1 i j i j cov(V ,V ) = 0. In the PDS-type sampling, given the Monte Carlo realization of V , to generate i j i V , one uses the previously generated S ,...,S and only samples from S ,...,S . In this case j 1 i i+1 j V and V are dependent. i j Using conditional covariance formula and arguments similar to the ones used in the path independent case, it can be shown that cov(V ,V ) > 0, i (cid:54)= j. More specifically, it can be shown i j that cov(V ,V ) > 0 for unilateral and bilateral transactions under the first condition of Propoi j sition 1 and Proposition 2’s condition, respectively. That is, for the above mentioned covariance function to be positive, we need the numeraire money market account with deterministic short rate in the unilateral case. The bilateral case, additionally, requires monotonicity of the payoff function and its dependence on a single risk factor. To compare the efficiency of the DJS and PDS-based estimators of θ in the path dependent case, computing time is also to be considered in parallel with variance of the estimators.7 More specifically, the estimator with the lower variance per replication×expected computing time, should be selected (see [15] for the formal formulation of this useful criterion in comparing alternative Monte Carlo estimators). Consider, for instance, arithmetic Asian options. Suppose 7In the path independent case computing time of DJS and PDS-based estimators of θ are roughly equal. 11

that the computational time to calculate θˆ is proportional to the number of random varib,m,n,k ables that are to be generated. Let ct(θˆ ) denote the computational effort associated with b,m,n,k θˆ . Note that, b,m,n,k ct(θˆ ) Var(θˆ ) b,1,n,d b,1,n,p ≈ n and ≈ n. (12) ct(θˆ ) Var(θˆ ) b,1,n,p b,1,n,d To see why (12) holds note that to calculate θˆ , n(n+1) random variables are to be generated b,1,n,d 2 while θˆ requires generating n random variables, (assuming that the calculation of E[ΠA|F ] b,1,n,p i does not require generating additional random variables). Also, note that as can be seen from (11), variance of the PDS-based estimator is of order n2 because of the covariance terms while the DJS-based estimator has a variance of order n . So, θˆ and θˆ have a similar b,m,n,d b,m,n,p performanceforfixedandsufficientlylargen. PDSandDJS-basedestimatorsofotherderivatives whose payoff depends on the path in a different form can be compared similarly. We emphasize that what we have referred to as the “DJS method” for the path dependent case can be considered as a generalized version of what the conventional CCR literature often refers to as “DJS”. While clarifying this point, Appendix G provides a further explanatory discussionoftheDJSmethod. AppendixH, whichdiscussessettingsthathavenotbeenformulated in our study, illustrates the potential applicability of the DJS method for these more general cases. 3.1.3 Summary of Section 3.1 We summarize the result of Section 3.1 as it is used in the sequel and is directly applied to the efficient CVA estimation.8 To compare the DJS and PDS-based estimators of EPE and CVA I I (both being viewed as weighted sums of expected exposures) variance and computing time of the Monte Carlo estimators are considered. The DJS method induces zero covariance between any two distinct time points of the simulated credit exposure process. So, it remains to look at this covariance function for the credit exposure process under the PDS method. When the dynamics of the risk factors are modeled by the class of continuous time stochastic processes considered in this paper, the covariance function of the credit exposure process under the PDS method becomes positive under conditions of Propositions 1 and 2 for unilateral and bilateral path independent derivatives transactions, respectively. Similar results hold for path dependent derivatives. That is, under conditions of Propositions 1 and 2, DJS-based estimators of EPE and CVA outperform the PDS-based estimators in terms of variance. For path independent derivatives PDS and DJS-based computing times are roughly equal. So, we recommend that 8Section9.3.2ofGregory[18]givesanoverviewofthePDS,referredtoaspathwisemethod,andDJS,referred to as direct simulation, and broadly concludes that the pathwise method is more suitable for path-dependent derivatives. Section 12.4.4 of [18] discusses the PDS and DJS methods in Monte Carlo estimation of CVA where the counterparty’s default is modeled by a structural-copula approach. In this context, [18] gives numerical results,Figure12.10and12.11of[18],illustratingsuperiorstatisticalefficiencyofDJSoverPDSforMonteCarlo estimation of CVA when the counterparty’s default is modeled by a structural-copula method. 12

the counterparty credit risk modeler uses DJS for path independent derivatives. DJS-based estimators usually have larger computing times for path dependent derivatives, and so the computing time should be considered in parallel with the variance. The useful criterion that does this is ‘variance per replication×expected computing time’ as introduced in the previous section. There are widely traded path dependent derivatives for which PDS and DJS-based estimators of EPE and CVA perform approximately equally. For instance, for arithmetic Asian options the DJS and PDS-based estimators of EPE and CVA perform similarly. There are contracts whose payoff function does not exactly match the mathematical conditions of Proposition 1 and 2. For those contracts, a small simulation study could compare the variance of the DJS and PDS-based estimators of EPE and CVA. The Appendix contains numerical examples for EPE estimation of a single interest rate swap, where we conclude that DJS outperforms PDS by at least an order of 10 in terms of variance while the computing times of both Monte Carlo estimators are roughly equal. Hereafter, we assume that the credit exposure process V satisfies cov(V ,V ) = 0 and u s cov(V ,V ) > 0whensimulatedundertheDJSandPDSmethods,respectively,forany0 < u < t. u s 3.2 Efficient Monte Carlo EPE Estimation: Biased Estimators In this subsection, we suppress the subscript b in θˆ and instead write θˆ for notational b,m,n,k m,n,k simplicity. Wewouldliketofindthenumberofvaluationpoints,n,andthenumberofsimulation runs at each valuation point, m, to minimize MSE(θˆ ), m,n,k MSE(θˆ ) = Var(θˆ )+(E[θˆ ]−EPE)2. m,n,k m,n,k m,n,k given a fixed computational budget, denoted by s, that is proportional to, mn. Also, k = p and d refer to PDS and DJS-based simulation of the credit exposure process on a time grid 0 ≡ t < t < ... < t ≡ T. That is, as shown in the previous section, under PDS sampling and 0 1 n DJS sampling, cov(V ,V ) > 0 and cov(V ,V ) = 0, respectively, for any i (cid:54)= j, i,j = 1,...,n. i j i j To formulate and solve this optimization problem, we specify the order of the variance and bias of the Monte Carlo estimator of EPE, θˆ . Note that from basic results on endpoint m,n,k Reimann sum approximation of integrals, time-discretization bias is of order 1/n. We are not concerned with deriving sharp estimates of the orders of variance. In fact, our numerical examples indicate that choosing approximately optimal m and n using even very rough approximates for the orders of variance and bias leads to substantial MSE reduction compared to industry practice. Suppose that the time grid is equidistant, i.e., ∆ ≡ ∆ = T. We assume that E[V2] < ∞ i n t for all t ∈ [0,T]. First, we note that 1 Var(θˆ ) = O( ). (13) m,n,d mn To see this,9 consider M > 0 such that E[V2] ≤ M for t ∈ (0,T]. Note that, t 9The Landau symbol, O, in f(x,y)=O(g(x,y)) means that f(x,y)/g(x,y) stays bounded in some limit, say x,y→0 or x,y→∞. 13

Var(θˆ ) = ∆2 (cid:88) n Var(V i ) ≤ ( T )2 (cid:88) n E(V i 2) ≤ MT2 . m,n,d m n m mn i=1 i=1 Now, consider the variance of the PDS-based estimator, θˆ , m,n,p n Var(θˆ ) = ∆2 (cid:88) Var(V i ) +∆2 2 (cid:88)(cid:88) cov(V ,V ). m,n,p i j m m i=1 i<j As shown before, the first term above is O( 1 ). Also, under PDS sampling, the credit exposure mn process is simulated according to its finite dimensional distributions for which the covariance terms are positive. So, the second term is O(1). This gives, m 1 1 Var(θˆ ) = O( + ). (14) m,n,p mn m PDS-Based Biased Efficient Estimator of EPE We choose the number of valuation points, n, and number of simulation runs at each valuation point, m, to minimize the mean square error of the PDS-based estimator, θˆ , under a fixed computational budget proporm,n,p tionaltomn. Approximatingthevarianceofθˆ using(14)leadstothefollowingoptimization m,n,p problems, (cid:16)c c c (cid:17) p,1 p,2 2 min + + subject to s = c mn, (15) m,n mn m n2 3 for some constants, c ,c ,c , and c . MSE of θˆ is minimized at, p,1 p,2 2 3 m,n,p 2 1 m = cs3 and n = c˜s3, (16) for constants c and c˜. DJS-Based Biased Efficient Estimator of EPE Let c denote a constant. Given (13), we d approximateVar(θˆ )with c d intheMSEminimizationproblemfortheDJS-basedestimator, m,n,d mn (cid:16) c c (cid:17) d 2 min + subject to s = c mn, m,n mn n2 3 to which the trivial optimal solution is m = 1 and n = cˆs fo some constant cˆ. We note that estimating the various constant parameters appearing in all the above mentioned MSE minimization problems is not possible in practice. In our numerical examples we simply set all these constant parameters equal to 1. 14

Remark We do not claim originality in setting up an MSE minimization problem to derive an optimumbalancebetweenvarianceandbiassquared;thiscanbeseeninChapter6ofGlasserman [13] and the references there, particularly the paper by Duffie and Glynn [7]. Our contribution is that in our proposed efficient Monte Carlo CCR framework, choosing approximately optimal m and n via solving MSE minimization problems achieve substantial MSE reduction in Monte Carlo estimation of EPE, CVA, (and as will be seen in Section 5), and EEPE. This has neither appeared in the CCR literature nor been applied by practitioners. Moreover, our result that the efficient DJS-based estimator requires all its computational budget allocated to the number of valuation points is surprising. Consider a well defined continuous time stochastic process S whose finite dimensional distributions satisfy cov(S ,S ) > 0 for any 0 < u < t. Suppose that u t (cid:82)T Monte Carlo is to be used to estimate θ = E[ S dt] for a given T > 0. The porposed efficient 0 t Monte Carlo estimator of θ employs the notion of marginal matching (as opposed to simulating the process based on the finite dimensional distribution of S) and uses 1 simulation run at each time point in a discrete time grid given a fixed computational budget all allocated to making the grid as fine as possible. 3.3 Efficient Monte Carlo EPE Estimation: Unbiased Estimators In this section we derive unbiased estimators of EPE. Specifically, we eliminate the time discritization bias at the expense of introducing additional randomness. To control the variance that would be increased as the result of this new source of randomness, we use stratified sampling. Let τ denote a [0,T] Uniform random variable that is independent of the credit exposure, V. We have, EPE = TE[V ], (17) τ which simply follows from conditioning on τ, i.e., using E[V ] = E[E[V |τ]], independence of τ τ V and τ, and noting that f(t) = 1, t ∈ [0,T], is the probability density function of τ. Now, T consider the following identity, n n (cid:88) (cid:88) EPE = TE[V ] = T E[V |τ ∈ A ]p = E[V |τ ∈ A ]∆ , (18) τ τ i i τ i i i=1 i=1 where A = [0,t ), p ≡ P(τ ∈ A ) = ∆i, on the time grid, 0 ≡ t < t < ... < t ≡ T, and i i i i T 0 1 n ∆ = t −t . Assuming t = iT/n for all i = 1,...,n, our proposed unbiased estimators of EPE i i i−1 i use the identity (18) by estimating the conditional expectations, E[V |τ ∈ A ], i.e., τ i n (cid:88) θˆ = V¯ ∆ , (19) u,m,n,k τi i i=1 where τ ≡ τ|τ ∈ A , V¯ = (cid:80)m V /m, and τ ,...,τ are i.i.d. copies of τ . That is, to i i τi j=1 τij i1 im i draw a single realization of V , we first sample from τ conditional on τ ∈ A . Note that τi i τ is a [t ,t ] Uniform random variable. Next, given this realization of τ , we generate V . i i−1 i i τi 15

The subscript k = p and d refer to PDS and DJS sampling, respectively.10That is, PDS-based simulation in calculating θˆ implies that cov(V ,V ) > 0 for i (cid:54)= j, i,j = 1,...,n, and u,m,n,p τi τj DJS-based simulation in calculating θˆ implies that cov(V ,V ) = 0 for i (cid:54)= j. This u,m,n,d τi τj immediately implies Var(θˆ ) ≤ Var(θˆ ). Consider a more general setting that allows u,m,n,d u,m,n,p different numbers of simulation runs for each stratum. That is, let m denote the number of i runs used to estimate E[V |τ ∈ A ] and N = m +...+m denote the total number simulation τ i 1 n runs. Note that our setting with equidistant strata and m ≡ m, for i = 1,...,n coincides i with proportional stratified sampling which uses m = Np , (see [26] for results on proportional i i stratification). This is because τ is a [0,T] Uniform random variable. In this paper we do not address further possible improvements of our unbiased stratified sampling-based estimators of EPE by attempting to find optimal m ,...,m and n under fixed computational budgets. Our 1 n numerical examples indicate that using our unbiased stratified sampling-based estimators by setting m ≡ m and choosing m and n as specified in subsection 3.2 leads to substantial MSE i reduction when compared to crude biased Monte Carlo estimators of EPE. Comparing DJS-based Biased and Unbiased estimators Proposition 3 below shows that θˆ and the biased DJS-based estimator of EPE, θˆ , are asymptotically equivalent u,m,n,d b,m,n,d in terms of MSE. This equivalence is further confirmed by our numerical experiments (see the next subsection) in practical settings with fixed and finite computational budgets proportional to mn. Proposition 3. Consider the credit exposure process, {V ;t ≥ 0}, defined on a given filtered t probabilityspace(Ω,F,(F ) ,P). SupposethatbiasedandunbiasedMonteCarloestimators t 0≤t≤∞ of EPE calculated under DJS-sampling, n n (cid:88) (cid:88) θˆ = V¯∆ , and θˆ = V¯ ∆ . (20) b,m,n,d i i u,m,n,d τi i i=1 i=1 are defined on an equi-distant time grid, 0 ≡ t < t < ... < t ≡ T, where ∆ ≡ t −t = 0 1 n i i i−1 T/n ≡ ∆, τ ≡ τ|τ ∈ A and A = [t ,t ). Let V¯ and V¯ denote the averages of m Monte i i i i−1 i i τi Carlo realizations of V , and V , respectively. That is, the total number of simulation runs is i τi N = mn. We assume that E[V2] < ∞, for all i = 1,...,n. Asymptotic performance of θˆ i b,m,n,d and θˆ is equivalent in the following sense, u,m,n,d (cid:90) T lim nMSE(θˆ ) = nVar(θˆ ) = c Var(V )dt, (21) b,m,n,d u,m,n,d t n→∞ 0 where c is a constant. 10RecallthatthebiasedestimatorsofEPE,θˆ ,k=p,d,arebasedonRightReimansumapproximationof m,n,k the integral of the expected exposures in the EPE formula. Our proposed unbiased estimators θˆ , k=p,d, u,m,n,k can simply be viewed as a Reiman sum approximation of the EPE where each expected exposure is evaluated at a randomly selected point within each subinterval. 16

Comparing PDS-based Biased and Unbiased estimators Analytically comparing the MSE(θˆ ) and Var(θˆ ) is quite difficult due to the presence of the covariance terms. b,m,n,p u,m,n,p Our numerical examples presented in the next subsection show that the unbiased PDS-based estimator of EPE, θˆ , outperforms the efficient biased PDS-estimator, θˆ , introduced u,m,n,p b,m,n,p in the previous section. 4 Efficient Monte Carlo Estimation of Independent CVA Independent CVA can be viewed as the weighted sum of expected exposures with the weights being default probabilities. Therefore, our results from Section 3 on efficient estimation of EPE immediately apply here (note that for EPE, the weights are subinterval lengths). To summarize our results on efficient Monte Carlo CVA estimation, we suppress the dependence of CVA on I the stochastic discount factor by assuming zero short rate, (cid:90) T CVA = E[E[V 1{τ ≤ T}|τ]] = E[V ]dF , (22) I τ t t 0 whereF denotesthecumulativedistributionfunctionofτ,whichisassumedtobeknown(market observable) from, for instance, credit default swap spreads of the counterparty, (see, e.g., [19]). The Remark at the end of this section discusses Monte Carlo CVA estimation with stochastic discounting. Appendix I discusses PDS and DJS-based estimation of CVA sensitivities with respect to the initial values of underlying risk factors using finite-difference approximations. Efficient Biased Estimators of CVA We can employ our MSE minimization formulation I to first specify the approximately optimal n and m under a fixed computational budget, and then estimate CVA with I n (cid:88) ξ = V¯∆F , (23) b,k i i i=1 where k = p,d denotes PDS and DJS sampling, respectively, V¯ = (cid:80)m V /m as defined in i j=1 ij Section 3, and ∆F ≡ F(t )−F(t ). (We have suppressed the dependence of ξ on m and n, i i i−1 b,k i.e., ξ ≡ ξ .) b,k b,m,n,k Efficient Unbiased Estimators of CVA Note that I n (cid:88) E[V 1{τ ≤ T}] = E[V |τ ∈ A ]P(τ ∈ A ), (24) τ τ i i i=1 where stratum i is A = [t ,t ). Let m , i = 1,2,...,n denote the number of simulation runs i i−1 i i used to estimate E[V ], where V ≡ V , t ≡ 0, and t = T. Also, N = (cid:80)n m denotes the i i ti 0 n i=1 i total number of simulation runs used in estimating CVA . Using τ as the stratification variable I and the identity (24), the stratified sampling estimator of CVA is I 17

n (cid:88) ξ = V¯ p , (25) u,k τi i i=1 where k = p,d denotes PDS and DJS sampling, respectively. Also, p ≡ P(τ ∈ A ) = ∆F , i i i τ ≡ τ|τ ∈ A , and V¯ = (cid:80)mi V /m . That is, to draw a single realization of V , we first i i τi j=1 τij i τi sample from τ conditional on τ ∈ A ; next, given this realization of τ , we generate V . In terms i i τi of computing time, ξ requires generating N realizations of V and ξ requires N additional b,k i u,k samples from the truncated τ based on the strata defined above. Note that since generating V is computationally much more intensive than the truncated τ, ξ outperforms ξ merely i b,k u,k marginally in terms of the computational time. As mentioned before, proportional stratified sampling sets m = NP(τ ∈ A ). For CVA i i I estimation, even if we assume A ’s to be equidistant strata, P(τ ∈ A )’s are not equal in general. i i Therefore, proportional stratified sampling does not lead to an equal number of simulation runs at all the valuation points, as is the case for EPE estimation (see Section 3.3). We have empirically observed that ξ outperforms ξ in terms of mean square error.11 u,p b,p Proposition4belowshowsanasymptoticequivalencebetweenξ andξ . TheproofofPropob,d u,d sition 4 is similar to Proposition 3, and so it is omitted. Also, similar to our numerical examples to be shown in Section 6.1, we have observed that DJS-based biased and unbiased estimators of CVA are equivalent in terms of MSE for large n, which confirms the result of Proposition 4 I below. Proposition 4. Consider the proposed estimators of CVA , ξ and ξ as defined in (23) I b,d u,d and (25), respectively. Suppose that proportional sampling is used for both estimators, i.e., m = Np , and (cid:80)n m = N, i = 1,...,n. We assume that E[V2] < ∞, i = 1,...,n. Note that i i i=1 i i DJS sampling gives cov(V ,V ) = 0 for all i (cid:54)= j and i,j = 1,...,n. Then the following holds, i j (cid:90) T lim nVar(ξ ) = nMSE(ξ ) = c Var(V )dF(t) (26) u,d b,d t n→∞ 0 where c is a constant and F is the cumulative distribution function of τ. That is, ξ and ξ b,d u,d perform similarly in terms of asymptotic MSE. SimilartoourproposedefficientEPEestimationframework, whentheDJSmethodischosen for sampling from the credit exposure process, we recommend using efficient biased CVA estimationwheretheapproximatelyoptimalmandnarederivedviasolvingtheMSEminimization problem discussed in Section 3.2. Note that the surprising approximately optimal solution sets m = 1 and allocates the total computational budget to the number of valuation points. When the PDS method simulates the credit exposure process, we suggest our unbiased stratified sampling based estimator of CVA , where the number of strata and simulation runs at each stratum I are specified based on the solution to the MSE minimization problem as in Section 3.2. 11Assuming that τ is an exponential random variable, the results of our numerical examples for ξ and ξ u,p b,p are similar to those to be shown in Section 6.1, and so are omitted. 18

Remark This remark clarifies and discusses the potential application of the DJS method in Monte Carlo estimation of (cid:90) T CVA = E[D V ]dF , I t t t 0 in the presence of stochastic discounting.12 Let D t = e− (cid:82) 0 trudu, 0 < t ≤ T, and assume the risk neutral dynamics of the short rate is specified by a widely used short rate model, (see, e.g., Chapter 21 of [2]). Also, suppose that the valuation points of the Monte Carlo CCR “engine” are specified by the equidistant time grid 0 ≡ t < t < ... < t ≡ T; t −t = ∆ and that the 0 1 n i i−1 time integral of the short rate in D is time-discretized on a finer grid 0 ≡ t˜ < t˜ < ... < t˜ ≡ T; 0 1 n˜ t˜ − t˜ = ∆˜ < ∆. Assume that ∆ = c∆˜; c being a positive integer. The DJS-based estij j−1 mator of CVA seeks to make D V ’s independent at any two distinct time points on the gird. I i i Suppose that the V ≡ V is a function of S ,...,S , which are dependent random variables, i ti 1 i i = 1,...,n. Note that {r ≡ r ;j = 1,..,n˜} represents the finite dimensional distribution of r i t˜ j on the finer grid. In the DJS-based single-simulation-run estimator of CVA , we are to have I cov(D V ,D V ) = 0 for any i (cid:54)= j, i,j = 1,..,n. Achieving this zero correlation requires adi i j j ditional computational time due to the generation of additional random variables. To see this, consider two distinct valuation points on the grid, t < t . A single Monte Carlo realization of k l D V requires generating (r ,r ,...,r ) and (S ,...,S ). To sample from D V , the DJS method, k k 1 2 ck 1 k l l instead of using the Monte Carlo realizations of the risk factor in the simulation of D V , starts k k anew and generates (r ,r ,...,r ,...,r ) and (S ,...,S ,...,S ), and so, cov(D V ,D V ) = 0. 1 2 ck cl 1 k l k k l l One should contrast this with the PDS method. In sampling from D V , the PDS method generl l ates only (r ,...,r ) and (S ,...,S ) along with the random variables generated previously ck+1 cl k+1 l in the simulation of D V . So, the DJS method is applicable in the context of Monte Carlo CVA k k estimation with stochastic discounting at an additional computational cost. Suppose the CCR modeler analytically, using Proposition 1 and 2 or their equivalents, or numerically, via a small simulation study, concludes that the PDS-based estimator has a higher variance, due to some positive covariance terms, cov(D V ,D V ) > 0 for i (cid:54)= j, i,j = 1,..,n, compared to that of the i i j j DJS-based estimator. Then, the DJS-based method is chosen only when the computing time has also been taken in to account. That is, the modeler chooses the DJS-based estimator of CVA only if it outperforms the PDS-based estimator in term of the following useful criterion, I variance per replication×expected computing time, as we discussed in Section 3.1.2. 5 Efficient Monte Carlo Estimation of EEPE InthissectionwediscussefficientMonteCarloestimationofeffectiveexpectedpositiveexposure, EEPE, (cid:90) T EEPE = max E[V ]dt, u 0 0≤u≤t 12We continue to assume that the recovery rate is constant, and to simplify the notation, we suppress the dependence of CVA to it. 19

where {V ;t ≥ 0} denotes the credit exposure process, and T denotes the expiration time of t the transaction with the longest maturity in a portfolio of OTC derivatives held by a financial institution with its counterparty. Consider the time grid, 0 ≡ t < t < ... < t ≡ T. Set 0 1 n ∆ ≡ t −t , i = 1,...,n. Monte Carlo estimators of EEPE are, i i i−1 n (cid:88) θˆ = max{V¯ }∆ , (27) m,n,k j i 1≤j≤i i=1 where V¯ denotes the m-simulation run average of the i.i.d. random variables, V ,...,V . The j j1 jm subscriptk = pandddenotePDSandDJSsampling,respectively. Thatis,underk = p(k = d), V ’s are positively correlated (uncorrelated). Consider the mean square error of θˆ , j m,n,k (cid:32) n (cid:33) (cid:32) n (cid:33)2 (cid:88) (cid:88) MSE(θˆ ) = Var max{V¯ }∆ + E[max{V¯ }]∆ −EEPE . (28) m,n,k j i j i 1≤j≤i 1≤j≤i i=1 i=1 Bias Decomposition It is useful to differentiate the following two sources of bias, (cid:32) n n (cid:33) (cid:32) n (cid:33) (cid:88) (cid:88) (cid:88) E[max{V¯ }]∆ − max E[V¯ ]∆ − EEPE− max E[V¯ ]∆ . (29) j i j i j i 1≤j≤i 1≤j≤i 1≤j≤i i=1 i=1 i=1 That is, the first part of the bias is due to the presence of the maximum operator and the second part is time-discretization bias. Note that for a fixed n, variance of θˆ converges to zero as m,n,k m → ∞. Now, consider Proposition 5 below whose proof is in the Appendix. Proposition 5. Let {V ; t ≥ 0} denote the credit exposure process. Let, t M ≡ max{V¯ ,...,V¯ }, n,m,k 1 n where V ≡ V on the time grid 0 ≡ t < t < ... < t ≡ T, and V¯ = (cid:80)m V /m, V ,...,V are i ti 0 1 n i j=1 ij i1 im i.i.d random variables. Also, k = d and k = p refer to the cases where V are uncorrelated and i positively correlated, respectively, resulting from DJS and PDS-based simulation of V. Assume that E[V2] < ∞ for all i = 1,...,n. Let M ≡ max{E[V ],...,E[V ]}. Then, as m → ∞, i n 1 n M → M a.s., (30) n,m,k n where a.s. stands for almost surely. Note that dominated convergence theorem and Proposition 5 give E[M ] → M as n,m,k n m → ∞.13 So, the first part of the bias n n (cid:88) (cid:88) E[max{V¯ }]∆ − max E[V¯ ]∆ j i j i 1≤j≤i 1≤j≤i i=1 i=1 13Note that M ≤ (cid:80)n V¯ and Proposition 5 assumes E[V¯]=E[V ]<∞. n,m,k i=1 i i i 20

converges to zero as m → ∞. That is, θˆ and θˆ are consistent estimators of EEPE m,n,d m,n,p dst for a fixed n. In what follows we first show that for a fixed n and sufficiently large m, θˆ outperforms m,n,d θˆ in terms of variance. Next, after specifying approximates for the order of variance and m,n,p bias of θˆ , we formulate an MSE minimization problem over m and n given a fixed compum,n,k tational budget. Our numerical results indicate that our proposed estimators of EEPE, which use approximately optimal m and n, lead to substantial MSE reduction when compared to the crude estimators. 5.1 Comparing PDS and DJS-based Monte Carlo Estimators of EEPE We are to compare the variance of θˆ and θˆ for a fixed n and sufficiently large m. Set m,n,p m,n,d n n (cid:88) (cid:88) θ ≡ max E[V ]∆ , θˆ ≡ max{V¯ }∆ , j i m,n,k j i 1≤j≤i 1≤j≤i i=1 i=1 where k = p (k = d) refer to the cases where V and V for any i (cid:54)= j, are positively correlated i j (uncorrelated). In what follows we find it useful to append a second subscript m to V¯ to i emphasize that the average is based on m i.i.d random variables and a third subscript k = d or p to indicate DJS or PDS. Denote by δ ≡ E[V ]−E[V ] and δ ≡ min{|δ | : i (cid:54)= j,i,j = 1,...,n}. Without loss of i,j i j i,j generality, assume δ > 0. Let σ denote the standard deviation of V − V under estimai,j,k i j tion method type k and σ ≡ max{σ : i,j = 1,...,n,k = d,p}. For i = 1,...,n, let τ max i,j,k i denote the index for which max{E[V ],...,E[V ]} is attained and τ be the index for which 1 i i,m,k max{V¯ ,...,V¯ } is achieved. It then follows from these definitions that 1,m,k i,m,k n n (cid:88) (cid:88) θ = E[V ]∆ and θˆ = V¯ ∆ . (31) τi i m,n,k τ i,m,k ,m,k i i=1 i=1 For k = d or p and i = 2,...,n, the probability that simulations do not yield the right τ can be i 21

bounded from above as follows (cid:88) P(τ (cid:54)= τ ) ≤ P(V¯ −V¯ < 0) i,m,k i τi,m,k j,m,k j(cid:54)=τi,j=1,...,i (cid:88) = P(V¯ −V¯ −δ < −δ ) τi,m,k j,m,k τi,j τi,j j(cid:54)=τi,j=1,...,i (cid:88) < P(V¯ −V¯ −δ < −δ) τi,m,k j,m,k τi,j j(cid:54)=τi,j=1,...,i (cid:88) < P(|V¯ −V¯ −δ | > δ) τi,m,k j,m,k τi,j j(cid:54)=τi,j=1,...,i σ2 ≤ (cid:88) τi,j,k (32) mδ2 j(cid:54)=τi,j=1,...,i (i−1)·σ2 ≤ max, mδ2 where (32) follows from the Chebyshev’s inequality. Consider the event B = {τ = τ = τ ,for all i = 1,...,n}. It makes sense to call B m i i,m,d i,m,p m the desirable event and Bc the undesirable event. Let θˆ denote θˆ conditional on m m,n,k,Bm m,n,k the event B . We have that m Var(θˆ ) < Var(θˆ ). (33) m,n,d,Bm m,n,p,Bm This order can be established by first noting that n (cid:88) θˆ = V¯ ∆ . (34) m,n,k,Bm τi,m,k i i=1 Then since V and V , for any i (cid:54)= j, are positively correlated (uncorrelated) under PDS (DJS) i j sampling, the variance of expression (34) is lower under DJS than under PDS. Note that: n n (cid:88) (cid:88) P(Bc ) ≤ P(τ (cid:54)= τ )+ P(τ (cid:54)= τ ) (35) m i,m,d i i,m,p i i=2 i=2 (cid:88) n (i−1)·σ2 < 2 max, (36) mδ2 i=2 The above argument leads to the following result. Proposition 6. Consider the desirable event B as defined above. First, conditional on this m event, (33) holds. Secondly, the desirable event occurs asymptotically almost surely as m → ∞. That is, lim P(B ) = 1. More specifically, P(Bc ) goes to zero at rate 1/m as m → ∞. m→∞ m m Proposition6suggeststhatforsufficientlylargem,Var(θˆ ) ≤ Var(θˆ ). Ournumerical m,n,d m,n,p examples in Section 6.2 use m > 400; they all indicate that Var(θˆ ) ≤ Var(θˆ ). m,n,d m,n,p 22

5.2 Efficient Monte Carlo Estimation of EEPE Similar to our approach in subsection 3.2, we would like to find the number of valuation points, n, and the number of simulation runs at each valuation point, m, to minimize MSE(θˆ ) m,n,k given a fixed computational budget, s, that is proportional to, mn. To do so, we need to specify the order of the variance and bias of the Monte Carlo estimator of EEPE, θˆ . We m,n,k are not concerned with deriving sharp estimates of the orders of variance and bias. In fact, our numerical examples indicate that choosing approximately optimal m and n using even very rough approximates for the orders of variance and bias lead to substantial MSE reduction. The following is used to formulate our MSE minimization problem: for k = p or d, c c c c Var(θˆ ) ≈ 1,k + 2,k and Bias(θˆ ) ≈ 3,k + 3 , (37) m,n,k m,n,k mn m m n for some constants c ,c ,c . The above approximation of the order of bias uses (29) and 1,k 2,k 3 Proposition 5. Note that our rough approximate of the order of variance, applicable to both θˆ and θˆ , is similar to (14). This is because of the presence of the maximum operators m,n,d m,n,p that leads to positive covariance terms. To see this, let V¯ denote the m-simulation-run average i of the i.i.d random variables, V ,...,V , i = 1,...,n, and consider an equidistant time grid with i1 im n time points, ∆ = T/n. Note that, Var(θˆ ) = ∆2Var (cid:0) V¯ +max{V¯ ,V¯ }+...+max{V¯ ,...,V¯ } (cid:1) , m,n,k 1 1 2 1 n isequalto∆2 = T2 timesthesumofnnon-zerovariancetermsandn(n−1)/2positivecovariance n2 terms both for k = d and k = p. This leads to a result similar to (14). Given (37), we recommend solving the following MSE minimization problem to specify the approximately optimal m and n, (cid:16) c c c c (cid:17) min 1 + 2 +( 3 + 4 )2 subject to s = cmn, (38) m,n mn m m n for some constants c ,c ,c ,c , and c. 1 2 3 4 6 Numerical Examples 6.1 EPE estimation In this section we use simple numerical examples to illustrate the efficiency of our proposed Monte Carlo estimators of EPE. We consider contract level exposure in a simple setting where V t ≡ S t denotes the value of a geometric Brownian motion at time t > 0. That is, S t = S 0 eXt with {X ; t ≥ 0} being a Brownian motion with drift µ, and volatility σ. This stylized example t enables us to calculate the MSE exactly. We consider six different Monte Carlo estimators of EPE in our numerical examples. Let θˆ and θˆ denote the “crude” and biased Monte Carlo estimators of EPE under PDS c,p c,d and DJS sampling, respectively. That is, 23

n (cid:88) θˆ = V¯∆ , (39) c,k i i i=1 where ∆ = t − t , 0 ≡ t < t < ... < t ≡ T, k = p,d, and V¯ is the m-simulation-run i i i−1 0 1 n i average of V . We shall shortly specify the choice of the valuation points. i Letθˆ andθˆ denotetheefficientandbiasedMonteCarloestimatorsofEPEunderPDS e,b,p e,b,d and DJS sampling, respectively. In particular, their statistical efficiency is a result of solving the MSE minimization problems in Section 3.2 to derive the (approximately) optimal number of points on the time grid, n, and simulation runs at each of these time points, m, given a fixed computational budget proportional to mn. Let θˆ and θˆ denote the unbiased stratified sampling-based Monte Carlo estimators of u,p u,d EPE under PDS and DJS sampling, respectively. That is, n (cid:88) θˆ = V¯ ∆ , (40) u,k τi i i=1 where V¯ = (cid:80)mi V /m with τ ≡ τ|τ ∈ A , A = [t ,t ], and k = p,d. τi j=1 τij i i i i i−1 i We set T = 1. The crude estimators of EPE are calculated based on 12 valuation points, n = 12, at 1, 2, 3, 4, 8, 12, 18, 21, 24, 36, 49 weeks and 1 year. We note that one year, T = 1, withthenumberofvaluationpointsfixedat12, isasettingwidelyusedbyfinancialinstitutions. Thereisnomathematicalbasisforthisarrangementofvaluationpoints. Itisbelievedthatsince some trades have “short” expiration times, having more valuation points earlier would increase the accuracy of the estimators of CCR measures. The time grid used to calculate our efficient estimators of EPE is equidistant, i.e., ∆ ≡ ∆ = T/n . Computational budget, s, is fixed at i 12,000 and 120,000, respectively. Therefore, the number of simulation runs at each valuation point m is 1,000 and 10,000, respectively, for the crude estimator subject to these two budget values. To calculate θˆ under these fixed computational budgets, the solution, (16) with both c e,b,p and c˜set to 1, to the MSE minimization problem of Section 3.2 is used. This gives, n = 23 and m = 524 for s = 12,000, and n = 50, and m = 2433 for s = 120,000. Similarly, to calculate θˆ , we use the solution to the MSE minimization problem, (3.2). That is, we set n = 12,000 e,b,d and m = 1 for s = 12,000, and n = 120,000 and m = 1 for s = 120,000. In calculating the stratified sampling estimators of EPE, θˆ and θˆ , we do not address the problem of deriving u,p u,d the optimal values of n, and m ,...,m . Instead, we simply use the setting of θˆ and θˆ , 1 n e,b,p e,b,d respectively. That is, to calculate θˆ , we set n = 23, m = 524, and n = 50, m = 2433, under u,p s = 12,000 and s = 120,000, respectively. And to calculate θˆ , we set n = 12,000, m = 1, 1 2 u,d and n = 120,000, m = 1, under s = 12,000 and s = 120,000, respectively. 1 2 Table 1 on page 26 illustrates that our proposed estimators of EPE lead to substantial MSE reduction when compared to the “crude” Monte Carlo estimators for a variety of parameter settings. Comparing the MSE of the PDS-based estimators, θˆ , θˆ , and θˆ , we find that c,p e,b,p u,p our proposed stratified sampling-based estimator of EPE leads to an MSE reduction by a factor of up to 100; this unbiased estimator also dominates the efficient biased estimator of EPE, 24

in some cases quite substantially (see the third and fourth sections of Table 1). Comparing MSE of the DJS-based Monte Carlo estimators of EPE, θˆ , θˆ , and θˆ , we observe that c,d e,b,d u,d the stratified sampling-based estimator of EPE and our efficient biased EPE estimator perform similarly, which suggests that the asymptotic equivalence result in Proposition 3 can hold for even a moderate number of valuation points. Both efficient DJS estimators lead to substantial MSE reduction when compared to the corresponding crude estimator of EPE. Finally, we note that the variance and MSE for the crude estimators do not change much as the computational budget increases from 12,000 to 120,000, whereas those of efficient estimators reduce by up to an order of ten. This contrast yields the simple, yet useful insight that the number of valuation points should vary as the computational budget varies. 6.2 EEPE estimation Our numerical examples presented below illustrate the efficiency of our proposed estimators of EEPE.14 As in Section 6.1, we consider the geometric-Brownian-motion stylized example, where V t ≡ S t = S 0 eXt with X being a Brownian motion with drift µ and volatility σ. Let θˆ c,p and θˆ c,d denotethe“crude”MonteCarloestimatorsofEEPEunderPDSandDJSsampling,respectively. That is, n (cid:88) θˆ = max V¯ ∆ , (41) c,k j i 1≤j≤i i=1 where k = p,d and ∆ = t −t , and the t ’s are 1, 2, 3, 4, 8, 12, 18, 21, 24, 36, 49 weeks i i i−1 i and 1 year, with t = T = 1 year. Also, the number of simulation runs n = 1,000 and 10,000 12 respectively, for s = 12,000 and 120,000. Let θˆ and θˆ denote the efficient Monte Carlo estimators of EEPE under PDS and DJS e,p e,d sampling, respectively, based on an equidistant time grid, i.e., expression (41) with ∆ ≡ ∆ = i T/n) and resulting from solving the MSE minimization problem (38) (with constants c , i = i 1,2,3, and c therein set to 1) in Section 5.2. In particular, under s = 12,000, the optimal n = 29 and m = 414, and under s = 120,000, the optimal n = 62 and m = 1935. Our various numerical examples result in findings similar to those for the EPE estimation. For example, Table 2 on page 27, all based on 104 replications, show that the variance of the DJS-based estimators are much lower than that of the corresponding PDS-based estimators. Also, our proposed estimators of EEPE substantially outperform the crude Monte Carlo estimators in terms of MSE; for instance, MSE is reduced by a factor of 100 in the fourth section of Table 2. 7 Conclusion Ithasbecomeincreasinglycrucialforfinancialinstitutionstoactivelymanagetheircounterparty credit risk. Proper counterparty credit risk management is challenging and computationally intensive. Monte Carlo simulation is often used for CCR pricing and measurement. Poor Monte 14We refer the reader to Section (F) of the Appendix for a discussion on EEPE and numerical illustrations dst of Propositions 5 and 6. 25

EPE Variance MSE CPU Time Parameters: S = 30,µ = .2,σ = .3, s = 12,000 0 θˆ 34.6559 0.047219 0.48478 0.00380 c,p θˆ 34.6522 0.005028 0.43768 0.00162 c,d θˆ 34.1802 0.077212 0.1117 0.00253 e,b,p θˆ 33.9955 0.004785 0.004786 0.00174 e,b,d θˆ 33.9964 0.072068 0.072064 0.00518 u,p θˆ 33.9956 0.004865 0.004866 0.00335 u,d Parameters: S = 30,µ = .2,σ = .3, s = 120,000 0 θˆ 34.652 0.004791 0.4372 0.03887 c,p θˆ 34.6521 0.000501 0.43303 0.01564 c,d θˆ 34.0798 0.016741 0.024026 0.02299 e,b,p θˆ 33.9948 0.000483 0.000483 0.02409 e,b,d θˆ 33.9957 0.015533 0.015533 0.04420 u,p θˆ 33.9945 0.000486 0.000486 0.03426 u,d Parameters: S = 30,µ = 1,σ = .3, s = 12,000 0 θˆ 57.7556 0.16106 23.5389 0.00389 c,p θˆ 57.7598 0.01628 23.4351 0.00189 c,d θˆ 54.1296 0.23369 1.6954 0.00270 e,b,p θˆ 52.9238 0.015853 0.015862 0.00189 e,b,d θˆ 52.9226 0.217 0.21698 0.00516 u,p θˆ 52.9198 0.015796 0.015796 0.00390 u,d Parameters: S = 30,µ = 1,σ = .3, s = 120,000 0 θˆ 57.7579 0.016112 23.4159 0.03891 c,p θˆ 57.7591 0.001616 23.4136 0.01661 c,d θˆ 53.4783 0.047841 0.35899 0.02412 e,b,p θˆ 52.9212 0.001563 0.001564 0.02627 e,b,d θˆ 52.9189 0.045783 0.045781 0.04657 u,p θˆ 52.9203 0.001565 0.001565 0.03598 u,d Table 1: Monte Carlo EPE estimates for different parameter settings 26

EEPE Variance MSE CPU Time Parameters: S = 30,µ = 1,σ = .25, s = 12,000 0 θˆ 57.2278 0.10931 22.4936 0.00259 c,p θˆ 57.2233 0.034659 22.3768 0.00168 c,d θˆ 53.4344 0.19358 1.0731 0.00191 e,p θˆ 53.4379 0.011188 0.89722 0.00211 e,d Parameters: S = 30,µ = 1,σ = .25, s = 120,000 0 θˆ 57.2277 0.010824 22.3945 0.02866 c,p θˆ 57.2262 0.003591 22.3734 0.01427 c,d θˆ 52.9363 0.03962 0.23301 0.01817 e,p θˆ 52.9354 0.001083 0.19367 0.01545 e,d Parameters: S = 30,µ = 1.5,σ = .25, s = 12,000 0 θˆ 81.0388 .24286 101.0233 0.00279 c,p θˆ 81.0309 0.0843 100.7055 0.00173 c,d θˆ 72.8899 0.3986 3.9705 0.00221 e,p θˆ 72.8885 0.024652 3.5914 0.00226 e,d Parameters: S = 30,µ = 1.5,σ = .25, s = 120,000 0 θˆ 81.0332 0.024156 100.692 0.02929 c,p θˆ 81.0302 0.008395 100.6154 0.0144 c,d θˆ 71.8779 0.083579 0.85454 0.01935 e,p θˆ 71.8807 0.002425 0.77819 0.01630 e,d Table 2: Monte Carlo EEPE estimates for different parameter settings Carlo CCR estimation can lead to overestimation or underestimation of CCR risk. We improve the existing widely used Monte Carlo CCR frameworks by substantially increasing the efficiency of Monte Carlo estimators of the key CCR measures: EPE, CVA, and EEPE. Our proposed efficient framework can be summarized in two steps. The counterparty credit risk modeler first needs to choose between the two credit exposure sampling methods: PDS or DJS.Introducingandusingthenotionofmarginalmatching, weidentifyconditionsunderwhich the so-called path dependent simulation (PDS) method, which simulates the credit exposure process based on the finite dimensional distributions of the underlying risk factors, leads to CCR estimators whose variance is substantially larger than the variance of the CCR estimators calculated based on the so-called direct jump to simulation date (DJS) method. Taking into accountthecomputationaltimeinparallelwiththemeansquareerror,wedemonstratethatDJS sampling is preferable to PDS sampling for path independent derivatives. For path dependent derivatives since the computational time of the DJS-based estimator usually exceeds that of the PDS-based estimator, the two sampling methods could become approximately equivalent in some cases. Next, in the second step, the CCR modeler needs to choose the number of valuation points 27

andsimulationrunsateachvaluationpointforefficientEPE,CVA ,andEEPEestimations. We I show that the mean square error (MSE) of the crude Monte Carlo estimators of EPE, CVA, and EEPE can be substantially reduced by solving approximate MSE minimization problems that specify how to achieve an approximately optimal balance between bias squared and variance. These MSE minimization problems can be easily solved after approximate orders of variance and bias for each of the above mentioned CCR measures under the PDS and DJS methods are derived. For efficient EPE and CVA estimation, if the PDS method has been chosen in the first step I above, we recommend employing our unbiased stratified sampling-based estimators of EPE and CVA . These unbiased estimators of EPE and CVA use stratified sampling with the number I I of strata and simulation runs (allocated to each stratum) being chosen based on the solution to the above mentioned MSE minimization problems. Our numerical examples indicate that unbiased PDS-based estimators of EPE and CVA are preferable to the efficient biased PDS-based estimators of EPE and CVA. Our analytical results suggest that unbiased and biased DJS-based estimators of EPE and CVA are asymptotically equivalent. Our numerical examples confirm this approximate asymptotic equivalence. Finally, an interesting case arises when the CCR modeler chooses the DJS method for estimating EPE and CVA . In this case, our proposed efficient EPE and CVA estimators use 1 I I simulation run at each valuation point and the total computational budget is allocated to making the discrete time grid (the set of valuation points) as fine as possible. Our various numerical examples illustrate that employing this two-step Monte Carlo CCR framework will substantially increase the efficiency of the existing Monte Carlo CCR “engines”. Appendix A Proof of Proposition 1 We first show that cov(V ,V ) > 0 holds under condition 1. That is, V = C = B E[B−1Π |F ]. u t t t t T T t For any 0 < u < t we have cov(V ,V ) = E[cov(V ,V |F )]+cov(V ,E[V |F ]), u t u t u u t u where the last equality follows from the conditional covariance formula (see Chapter 3 of [25]). It is easy to check that the first term on the right hand side above is zero. Consider the second term and note that (cid:20) (cid:21) (cid:20) (cid:21) Π Π T T E[V |F ] = E B E[ |F ]|F = E B |F , t u t t u t u B B T T and so we conclude that for any 0 < u < t, cov(V ,V ) = B B B−2Var(E[Π |F ]) > 0. u t u t T T u 28

The second part of the proof, which is based on condition 2, uses standard results on changes of numeraire techniques and Chebyshev’s algebraic inequality (see, for instance, Proposition 2.1 in [9]). Let C denote the time-t value of a derivative contract specified in Assumption 2. Recall t Theorem 1 and Theorem 2 of [10], and note that V = C = S (t)E [(α +α Z(T))+|F ] t t 1 Q1 1 2 t where Z = S 2 /S 1 and the subscript Q1 refers to expectation under QS1, i.e. S 1 is the numeraire. Note that cov(V ,V ) = cov(V ,E[V |F ]) u t u t u = cov (cid:0) S (u)E (cid:2) (α +α Z(T))+|F (cid:3) ,E (cid:2) S (t)(α +α Z(T))+|F (cid:3)(cid:1) , 1 Q1 1 2 u Q1 1 1 2 u Consider the second term on the right hand side above, and suppose that the transition law of S accepts the following specification 1 S (t) =d βS (u)S (τ), (42) 1 1 1 where β > 0 is a constant and t−u ≡ τ (see the last part of the proof on transition law of the numeraire). We, then, have E (cid:2) S (t)(α +α Z(T))+|F (cid:3) = βS (u)E (cid:2) S (τ)(α +α Z(T))+|F (cid:3) . Q1 1 1 2 u 1 Q1 1 1 2 u This gives cov(V ,V ) = βS2(u)cov (cid:0) E (cid:2) (α +α Z(T))+|F (cid:3) ,E (cid:2) S (τ)(α +α Z(T))+|F (cid:3)(cid:1) . u t 1 QT 1 2 u QT 1 1 2 u Note that both expectations above are monotone functions of Z(u) and so using Chebyshev’s algebraic inequality (see, for instance, Proposition 2.1 in [9]) gives cov(V ,V ) > 0. u t We now show that (42) holds for the class of stochastic processes considered in this paper for modeling the dynamics of risk factors. The dynamics of risky asset S > 0 selected as the 1 numeraire in the proof above is assumed to be modeled by a GBM or a square-root diffusion. In cases where the numeraire is a maturity-T zero-coupon bond, the second part of the proof above is modified as follows. We assume that the zero-coupon bond is modeled such that it possesses an affine term structure, i.e. it has the form S 1 (t) ≡ p(t,T) = eA−Brt, where A ≡ A(t,T) and B ≡ B(t,T) are deterministic functions and the short rate r is modeled by a Gauss-Markov process or a square-root diffusion as specified before. It is, then, not difficult to see that p(t,T) = eA−Brt =d eβ1rueβ2rτ, where β 1 and β 2 are constants and 0 < u < t < T, t−u ≡ τ. Using a similar approach shown in the second part of the proof above, we arrive at cov(V ,V ) > 0. u t 29

B Proof of Proposition 2 Conditioning on FS ≡ F and using conditional covariance formula gives u u cov(V ,V ) = cov(V ,E[V |F ]) = cov(max{f(S ),0},E[max{f(S ),0}|F ]), u t u t u u t u for a well-defined function f. First consider the first term max{f(S ),0} inside the covariance u functionontherighthandsideabove. Notethatsincef isamonotonefunction,max{f(S ),0} ≡ u f˜(S ) is also a monotone function of S . Next, consider the second term E[max{f(S ),0}|F ]. u u t u Note that when S is a Gauss-Markov process, the transition law of S implies that for any 0 < u < t, S =d β S +β S , where S and S are independent random variables. Also, t 1 u 2 t−u u t−u when S is a GBM, for any 0 < u < t we have log(S ) =d log(S )+log(S ), where S and t u t−u u S are independent random variables. This follows from the independent and stationary int−u crements properties of Gauss-Markov processes15 and that their finite dimensional distributions √ are multivariate normal. When S is a square-root diffusion, dS = α(b−S )dt+σ S dB with t t t t B a standard one-dimensional Brownian motion and positive constants α and b, it can be shown that S given S is distributed as a positive constant times times a noncentral chi-square rant u dom variable with degrees of freedom that depends on α, σ, and b, and noncentrality parameter which is an increasing function of S , (see Chapter 3 of [13]). u This implies that under the class of risk factor models considered in the paper, E[max{f(S ),0}|F ] is a monotone function of S . To see this, consider the case where f t u u is an increasing function. Increasing S will increase S ; this increases max{f(S ),0}. So, u t t E[max{f(S ),0}|F ] ≡ h˜(S ) also becomes an increasing function of S . A similar argut u u u ment can be used when f is a decreasing function. Consequently, we can write cov(V ,V ) = u t cov(f˜(S ),h˜(S )), where f˜and h˜ are both either increasing or decreasing functions of S . Using u u u Chebyshev’s algebraic inequality gives cov(V ,V ) = cov(f˜(S ),h˜(S )) > 0. u t u u C Proof of Proposition 3 In this proof, for notational simplicity, we suppress the dependence of θˆ and θˆ on m b,m,n,d u,m,n,d and n. Note that, T (cid:88) n (cid:88) n (cid:90) T nMSE(θˆ ) = Var(V )∆ +n( E[V ]∆ − E[V ]dt)2, (43) b,d i i i i t m 0 i=1 i=1 where the first term on the right hand side of the above equality uses Var(V¯) = Var(V )/m. So, i i nMSE(θˆ ) converges to c (cid:82)T Var(V )dt as n → ∞. b,d 0 t Now, consider Var(θˆ ), and let I ≡ I (τ) ∈ {1,...,n} denote the index of the stratum u,d n n containing τ. Set p = P(τ ∈ A ) = ∆i. From standard results on stratified sampling we have, i i T 15Note that when S is a GBM, logarithm of S is a Brownian motion, which is a Gauss-Markov process. 30

T2 (cid:88) n T2 Var(θˆ ) = Var(V |τ ∈ A )p = E[Var(V |I )]. (44) u,d τ i i τ n mn mn i=1 (cid:82)T Since Var(V )dt = TE[Var(V |τ)], to complete the proof, it suffice to show that, as n → ∞, 0 t τ E[Var(V |I )] −→ E[Var(V |τ)]. (45) τ n τ Fromtheformulafortheconditionalvariance, toshowtheconvergencein(45), itsufficetoshow that, as n → ∞, E (cid:2) (E[V |I ])2(cid:3) −→ E (cid:2) (E[V |τ])2(cid:3) . (46) τ n τ Set X = E[V |τ] and X = E[V |I ]. Note that X is a martingale because as n increases τ n τ n n I generate increasing family of sigma-algebras. We can use martingale convergence theorem n (see Chapter 4 of [8]) to conclude that X converges to X almost surely as n → ∞. Using n continuous mapping theorem and dominated convergence theorem (see Chapter 1 of [8]) we conclude that, E[X2] converges to E[X2] almost surely, and so (46) holds. This completes the n proof of Proposition 3. 16 D Proof of Proposition 5 We first consider M . Let us assume that M = E[V ] without loss of generality. Note that, 2,m,k 2 2 max{V¯ ,V¯ }−E[V ] = V¯ 1{V¯ > V¯ }+(V¯ −E[V ])1{V¯ > V¯ }−E[V ]1{V¯ > V¯ }. (47) 1 2 2 1 1 2 2 2 2 1 2 1 2 First, consider the indicator random variable, 1{V¯ > V¯ }; the dependence of V¯ on m is sup- 1 2 i pressed for notational simplicity. Set Wk ≡ V −V , where k = d,p refer to the cases where V 1 2 1 andV areuncorrelatedandpositivelycorrelated,respectively. Notethat1{V¯ > V¯ } ≤ 1{W ¯ k > 2 1 2 E[Wk]}; Wk,...,Wk are i.i.d random variables and W ¯ k is their average. It is well known that 1 m 1{V¯ > V¯ } → 0 a.s. if and only if for all (cid:15) > 0, 1 2 P(1{V¯ > V¯ } > (cid:15) i.o.) = 0, (48) 1 2 where i.o. stands for infinitely often. To see that (48) holds, note that, P(|W ¯ k −E[Wk]| > (cid:15)˜) P(1{V¯ > V¯ } > (cid:15)) ≤ P(1{W ¯ k > E[Wk]} > (cid:15)) ≤ , (49) 1 2 (cid:15)2 for all (cid:15)˜ > 0. To derive the last inequality above the Chebyshev’s inequality is used. Then, (49), almost sure convergence of W ¯ k → E[Wk] following from the strong law of large numbers (SLLN), and Kolmogorov’s 0-1 law, (see Theorem 8.1 of [8]), give (48). 16The probabilistic arguments used in second part of the proof are similar to the ones used in the proof of Lemma 4.1 in [14]. 31

EPE Variance CPU Time Parameters: a = 0.01,b = 0.0004,σ = 0.3,r = 0.02 0 θˆ 36.0561 0.001335 0.00953 b,m,n,p θˆ 36.0553 0.000196 0.01018 b,m,n,d Parameters: a = 0.01,b = 0.0002,σ = 0.3,r = 0.04 0 θˆ 35.9658 0.001416 0.01012 b,m,n,p θˆ 35.9638 0.000206 0.01146 b,m,n,d Parameters: a = 0.01,b = 0.0003,σ = 0.4,r = 0.03 0 θˆ 37.009 0.002546 0.00909 b,m,n,p θˆ 37.012 0.000380 0.01025 b,m,n,d Table 3: Monte Carlo EPE estimates for a single interest rate swap with different parameter settings Now, consider the first term on the right side of (47). Given that V¯ and 1{V¯ > V¯ }, almost 1 1 2 surely converge to E[V ] and zero, respectively, it is not difficult to show that 1 V¯ 1{V¯ > V¯ } → 0 a.s.. 1 1 2 To see this, it suffices to write V¯ 1{V¯ > V¯ } = (V¯ −E[V ])1{V¯ > V¯ }+E[V ]1{V¯ > V¯ }, 1 1 2 1 1 1 2 1 1 2 and use SLLN for the sequence of indicator random variables and V¯ . The last term on the right 1 side of (47) converges to zero a.s. based on (48). Analogous arguments led to (48) show that the second term on the right side of (47) converges to zero a.s. This completes the proof for n = 2. Induction and analogous arguments are employed for the general case. Suppose that M → M , a.s. as m → ∞. Assume that M = E[V ]. Then, we need n−1,m,k n−1 n n to show that a similar almost sure convergence holds for M . To see this, it suffices to note n,m,k that, for all (cid:15) > 0 and (cid:15)˜> 0, P (cid:0) (V¯ −E[V ])−(max{V¯ ,...,V¯ }−E[V ]) > (cid:15)˜ (cid:1) P(1{V¯ > max{V¯ ,...,V¯ }} > (cid:15)) ≤ 1 1 2 n n 1 2 n (cid:15)2 P (cid:0) |V¯ −E[V ]| > (cid:15)˜ (cid:1) P (cid:0) |max{V¯ ,...,V¯ }−E[V ]| > (cid:15)˜ (cid:1) 1 1 2 n n ≤ + , (cid:15)2 (cid:15)2 which is then used to show that 1{V¯ > max{V¯ ,...,V¯ }} −→ 0 a.s. 1 2 n This completes the proof. 32

E Numerical Examples for Interest Rate Swaps Table 3 on the preceding page above compares the variance and computing time of the DJS and PDS-basedMonteCarloestimatorsofEPEforasingleinterestrateswapcontractintheVasicek shortrateframework(see,forinstance,[2]). Specifically,theshortrate,denotedbyr,ismodeled by dr = (b−ar )dt+σdB , where B is a standard Brownian motion and a > 0. Recall that t t t standard results from affine term structure modeling expresses the time-t value of a T-maturity zero coupon bond as p(t,T) = eA(t,T)−B(t,T)rt, where A and B are given deterministic functions of time, a, b, and σ. The interest rate swap in our examples is a forward swap settled in arrears (see, for instance, [2]) with quarterly payments, a principal value of 100, and maturity 1 year. We let m = 104 and n = 12. In particular, the 12 valuation points are equally spaced, i.e., monthly, within the one-year interval. Variances of the DJS and PDS-based estimators θˆ b,m,n,d and θˆ (column 2 of Table 3 on the previous page) are estimated based on 1000 simulation b,m,n,p runs. Table 3 presents instances of our numerical examples for interest rate swap where the DJS-based estimator of EPE outperforms the PDS-based estimator in terms of variance (by an order of at least 10) while the computing times are approximately equal. F Numerical Examples for EEPE estimation dst The numerical results presented in this section demonstrate the consistency of PDS and DJS estimators for EEPE and the asymptotic efficiency of DJS over PDS. In particular, they dst support our Propositions 5 and 6. Similar to our numerical examples of Sections 3 and 5, we set V ≡ S with S being a geometric Brownian motion with initial value S = 30, drift µ = 0.01, and volatility σ = 1 here. 0 We compare the crude PDS and DJS estimators θˆ and θˆ as defined in Section 6.2. Each c,p c,d estimation procedure is replicated 10,000 times to produce the estimates. In Table 4 on the following page, in addition to presenting the estimator value, variance, MSE (which, unlike that in Section 6.2, is defined with respect to the estimand of EEPE ), dst and CPU time, we also include a column named “WrongOrderProb”, which gives the estimate for the probability that the indices at which the running maximums are achieved ever go wrong, i.e., using the notation in Section 5.1 of the main paper, P(τ (cid:54)= τ ,for some i), k = p or i,m,k i d, corresponding to PDS and DJS respectively. The sum of these two probabilities provides an upper bound for P(Bc ) in the statement of Proposition 6. As this table shows, this upper m bound converges to zero as m increases, which implies lim P(Bc ) = 0. Also, the bias of m→∞ m both estimators vanishes as m increases; this is consistent with Proposition 5. G The DJS Method This section gives an explanatory discussion of what we have referred to as the “DJS method” throughoutourstudy. TheunderlyingideaoftheDJSmethodisrathersimple. Supposethatthe twogivenrandomvariablesV andV aredependent;MonteCarloconsidersthejointdistribution i j 33

EEPE Variance MSE CPU Time WrongOrderProb dst Number of simulation runs m = 50 θˆ 41.0514 20.6228 20.7011 0.000329 0.9613 c,p θˆ 42.4019 7.901 10.5697 0.000431 1 c,d Number of simulation runs m = 500 θˆ 40.7481 2.1847 2.1849 0.00156 0.2789 c,p θˆ 40.8833 0.77441 0.78761 0.00121 0.919 c,d Number of simulation runs m = 5000 θˆ 40.7756 0.21446 0.2146 0.0173 0 c,p θˆ 40.7708 0.07379 0.07379 0.0090 0.1987 c,d Number of simulation runs m = 50000 θˆ 40.7708 0.02156 0.02156 0.1708 0 c,p θˆ 40.7682 0.00711 0.00711 0.0932 0.0001 c,d Table 4: Monte Carlo EEPE estimates for based on increasing number of simulation runs dst function of these two random variables to generate copies of V and V . This is what the paper i j refers to as the “PDS method”. Since the functional form of the CCR measures considered in our study depends only on the marginal distributions of V and V as characterized by the i j notion of Marginal Matching, the dependence structure of V and V , specified by their joint i j distribution function, need not be preserved when Monte Carlo is used to estimate these classes of CCR measures. That is, V and V can be sampled from independently and so cov(V ,V ) = 0. i j i j This is what the paper refers to as the “DJS method”. In the context of CCR, V , representing i the exposure process evaluated that the i’th time point on the grid, could be a function of m dependent random variables, S ,...,S and V could be a function of n dependent random 1 m j variablesS ,...,S ,...,S . IngeneratingV ,whilethePDSmethodusesthesameS ,...,S that 1 m n j 1 m have been used in simulating V , the DJS method, starting anew, uses new copies of S ,...,S i 1 m that are independent from the previous ones. This is an additional computational burden, and so in cases where the DJS method proves to be statistically more efficient than the PDS method in terms of the variance of the Monte Carlo CCR estimators, (by using the analytical results of Proposition 1 and 2 or by conducting a small simulation study), computing time should also be taken in to account. The useful criterion that does this is, variance per replication×expected computing time. This criterion, when used even as a rule of thumb, would specify whether to choose the DJS or the PDS method. The DJS method in the CCR literature, to the best of our knowledge, often only refers to cases where the computing time of the DJS-based estimator is equivalent to that of the PDSbased estimator. The Path-Independent Case of Section 3.1.1 is the classical situation: V , in i the above stylized example, is a function of a single random variable S , and V is a function of i j onlyS , whereS andS aredependentwithagivenjointdistributionfunction. Whatourstudy j i j 34

refers to as the “DJS method”, as illustrated above, would, then, be an extension of the DJS in its conventional sense. While we have not changed the well established CCR terminologies, “DJS and PDS”, we invite the reader to be mindful of the distinction between the DJS in its traditional sense and in the more general sense of this paper which could require additional computing time. H Time-Discretized SDEs and Monte Carlo-Based Pricing This section discusses cases that have not been formulated in the main body of the paper. We discussMonteCarloCCR measurementwhen theunderlyingriskfactors aremodeledby continuoustimestochasticdifferentialequations(SDEs)whosesimulationrequirestime-discretization. We also discuss Monte Carlo CCR when in addition to generating the underlying risk factors up to each valuation point, Monte Carlo is also used for pricing the portfolio constituents. The potential applicability of the so called DJS method in this paper’s setting is further clarified for these more general cases. Time-discretized SDEs and Monte Carlo-based pricing create additional sources of bias (for the estimators of the CCR measures) that have not been mathematically formulated in our study. This could be the subject of future research. We are to show and emphasize that the proposed framework does lead to basic Monte Carlo CCR efficiency improvements in these cases; the proposed framework assumes that the factors influencing these additional sources of bias are held constant and then specifies how the remaining factors determiningthestatisticalefficiencyandcomputingtimeoftheCCRestimatorsshouldbedealtwith. The simple stylized example below is used throughout this section to facilitate our discussion of the above-mentioned topics. Example: Arithmetic Asian Options Let{S ;t ≥ 0}denoteageometricBrownianmotion t satisfying the following SDE, dS = µS dt+σS dB , (50) t t t t where {B ;t ≥ 0} denotes a one-dimensional standard Brownian motion, and the drift µ and t volatility σ parameters are given constants. Consider the equi-distant time grid 0 ≡ t < t < 0 1 ... < t ≡ T, t −t ≡ h, i = 1,...,n. Let V denote the time-t ≥ 0 value of an arithmetic n i i−1 t Asian option with time-T payoff (S¯−K)+ for a given strike K, where S¯ = (S +...+S )/n 1 n amd S ≡ S . For simplicity assume that risk-free rate is zero, and so i ti V = E[(S¯−K)+|F ], i i where F = F denote the filtration generated by the GBM up to time t . We are to use Monte i ti i (cid:82)T Carlo to generate V ’s and estimate θ = E[V ]dt. i 0 t Time-Discretized SDEs Suppose that we cannot analytically solve the SDE in (50). That is, suppose that the GBM could not have been simulated exactly (the joint distribution of the simulated values would not coincide with that of the continuous time model at the simulated 35

datesonthetimegrid),andanEulertime-discretizationschemewouldhavetobeusedtosample from the SDE. Let Sˆ denote a time-discretized approximation to S on the above mentioned time grid,17 √ Sˆ = Sˆ +µSˆh+σSˆ hZ , (51) i+1 i i i i+1 with Z ,Z ,... independent standard normal random variables. Consider t and t on the time 1 2 j l grid where 0 < t < t < T. Generating Sˆ with the Euler scheme using (51) requires j i.i.d j l j standard normals Z ,...,Z . Now, consider the simulation of Sˆ. Consider two cases. If the set 1 j l of normals Z ,...,Z generated in calculating Sˆ are used again along with “new” i.i.d. standard 1 j j normal random variables Z ,...Z to generate Sˆ, then Sˆ and Sˆ become positively correlated j+1 l l j l according to (51). If, however, l new i.i.d. standard normal random variables Z ,...,Z are 1 l generated to simulate Sˆ, then, clearly, Sˆ and Sˆ become independent. The price of achieving, l j l cov(Sˆ ,Sˆ) = 0, istheadditionalcomputationaltimetogeneratej newstandardnormalrandom j l variable instead of using the old ones that have already been generated when simulating Sˆ . j Time-Discretized SDEs and the DJS Method Continuing to assume the Euler timediscretization is used to sample from the underlying GBM, we now return to the Asian option example. View V as a function of (S ,...,S ), i.e., V = E[(S¯−K)+|F ] = f(S ,...,S ). Assume i 1 i i i 1 i thatgivenF ,theconditionalexpectationE[(S¯−K)+|F ]canbecalculatedanalytically,(laterin i i this section, we discuss the implications of relaxing this assumption). Consider two Monte Carlo realizations of V and V . Similar to our analysis of the Path Dependent case of Section 3.1.2, j l V = f(S ,...,S )andV = f(S ,...,S ,...,S )becomeindependentifwestartanewingenerating j 1 j l 1 j l V . More specifically, assuming the Euler scheme in (51) is in place, if the set of Z ,Z ,..,Z l 1 2 l generatedincalculatingVˆ = f(Sˆ ,...,Sˆ)areindependentfromthesetofZ ,Z ,...,Z generated l 1 l 1 2 j in sampling from Vˆ = f(Sˆ ,...,Sˆ ), Vˆ and Vˆ become independent. We emphasize that it is j 1 j j l only in this sense that what we continue to refer to as the DJS method is relevant to our study; that is, in sampling from the underlying risk factors in a way that V’s evaluated at any two distinct time points on the grid become independent. Recall the notion of Marginal Matching and note that in sampling from V = f(S ,...,S ), j 1 j making Sˆ ,...,Sˆ independent from each other is neither relevant to our objective - having 1 j cov(Vˆ ,Vˆ) = 0 - nor possible in the sense that, depending on the functional form of f, the j l expectedvalueofVˆ underindependentSˆ ,...,Sˆ maynotbecomeequaltotheexpectedvalueof j 1 j Vˆ = f(Sˆ ,...,Sˆ ) when Sˆ ,...,Sˆ are sequentially sampled from based on the finite dimensional j 1 j 1 j distribution of the GBM time-discretized in (51). That is, our analysis and results are not affected by the inability to use a DJS type method in sampling from V = f(S ,...,S ) for a j 1 j given fixed j. Again, what the DJS method, whose applicability is specified by the notion of Marginal Matching, is concerned with in our study is making Monte Carlo realizations of V and j V independent - whether they have been simulated exactly or via an Euler time-discretization l 17To simplify our discussion of time-discretized SDEs, we have assumed that the CCR time grid containing the valuations points and the time grid based on which the SDEs are time-discretized coincide. Relaxing this assumption does not change the results of this discussion section. 36

scheme.18 Time-Discretized SDEs and the Additional Source of Bias Two types of bias have been formulated in the MSE minimization problems of our study: the time-discretization bias associated with the time integral of the exposure process in the definition of EPE, CVA, EEPE and also the bias induced by the maximum operator in the definition of the EEPE. The Eulerscheme-based time-discretized SDE, Sˆ, which approximates S, creates an additional source of bias that has not been mathematically formulated in our study. A typical criterion based on which the quality of the SDEs’ time discretization methods is measured is the convergence order of the form |E[g(Sˆ)]−E[g(S )]| for a given i with a given function g usually satisfying some i i smoothness conditions. See Chapter 6 of [13] and the references therein for typical weak and strong error criteria and for methods to improve the convergence order of Euler schemes based on these criteria. The simple Euler scheme (51) converges with respect to the above-mentioned criterion as the time step h decreases to zero. The proposed efficient CCR framework assumes that an Euler scheme (possibly an improved version of it by using the methods introduced in Chapter 6 of [13]) is already in place with a fixed time step h, and it then achieves an approximately optimal balance between the number of simulation runs and valuation points under the chosen Euler scheme. Monte Carlo-Based Pricing We now discuss the implications of Monte Carlo-based pricing in our setting. Suppose that given F , the conditional expectation E[(S¯−K)+|F ] could not be i i calculated analytically. Also, suppose that, motivated by the Longstaff and Schwartz American optionpricingappraoch[21], American(orLeast-Squares)MonteCarloistobeusedtoestimate θ. In the m-simulation-run-based commonly-used version of American Monte Carlo for CCR, one often proceeds as follows (see [6]).19 First generate the m sample paths on the grid 0 ≡ t < t < ... < t ≡ T, i.e., {S ,S ,...S }, k = 1,...,m, based on the finite dimensional 0 1 n k1 k2 kn distribution of the process S. Then, all Monte Carlo estimators of E[V ],...,E[V ] are estimated 1 n based on this single set of sample path data. To generate Monte Carlo realization of V = E[(S¯ − K)+|F ] = f(S ,...,S ), one first i i 1 i uses a regression-based method to approximate f with f˜ and uses all the simulated data {S ,S ,...S }, k = 1,...,m to estimate the regression coefficients. Then, using a given samk1 k2 kn pled path of S ,...,S , to generate a single approximate realization of V , which we denote by V˜, 1 i i i 18Importantclassesofmodelswhoseapplicationrequirestime-discretizationoftheunderlyingSDEsareLIBOR andSwapMarketmodels,(see,e.g.,Chapter25of[2],Chapter3of[13],andthereferencestherein). Forinstance, consider pricing Caplets in the LIBOR model where the LIBOR forward rates are GBM driven, (see Definition 25.5 of [2]), where the simulation of SDEs characterizing the underlying LIBOR forward rates (see Proposition 25.7of[2])requireEulertime-discretizationschemes,(seeSection25.5of[2]andSection3.7of[13]). Asillustrated above, at the cost of additional computing time compared with that of the PDS method, the DJS method can be applied in the LIBOR market models to make the simulated V’s, representing the value of Caps or Caplets, independent at any two distinct time points on the CCR grid. That is, the so called DJS method is potentially applicable. However, the CCR modeler, considering both the computing time and the variance of the PDS and DJS-based estimators, may conclude that the PDS method should be chosen. 19We assume that the reader is familiar with American Monte Carlo option pricing and its application in counterparty credit risk. 37

one sets V˜ = f˜(S ,...,S ). The two random variables V˜ = f˜(S ,...,S ) and V = f(S ,...,S ) are i 1 i i 1 i i 1 i notequalindistribution; itisinthissensethatwerefertoV˜ asanapproximaterealizationofV . i i Note that this regression-based approximation of f with f˜makes American Monte Carlo based estimators of θ biased. This bias has not been mathematically formulated in our framework. As mentioned above, the commonly-used version of American Monte Carlo in CCR uses the same single set of simulated sample paths in estimation of E[V ],...,E[V ]. To make American 1 n Monte Carlo estimators of E[V ],...,E[V ] uncorrelated, using the notion of Marginal Matching, 1 n each American Monte Carlo estimator of E[V ] should use a ‘new’ independent set of simulated i sample paths. More specifically, fix i, and suppose that E[V ] has been estimated based on i {S ,S ,...S }, k = 1,...,m as we summarized above. Now, using the notion of Marginal k1 k2 kn Matching (or what we continue to refer to as the DJS method), to estimate E[V ], instead of i+1 using the same set of sampled paths employed in the estimation of E[V ], a new independent set i {S ,S ,...S }, k = 1,...,m is being generated to simulate V˜ ’s. So, under the DJS method, k1 k2 kn i+1 American Monte Carlo realizations of V˜ and V˜ become independent and cov(V˜,V˜ ) = 0. i i+1 i i+1 That is, at the cost of additional computational time, the DJS method can be used to make American Monte Carlo-based realizations of V ’s independent.20 i Convergence Results of American Monte Carlo Algorithms It is well known that formulating and deriving precise convergence results on American option pricing algorithms, e.g., the Least Squares method of Longstaff and Schwartz, is quite difficult, (see Section 2.3 of [21] and Chapter 8 and Section 8.6 of [13]). Similar challenges exist when American Monte Carlo is used to estimate CCR measures since we need to consider limits as the number of valuation points (time points on the grid), the number of regression-based basis functions, and thenumberofsimulationpathsgotoinfinity. WewouldliketoemphasizethatwhenanAmerican CCR “engine” is in place, even in the absence of the formulation and optimization of the MSE of the CCR estimators in the presence of all sources of bias, our framework would lead to basic efficiency improvements. It characterizes and differentiates the DJS and PDS methods by considering computational time in parallel with variance and then achieves an approximate optimal balance between the number of simulation runs (the number of sampled paths) and the number of valuation points assuming that the functional form and the number of basis functions used in the regression have been fixed. I Estimating CVA Sensitivities We consider Monte Carlo estimation of the sensitivity of independent CVA with respect to the initial value of an underlying risk factor, (section 3 of [19] discusses this type of CVA sensitivity estimation with crude Monte Carlo; see Chapter 16 of [18] for a more complete overview of CVA ‘Greeks’ for hedging counterparty risk). We merely focus on comparing the DJS and PDS methods for CVA sensitivity estimation with respect to an underlying risk factor’s initial value. 20SimilarargumentsholdfornestedMonteCarlo. Thatis,thesecondlayerofMonteCarlousedforpricingthe contract at each valuation point makes the estimator of θ biased. And, the DJS method in the above mentioned sense can still be used to make V’s independent at any two distinct time points on the grid. 38

Let S denote an underlying risk factor’s initial value and consider Monte Carlo estimation of 0 dCVA I,whereCVA = (cid:82)T E[V ]dF ,usingafinite-differenceapproximationmethod. TheMonte dS0 I 0 t t Carlo estimator of CVA uses the equidistant discrete time grid, 0 ≡ t < t < ... < t ≡ T. I 0 1 n The central-difference estimator of dCVA I based on m simulation runs at each valuation point dS0 is, n (cid:88) θˆ = DˆV ∆F , (52) k i,m i i=1 where V¯(S +h)−V¯(S −h) DˆV = i 0 i 0 , (53) i,m 2h where V¯(S +h) is the m-simulation-run average of V ,...,V with the underlying initial value i 0 i1 im at S +h for some h > 0. It is well-known that DˆV is a biased estimator of dE[Vti ] . It is 0 i,m dS0 also well-known that using common random numbers in estimating V¯(S +h) and V¯(S −h) i 0 i 0 will reduce the order of variance of the finite difference estimator. We assume that common random numbers are used. That is, suppose that V is a function of S , i.e., V ≡ f(S ). Then, i i i i the same set of generated S ,...,S used in the simulation of V¯(S +h) are also being used i1 im i 0 in the simulation V¯(S −h).21 The subscript k in θˆ takes d and p when the DJS and PDS i 0 k methods are used, respectively. Recall that under DJS and the notion of Marginal Matching, any two distinct pairs of S ,...,S areuncorrelated,andsothevarianceofθˆ isequaltotheweightedsumofthevariances 1 n d of DˆV ≡ DˆV (S ), i = 1,...,n, and the covariance terms, cov(DˆV (S ),DˆV (S )), are zero i,m m i m i m j foranyi (cid:54)= j. So,tocomparethevarianceofθˆ andθˆ itremainstoconsiderthecovarianceterms d p under the PDS method. If using PDS results in positive covariance terms, then the DJS based estimator of dCVA I outperform the PDS-baed estimator in terms of variance. Propositions 1 dS0 and 2 have identified conditions under which the covariance function cov(V ,V ) is positive for u t contract level exposure for any 0 < u < t < T. The same method of proof can be employed to show that when using PDS, under certain assumptions cov(DˆV (S ),DˆV (S )) > 0 for any m i m j i (cid:54)= j. We give the proof here only for a simple stylized case. Using the terminology and notation developed in Section 3, consider the contract level credit exposure for a unilateral transaction, i.e., V = max{C ,0} = C for all 0 < t ≤ T, where the t t t time-t risk neutral value of the maturity-T derivatives contract is given by C = n E[ΠT|F ] t t nT t with n being the numeraire, F denoting the time-t filtration generated by the underlying risk t factors, and Π being the payoff function. Consider the setting where the numeraire is the money market account, B, with deterministic short rate. Also, consider the single-simulation-run case, i.e., m = 1, where we set DˆV ≡ DˆV . We are to show that, i,1 i cov(DˆV ,DˆV ) > 0 i j 21See Chapter 7 of [13] and the reference therein for finite-difference approximations. More specifically, see Section 7.1 on how the MSE of a finite-difference estimator can be approximately optimized. 39

for given i < j on the discrete time grid. Using the conditional covariance formula as in the proofs of Propositions 1 and 2 gives, cov(DˆV ,DˆV ) = cov(DˆV ,E[DˆV |F ]). i j i j i Note that V (S + h) = B B−1E [Π |F ] and V (S − h) = B B−1E [Π |F ], where the j 0 j n h+ T j j 0 j n h− T j subscript h and h denote the cases where the underlying risk factor’s initial prices are set at + − S +h and S −h, respectively. Then, we have, 0 0 (cid:18) (cid:19) E [Π |F ]−E [Π |F ] E[DˆV |F ] = B B−1 h+ T i h− T i , j i j n 2h and so, (cid:18) (cid:19) E [Π |F ]−E [Π |F ] cov(DˆV ,DˆV ) = B B B−2Var h+ T i h− T i > 0. i j i j n 2h Note that more general cases, e.g., stochastic short rate and bilateral transactions, require the application of changes of numeraire techniques and the Chebyshev’s algebraic inequality and constraints on the number of risk factors and monotonicity of the payoff function as in the proofs of Proposition 1 and 2; we do not address these cases here as the proof methodology has alreadybeendevelopedinthepaper. Themathematicalproofsofthepositivityofthecovariance function under PDS are at the contract level exposure. For more general cases, comparing the varianceofthePDSandDJS-basedCVAsensitivityestimatorsrequiressmallsimulationstudies before the CCR modeler makes a decision on whether to use DJS or PDS. References [1] “Capital requirements for bank exposures to central counterparties”, Basel Committee on Banking Supervision, July (2012) [2] T. Bjork, “Arbitrage theory in continuous time”, 3rd Edition, Oxford University Press, (2009) [3] Bohme, M., Chiarella, D., Harle, P., Neukirchen, M., Poppensieker, T., and Raufuss, A., “Day of reckoning: new regulation and its impact on capital market businesses”, McKinsey Working Papaers on Risk, (September 2011) [4] Broadie, M., Du, Y., Moallemi, C. C., “Risk estimation via regression”, Working paper, December (2011) [5] Canabarro, E., and Duffie, D., “Measuring and Marking Counterparty Risk”, Asset/Liability Management for Financial Institutions, Institutional Investor Books, (2003) 40

[6] Cesari, G., Aquilina, J., Charpillon, N., Filipovic, Z., Lee, G., Manda, I, “Modelling, Pricing, and Hedging Counterparty Credit Exposure”, Springer Finance, (2010) [7] Duffie, D., and Glynn, P., “Efficient Monte Carlo estimation of security prices”, Annals of Applied Porbability, vol. 5, pp. 897-905, (1995) [8] Durrett, R., “Probability: Theory and Examples”, 3rd Edition, Duxbury Advanced Series, (2005) [9] Egozcue, M., Garcia, F., L., Wong, W. K., “On some covariance inequalities for monotonic and non-monotonic functions” Journal of Inequalities in Pure and Applied Mathematics, vol.10, (2009) [10] Geman, H., El Karoui, N., Rochet., J. C., “Changes of Numeraire, changes of probability measure and option pricing”, Journal of Applied Probability, vol. 32, pp. 443-458, (1995) [11] Ghamami, S., “Book Review: Counterparty Credit Risk by Jon Gregory,” Quantitative Finance, vol. 13, no. 12, pp. 1863-1865, (2013) [12] Ghamami, S., Goldberg, L. R., “Stochastic intensity models of wrong way risk: wrong way CVA need not exceed independent CVA”, Journal of Derivatives, vol. 21, no. 3, pp. 24-35, (2014) [13] Glasserman, P., “Monte Carlo Methods in Financial Engineering”, Springer, (2004) [14] Glasserman, P., Heidelberger, P., and Shahabuddin, P., “Asymptotically optimal importance sampling and stratification for path dependent options”, Mathematical Finance, vol. 9, pp. 117-152, (1999) [15] Glynn, P. W., and Whitt, W., “The efficiency of simulation estimators”, Operations Researcg, vol. 40, pp. 505-520, (1992) [16] Gordy, M. B., and S. Juneja, “Nested Simulation in Portfolio Risk Measurement”, Management Science, vol. 56, no. 10, pp. 1833-1848, (2010) [17] Gregory, J., “Counterparty credit risk: the new challenge for global financial markets”, 1st edition, John Wiley and Sons, (2010) [18] Gregory, J., “Counterparty credit risk: a continuing challenge for global financial markets”, 2nd edition, John Wiley and Sons, (2012) [19] Hull, J., and White, A., “CVA and wrong way risk”, Financial Analyst Journal, vol. 68, no. 5, pp. 58-69, (2012) [20] Karatzas, I., and Shreve, S. E. , “Brownian Motion and Stochastic Calculus”, Springer, (1991) 41

[21] Longstaff, F. A., and Schwartz, E. S., “Valuing American options by simulation: A simple Least-Squares approach”, The Review of Financial Studies, vol. 14, pp. 113-147, (2001) [22] Musiela,M.,andRutkowski,M.,“MartingaleMethodsinFinancialModelling”,2ndEdtion, Springer, (2010) [23] Pykhtin, M., and Zhu, S., “A guide to modeling counterparty credit risk”, GARP Risk Review, July/August, 16-22, (2007) [24] Pykhtin, M., and Zhu, S., “Measuring counterparty credit risk for trading products under Basel II”, In the Basel Handbook, 2nd edition, Risk books, London, (2006) [25] Ross, S. M., “Introduction to Probability Models”, Academic Press, 10th edition, (2009) [26] Ross, S. M., “Simulation”, 4th edition, Academic Press, (2006) [27] Xing, J., Fu, M., Xiong, X., “Probabilisticerrorboundsforsimulationquantileestimators”, Management Science, vol. 14, no. 2, pp. 230-246, (2003) 42