Efficient Computation of Various Valuation Adjustments Under Local Lévy Models
EEfficient Computation of Various Valuation Adjustments UnderLocal L´evy Models
Anastasia Borovykh ∗ Andrea Pascucci † Cornelis W. Oosterlee ‡ May 7, 2019
Abstract
Various valuation adjustments, or XVAs, can be written in terms of non-linear PIDEs equivalent toFBSDEs. In this paper we develop a Fourier-based method for solving FBSDEs in order to efficiently andaccurately price Bermudan derivatives, including options and swaptions, with XVA under the flexibledynamics of a local L´evy model: this framework includes a local volatility function and a local jumpmeasure. Due to the unavailability of the characteristic function for such processes, we use an asymptoticapproximation based on the adjoint formulation of the problem.
Keywords:
Fast Fourier Transform, CVA, XVA, BSDE, characteristic function
After the financial crisis in 2007, it was recognized that Counterparty Credit Risk (CCR) poses a substantialrisk for financial institutions. In 2010 in the Basel III framework an additional capital charge requirement,called Credit Valuation Adjustment (CVA), was introduced to cover the risk of losses on a counterpartydefault event for over-the-counter (OTC) uncollateralized derivatives. The CVA is the expected loss arisingfrom a default by the counterparty and can be defined as the difference between the risky value and thecurrent risk-free value of a derivatives contract. CVA is calculated and hedged in the same way as derivativesby many banks, therefore having efficient ways of calculating the value and the Greeks of these adjustmentsis important.One common way of pricing CVA is to use the concept of expected exposure, defined as the mean ofthe exposure distribution at a future date. Calculating these exposures typically involve computationallytime-consuming Monte Carlo procedures, like nested Monte Carlo schemes or the more efficient least squaresMonte Carlo method (LSM)[17]. Recently the Stochastic Grid Bundling method (SGBM)[14] was introducedas an improvement of the standard LSM. This method was extended to pricing CVA for Bermudan options in ∗ Dipartimento di Matematica, Universit`a di Bologna, Bologna, Italy. ( e-mail :[email protected]) † Dipartimento di Matematica, Universit`a di Bologna, Bologna, Italy. ( e-mail :[email protected]) ‡ Centrum Wiskunde & Informatica, Amsterdam, The Netherlands and Delft University of Technology, Delft, The Nether-lands. ( e-mail :[email protected]) a r X i v : . [ q -f i n . M F ] M a y B and its counterparty C in the presence of CCR, bilateral collateralization, MVA, FVA and KVA, by setting up a hedging portfolioin which we focus on hedging the default risks and take into account the different rates associated withdifferent types of lending. We extend the Fourier-based method known as the BCOS method, developedin [23], to solve the BSDE under L´evy models with non-constant coefficients. As this method requires the2nowledge of the characteristic function of the forward process, which, in the case of the L´evy process withvariable coefficients, is not known, we will use an approximation of the characteristic function obtained bythe adjoint expansion method developed in [19], [18] and extended to the defaultable L´evy process with astate-dependent jump measure in [2]. Compared to other state-of-the-art methods for calculating XVAs, likeMonte Carlo methods and PDE solvers, our method is more efficient and/or flexible. The efficiency is bothdue to the availability of the characteristic function in closed form through the adjoint expansion methodand the fast convergence of the COS method. Furthermore we propose an alternative Fourier-based methodfor explicitly pricing the CVA term in case of unilateral CCR for Bermudan derivatives under the local L´evymodel. The advantage of this method is that is allows us to use the FFT, resulting in a fast and efficientcalculation. The Greeks, used for hedging CVA, can be computed at almost no additional cost.The rest of the paper is structured as follows. In Section 2 we introduce the L´evy models with non-constant coefficients. In Section 3 we derive the non-linear PIDE and corresponding BSDE for pricingcontracts under XVA. In Section 4 we propose the Fourier-based method for solving this BSDE and inSection 5.1 this method is extended to pricing Bermudan contracts. In Section 5.2 an alternative FFT-based method for pricing and hedging the CVA term is proposed and Section 6 presents numerical examplesvalidating the accuracy and efficiency of the proposed methods. We consider a defaultable asset S t whose risk-neutral dynamics are given by S t = { t<ζ } e X t ,dX t = µ ( t, X t ) dt + σ ( t, X t ) dW t + (cid:90) R qd ˜ N t ( t, X t − , dq ) ,d ˜ N t ( t, X t − , dq ) = dN t ( t, X t − , dq ) − a ( t, X t − ) ν ( dq ) dt, (1) ζ = inf { t ≥ (cid:90) t γ ( s, X s ) ds ≥ ε } , where d ˜ N t ( t, X t − , dq ) is a compensated random measure with state-dependent L´evy measure ν ( t, X t − , dq ) = a ( t, X t − ) ν ( dq ) . The default time ζ of S t is defined in a canonical way as the first arrival time of a doubly stochastic Poissonprocess with local intensity function γ ( t, x ) ≥
0, and ε ∼ Exp(1) and is independent of X t . This way ofmodeling default is also considered in a diffusive setting in [5] and for exponential L´evy models in [4]. Thus,our model includes a local volatility function, a local jump measure, and a default probability which isdependent on the underlying. We define the filtration at time t of the market observer to be G t = F Xt ∨ F Dt ,where F Xt is the filtration generated by X upto time t and F Dt := σ ( { ζ ≤ u } , u ≤ t ), for t ≥
0, is thefiltration of the default. Using this definition of default, the probability of default isPD( t ) := P ( ζ ≤ t ) = 1 − E (cid:104) e − (cid:82) t γ ( s,X s ) ds (cid:105) . (2)3e assume furthermore (cid:90) R e | q | a ( t, x ) ν ( dq ) < ∞ . Imposing that the discounted asset price ˜ S t := e − rt S t is a G -martingale under the risk-neutral measure, weget the following restriction on the drift coefficient: µ ( t, x ) = γ ( t, x ) + r − σ ( t, x )2 − a ( t, x ) (cid:90) R ν ( dq )( e q − − q ) , (3)with r being the risk-free (collateralized) rate. In the whole of the paper we assume deterministic, constantinterest rates, while the derivations can easily be extended to time-dependent rates. The integro-differentialoperator of the process is given by (see e.g. [20]) Lu ( t, x ) = ∂ t u ( t, x ) + µ ( t, x ) ∂ x u ( t, x ) − γ ( t, x ) u ( t, x ) + σ ( t, x )2 ∂ xx u ( t, x )+ a ( t, x ) (cid:90) R ν ( dq )( u ( t, x + q ) − u ( t, x ) − q∂ x u ( t, x )) . (4) Consider a bank B and its counterparty C , both of them might default. Assume they enter into a contractpaying Φ( S t ) at maturity. Let φ ( x ) = Φ( e x ), and assume the risk-neutral dynamics of the underlying as in(1) with the drift given by (3). Define ˆ u ( t, x ) to be the value to the bank of the (default risky) portfoliowith valuation adjustments referred to as XVA and u ( t, x ) to be the risk-free value. Note that the differencebetween these two values is called the total valuation adjustment and in our setting this consists ofTVA := ˆ u ( t, x ) − u ( t, x ) = CVA + DVA + KVA + MVA + FVA . (5)The risk-free value u ( t, x ) solves a linear PIDE: Lu ( t, x ) = ru ( t, x ) ,u ( T, x ) = φ ( x ) , where L is given in (4). Assuming the dynamics in (1), this linear PIDE can be solved with the methodspresented in [2]. In [3], the authors derive an extension to the Black-Scholes PDE in the presence of a bilateral counterpartyrisk in a jump-to-default model with the underlying being a diffusion, using replication arguments that includethe funding costs. In [15] this derivation is extended to a multivariate diffusion setting with stochastic ratesin the presence of CCR, assuming that both parties B and C are subject to default. To mitigate the CCR,both parties exchange collateral consisting of the initial margin and the variation margin. The parties areobliged to hold regulatory capital, the cost of which is the KVA and face the costs of funding uncollateralizedpositions through collateralized markets, known as FVA. Both [3] and [15] extend the approach of [21], in4hich unilateral collateralization was considered. We extend their approach to derive the value of ˆ u ( t, x )when the underlying follows the jump-diffusion defined in (1). We assume a one-dimensional underlyingdiffusion and consider all rates to be deterministic and, for ease of notation, constant. We specify differentrates, defined in Table 3.1, for different types of lending.Rate Definition Rate Definition r the risk-free rate r R the rate received on funding secured by theunderlying asset r D the dividend rate in case the stock paysdividends r F the rate received on unsecured funding r B the yield on a bond of the bank B r C the yield on the bond of the counter-party Cλ B λ B := r B − r λ C λ C := r C − rλ F λ F := r F − r R B the recovery rate of the bank R C the recovery rate of the counterpartyTable 3.1: Definitions of the rates used throughout the paper.Assume that the parties B and C enter into a derivative contract on the spot asset that pays the bank B the amount φ ( X T ) at maturity T . The value of this derivative to the bank at time t is denoted byˆ u ( t, x, J B , J C ) and depends on the value of the underlying X and the default states J B and J C of thebank B and counterparty C , respectively. Define I T C to be the initial margin posted by the bank to thecounterparty, I F C the initial margin posted by the counterparty to the bank and I V ( t ) to be the variationmargin on which a rate r I is paid or received. The initial margin is constant throughout the duration of thecontract. Let K ( t ) be the regulatory capital on which a rate of r K is paid/received.The cashflows are viewed from the perspective of the bank B . At the default time of either the counter-party or the bank, the value of the derivative to the bank ˆ u ( t, x ) is determined with a mark-to-market rule M , which may be equal to either the derivative value ˆ u ( t, x, ,
0) prior to default or the risk-free derivativevalue u ( t, x ), depending on the specifications in the ISDA master agreement. Denote by τ B and τ C the ran-dom default times of the bank and the counterparty respectively. We will use the notation x + = max( x, x − = min( x, I V + I F C . If the outstanding value M − ( I V + I F C ) is negative, the bank has to pay the fullamount ( M − I V − I F C ) − , while if the contract has a positive value to the bank, it will recover only R C ( M − I V − I F C ) + . Using a similar argument in case the bank defaults, we find the following boundaryconditions: θ Bt := ˆ u ( t, x, ,
0) = I V ( t ) − I T C + ( M − I V ( t ) + I T C ) + + R B ( M − I V ( t ) + I T C ) − ,θ Ct := ˆ u ( t, x, ,
1) = I V ( t ) + I F C + R C ( M − I V ( t ) − I F C ) + + ( M − I V ( t ) − I F C ) − , so that the portfolio value at default is given by θ τ = 1 τ C <τ B θ Cτ + 1 τ B <τ C θ Bτ , τ = min( τ B , τ C ). Further we introduce the default risky, zero-recovery, zero-coupon bonds (ZCBs) P B and P C with respective maturities T B and T C with face value one if the issuer has not defaulted, and zerootherwise. Assume the dynamics for P Bt and P Ct to be given by P Bt = { τ B >t } e r B t and P Ct = { τ C >t } e r C t ,so that dP Bt = r B P Bt dt − P Bt − d J Bt ,dP Ct = r C P Ct dt − P Ct − d J Ct , with J Bt = τ B ≤ t and J Ct = τ C ≤ t , where the default times τ B and τ C are defined in a canonical way asthe first arrival time of a doubly stochastic Poisson process with intensity functions γ B and γ C , respectively(see also the definition of the defaultable asset in (1)). We define the market interest rates for B and C tobe r B = r + γ B and r C = r + γ C , so that by the usual arguments (see, for instance, [16, Section 2.2]) thediscounted bonds e − rt P Bt and e − rt P Ct are martingales under the risk-neutral measure.We construct a hedging portfolio consisting of the shorted derivative, α C units of P C , α B units of P B and g units of cash: Π( t ) = − ˆ u ( t, x ) + α B ( t ) P Bt + α C ( t ) P Ct + g ( t ) . In other words, since we assume both the underlying asset process and the tradeable bonds P B and P C tobe risk-neutral, we focus on hedging the risk arising from the defaults of both B and C by means of thedefault-risky bonds.If the value of the derivative is positive to B , it will incur a cost at the counterparties’ default. To hedgethis, B shorts P C , i.e. α C ≤
0. If we assume B can borrow the bond close to the risk-free rate r (i.e. nohaircut) through a repurchase agreement, it will incur financing costs of rα C ( t ) P Ct dt . The cashflows fromthe collateralization follow from the rate r T C received and r F C paid on the initial margin and the rate r I paid or received on the collateral, depending on whether I V >
0, and the bank receives collateral, or I V < r K K ( t ).Finally, the rates r and r F are respectively received or paid on the surplus cash in the account. This cashconsists of the gap between the shorted derivative value and the collateral and the cost of buying α B bonds P B in order for B to hedge its own default, i.e. − ˆ u ( t, x ) − I V ( t ) + I T C − α B ( t ) P Bt . Thus, the total changein the cash account is given by dg ( t ) =[ − rα C ( t ) P Ct + r T C I T C − r F C I F C − r I I V ( t ) − r K K ( t )+ r ( − ˆ u ( t, x ) − I V ( t ) + I T C − α B ( t ) P Bt ) + λ F ( − ˆ u ( t, x ) − I V ( t ) + I T C − α B ( t ) P Bt ) − ] dt. Note that this is in contrast with the change in cash in a portfolio without the XVA arising from the differenttypes of funding, i.e. where we assume the cash in the portfolio simply earns the risk-free rate dg ( t ) = − r ˆ u ( t, x ) dt. Assuming the portfolio is self-financing we have d Π( t ) = − d ˆ u ( t, x ) + α B ( t ) dP Bt + α C ( t ) dP Ct + dg ( t ) . u ( t, x ) gives us: d ˆ u ( t, x ) = L ˆ u ( t, x ) dt + σ ( t, x ) ∂ x ˆ u ( t, x ) dW t + (cid:90) R (ˆ u ( t, x + q ) − ˆ u ( t, x )) d ˜ N ( t, x, dq ) − ( θ B − ˆ u ( t, x )) d J Bt − ( θ C − ˆ u ( t, x )) d J Ct , with the operator L as in (4). Thus, we find, d Π = − L ˆ u ( t, x ) dt − σ ( t, x ) ∂ x ˆ u ( t, x ) dW t − (cid:90) R (ˆ u ( t, x + q ) − ˆ u ( t, x )) d ˜ N ( t, x, dq )+ ( θ B − ˆ u ( t, x )) d J Bt + ( θ C − ˆ u ( t, x )) d J Ct − α B ( t ) P Bt − d J Bt − α C ( t ) P Ct − d J Ct + [ α B ( t ) λ B P Bt + α C ( t ) λ C P Ct + ( r T C + r ) I T C − r F C I F C − ( r I + r ) I V ( t ) − r K K ( t ) + r ˆ u ( t, x ) + λ F ( − ˆ u ( t, x ) − I V ( t ) + I T C − α B ( t ) P Bt ) − ] dt. By choosing α B = − θ B − ˆ u ( t, x ) P B , α C = − θ C − ˆ u ( t, x ) P C , we hedge the jump-to-default risk in the hedging portfolio, i.e., d Π = − L ˆ u ( t, x ) dt + σ ( t, x ) ∂ x ˆ u ( t, x ) dW t − (cid:90) R (ˆ u ( t, x + q ) − ˆ u ( t, x )) d ˜ N ( t, X t − , dq )+ [ − ( θ B − ˆ u ( t, x )) λ B − ( θ C − ˆ u ( t, x )) λ C + ( r T C + r ) I T C − r F C I F C − ( r I + r ) I V ( t ) − r K K ( t ) + r ˆ u ( t, x ) + λ F ( θ B − I V ( t ) + I T C ) − ] dt. Then, using the fact that the portfolio has to satisfy the martingale condition in the risk-neutral world, i.e. E [ d Π] = 0, we find the non-linear pricing PIDE to be L ˆ u ( t, x ) = f ( t, x, ˆ u ( t, x )) , (6)where we have defined f ( t, x, ˆ u ( t, x )) = − ( θ B ( t ) − ˆ u ( t, x )) λ B − ( θ C ( t ) − ˆ u ( t, x )) λ C + ( r T C + r ) I T C − r F C I F C − ( r I + r ) I V ( t ) − r K K ( t ) + r ˆ u ( t, x ) + λ F ( θ B − I V ( t ) + I T C ) − . In this section we will cast the PIDE in (6) in the form of a Backward Stochastic Differential Equation. Inthe methods where we make use of BSDEs we assume γ ( t, x ) = 0. We begin by recalling the non-linearFeynman-Kac theorem in the presence of jumps, see Theorem 4.2.1 in [7]. Theorem 1 (Non-linear Feynman-Kac Theorem) . Consider X t as in (1) . We assume µ , σ and a to beLipschitz continuous in x and additionally | a ( t, x ) | ≤ K . Consider the BSDE Y t = φ ( X T ) + (cid:90) Tt f (cid:18) s, X s , Y s , Z s , a ( s, X s − ) (cid:90) R V s ( q ) δ ( q ) ν ( dq ) (cid:19) ds − (cid:90) Tt Z s dW s (cid:90) Tt (cid:90) R V s ( q ) d ˜ N s ( s, X s , q ) , (7) where the generator f is continuous and satisfies the Lipschitz condition in the space variables, δ is a mea-surable, bounded function and the terminal condition φ ( x ) is measurable and Lipschitz continuous. Considerthe non-linear PIDE Lu ( t, x ) = f ( t, x, u ( t, x ) , ∂ x u ( t, x ) σ ( t, x ) , a ( t, x ) (cid:82) R ( u ( t, x + q ) − u ( t, x )) δ ( q ) ν ( dq )) ,u ( T, x ) = ψ ( x ) . (8) If the PIDE in (8) has a solution u ( t, x ) ∈ C , , the FBSDE in (7) has a unique solution ( Y t , Z t , V t ( q )) thatcan be represented as Y t,xs = u ( s, X t,xs ) ,Z t,xs = ∂ x u ( s, X t,xs ) σ ( s, X t,xs ) ,V t,xs ( q ) = u ( s, X t,xs + q ) − u ( s, X t,xs ) , q ∈ R , for all s ∈ [ t, T ] , where Y is a continuous, real-valued and adapted process and where the control processes Z and V are continuous, real-valued and predictable. In our case, the BSDE corresponding to the PIDE in (6) reads Y t = φ ( X T ) + (cid:90) Tt f ( s, X s , Y s ) ds − (cid:90) Tt Z s dW s − (cid:90) Tt (cid:90) R V s ( q ) d ˜ N ( s, X s , dq ) , (9)where we have defined the driver function to be f ( t, x, y ) = − λ B ( θ B − y ) − λ C ( θ C − y ) + ( r T C + r ) I T C − r F C I F C − ( r I + r ) I V ( t ) − r K K ( t ) + ry + λ F ( θ B − I V ( t ) + I T C ) − . Following [11], one can derive that the KVA is a function of trade properties (i.e. maturity, strike) and/orthe exposure at default, which in turn is a function of the portfolio value, so that the cost of holding thecapital can be rewritten as r K K ( t ) = r K c ˆ u ( t, x ) , with c being a function of the trade properties. Thecollateral is paid when the portfolio has a negative value, and received when the collateral has a positivevalue. Assuming the collateral is a multiple of the portfolio value we have I V ( t ) = c ˆ u ( t, x ), where c issome constant. Then, the driver function is simply a function of the portfolio value. Remark 2.
Note that in the case of ‘no collateralization’ or ‘perfect collateralization’, the driver functionreduces to f ( t, ˆ u ( t, x )) = r u ( t ) max(ˆ u ( t, x ) , r u here left unspecified. In this case the BSDEis similar to the one considered in [22]. 8 Solving FBSDEs
In this section we extend the BCOS method from [23] to solving FBSDEs under local L´evy models withvariable coefficients and jumps (without default, i.e. γ ( t, x ) = 0). The conditional expectations resulting fromthe discretization of the FBSDE are approximated using the COS method. This requires the characteristicfunction, which we approximate using the Adjoint Expansion Method of [19] and [2]. Consider the forward process X t as in (1) and the BSDE Y t as in (9) with a more general driver function f ( t, x, y, z ). Define a partition 0 = t < t < ... < t N = T of [0 , T ] with a fixed time step ∆ t = t n +1 − t n , for n = N − , ...
0. Rewriting the set of FBSDEs we find, X n +1 = X n + (cid:90) t n +1 t n µ ( s, X s ) ds + (cid:90) t n +1 t n σ ( s, X s ) dW s + (cid:90) t n +1 t n (cid:90) R qd ˜ N s ( s, X s − , dq ) ,Y n = Y n +1 + (cid:90) t n +1 t n f ( s, X s , Y s , Z s ) ds − (cid:90) t n +1 t n Z s dW s − (cid:90) t n +1 t n (cid:90) R V s ( q ) d ˜ N s ( s, X s − , dq ) . (10)One can obtain an approximation of the process Y t by taking conditional expectations with respect to theunderlying filtration G n , using the independence of W t and ˜ N t ( t, X t − , dq ) and by approximating the integralsthat appear with a theta-method, as first done in [24] and extended to BSDEs with jumps in [23]: Y n ≈ E n [ Y n +1 ] + ∆ tθ f ( t n , X n , Y n , Z n ) + ∆ t (1 − θ ) E n [ f ( t n +1 , X n +1 , Y n +1 , Z n +1 )] . Let ∆ W s := W s − W n for t n ≤ s ≤ t n +1 . Multiplying both sides of equation (10) by ∆ W n +1 , takingconditional expectations and applying the theta-method gives Z n ≈ − θ − (1 − θ ) E n [ Z n +1 ] + 1∆ t θ − E n [ Y n +1 ∆ W n +1 ]+ θ − (1 − θ ) E n [ f ( t n +1 , X n +1 , Y n +1 , Z n +1 ) ∆ W n +1 ] . Since in our scheme the terminal values are functions of time t and the Markov process X , it is easily seenthat there exist deterministic functions y ( t n , x ) and z ( t n , x ) so that Y n = y ( t n , X n ) , Z n = z ( t n , X n ) . The functions y ( t n , x ) and z ( t n , x ) are obtained in a backward manner using the following scheme y ( t N , x ) = φ ( x ) , z ( t N , x ) = ∂ x φ ( x ) σ ( t N , x ) , for n = N − , ..., y ( t n , x ) = E n [ y ( t n +1 , X n +1 )] + ∆ tθ f ( t n , x ) + ∆ t (1 − θ ) E n [ f ( t n +1 , X n +1 )] , (11) z ( t n , x ) = − − θ θ E n [ z ( t n +1 , X n +1 )] + 1∆ t θ − E n [ y ( t n +1 , X n +1 )∆ W n +1 ] (12)+ 1 − θ θ E n [ f ( t n +1 , X n +1 )∆ W n +1 ] , f ( t, X t ) := f ( t, X t , y ( t, X t ) , z ( t, X t )) . In the case θ > y ( t n , x ) in (11) and we use P Picard iterationsstarting with initial guess E n [ y ( t n +1 , X n +1 )] to determine y ( t n , x ). Is it well-known (see, for instance, [16, Section 2.2]) that the risk-free pre-default price u ( t, x ) of a Europeanoption on the defaultable asset S t with maturity T and payoff φ ( X T ) is given by u ( t, x ) = { ζ>t } e − r ( T − t ) E (cid:104) e − (cid:82) Tt γ ( s,X s ) ds φ ( X T ) | X t (cid:105) , t ≤ T, in the measure corresponding to the dynamics in (1). Thus, in order to compute the price of an option, wemust evaluate functions of the form v ( t, x ) := E (cid:104) e − (cid:82) Tt γ ( s,X s ) ds φ ( X T ) | X t = x (cid:105) . (13)Under standard assumptions, by the Feynman-Kac theorem, v can be expressed as the classical solution ofthe following Cauchy problem Lv ( t, x ) = 0 , t ∈ [0 , T [ , x ∈ R ,v ( T, x ) = φ ( x ) , x ∈ R , (14)with L as in (4).The function v in (13) can be represented as an integral with respect to the transition distribution of thedefaultable log-price process log S t : v ( t, x ) = (cid:90) R φ ( y )Γ( t, x ; T, dy ) , where Γ( t, x ; T, dy ) is the Green’s function of the PIDE in (14) and we say that its Fourier transformˆΓ( t, x ; T, ξ ) := F (Γ( t, x ; T, · ))( ξ ) := (cid:90) R e iξy Γ( t, x ; T, dy ) , ξ ∈ R , is the characteristic function of log S . Following [19] and [2] we expand the state-dependent coefficients s ( t, x ) := σ ( t, x )2 , µ ( t, x ) , γ ( t, x ) , a ( t, x ) , around some point ¯ x . The coefficients s ( t, x ), γ ( t, x ) and a ( t, x ) are assumed to be continuously differentiablewith respect to x up to order n ∈ N .Introduce the n th-order approximation of L in (4): L n = L + n (cid:88) k =1 (cid:16) ( x − ¯ x ) k µ k ( t ) + ( x − ¯ x ) k s k ( t ) ∂ xx − ( x − ¯ x ) k γ k ( t )10 (cid:90) R ( x − ¯ x ) k a k ( t ) ν ( dq )( e q∂ x − − q∂ x ) (cid:17) , where L = ∂ t + µ ( t ) ∂ x + s ( t ) ∂ xx − γ ( t ) + (cid:90) R a ( t ) ν ( dq )( e q∂ x − − q∂ x ) , and s k = ∂ kx s ( · , ¯ x ) k ! , γ k = ∂ kx γ ( · , ¯ x ) k ! , µ k ( dq ) = ∂ kx µ ( · , ¯ x ) k ! , a k = ∂ kx a ( · , ¯ x ) k ! k ≥ . The basepoint ¯ x is a constant parameter which can be chosen freely. In general the simplest choice is ¯ x = x (the value of the underlying at initial time t ).Assume for a moment that L has a fundamental solution G ( t, x ; T, y ) that is defined as the solution ofthe Cauchy problem L G ( t, x ; T, y ) = 0 t ∈ [0 , T [ , x ∈ R ,G ( T, · ; T, y ) = δ y . In this case we define the n th-order approximation of Γ asΓ ( n ) ( t, x ; T, y ) = n (cid:88) k =0 G k ( t, x ; T, y ) , where, for any k ≥ T, y ), G k ( · , · ; T, y ) is defined recursively through the following Cauchy problem L G k ( t, x ; T, y ) = − k (cid:80) h =1 ( L h − L h − ) G k − h ( t, x ; T, y ) t ∈ [0 , T [ , x ∈ R ,G k ( T, x ; T, y ) = 0 , x ∈ R . Notice that L k − L k − =( x − ¯ x ) k µ h ( t ) ∂ x + ( x − ¯ x ) k s k ( t ) ∂ xx − ( x − ¯ x ) k γ k ( t )+ (cid:90) R ( x − ¯ x ) k a k ( t ) ν ( dq )( e q∂ x − − q∂ x ) . Correspondingly, the n th-order approximation of the characteristic function ˆΓ is defined to beˆΓ ( n ) ( t, x ; T, ξ ) = n (cid:88) k =0 F (cid:0) G k ( t, x ; T, · ) (cid:1) ( ξ ) := n (cid:88) k =0 ˆ G k ( t, x ; T, ξ ) , ξ ∈ R . Now, by transforming the simplified Cauchy problems into adjoint problems and solving these in the Fourierspace we findˆ G ( t, x ; T, ξ ) = e iξx e (cid:82) Tt ψ ( s,ξ ) ds , ˆ G k ( t, x ; T, ξ ) = − (cid:90) Tt e (cid:82) Ts ψ ( τ,ξ ) dτ F (cid:32) k (cid:88) h =1 (cid:16) ˜ L ( s, · ) h ( s ) − ˜ L ( s, · ) h − ( s ) (cid:17) G k − h ( t, x ; s, · ) (cid:33) ( ξ ) ds, ψ ( t, ξ ) = iξµ ( t ) + s ( t ) ξ + (cid:90) R a ν ( t, dq )( e izξ − − izξ ) , ˜ L ( t,y ) h ( t ) − ˜ L ( t,y ) h − ( t ) = µ h ( t ) h ( y − ¯ x ) h − + µ h ( t )( y − ¯ x ) h ∂ y − γ h ( t )( y − ¯ x ) h + s h ( t ) h ( h − y − ¯ x ) h − + s h ( t )( y − ¯ x ) h − (2 h∂ y + ( y − ¯ x ) ∂ yy )+ (cid:90) R a h ( t )¯ ν ( dq ) (cid:0) ( y + q − ¯ x ) h e q∂ y − ( y − ¯ x ) h − q (cid:0) h ( y − ¯ x ) h − − ( y − ¯ x ) h ∂ y (cid:1)(cid:1) , where ¯ ν ( dq ) = ν ( − dq ). Remark 3.
After some algebraic manipulations it can be shown, see [2], that the characteristic functionapproximation of order n is a function of the formˆΓ ( n ) ( t, x ; T, ξ ) := e iξx n (cid:88) k =0 ( x − ¯ x ) k g n,k ( t, T, ξ ) , (15)where the coefficients g n,k , with 0 ≤ k ≤ n , depend only on t, T and ξ , but not on x . The approximationformula can thus always be split into a sum of products of functions depending only on ξ and functions thatare linear combinations of ( x − ¯ x ) m e iξx , m ∈ N . Remark 4 (Error estimates for the approximated characteristic function) . Similar to the derivation in [2],one can derive the error bounds for the characteristic function approximation. Let n = 0 , s ( t, x ), γ ( t, x ) and a ( t, x ) are continuously differentiable with bounded derivatives up to order n .For the n th-order approximation Γ ( n ) ( t, x ; T, ξ ), for any ¯ x ∈ R , (cid:12)(cid:12)(cid:12) Γ( t, x ; T, ξ ) − Γ ( n ) ( t, x ; T, ξ ) (cid:12)(cid:12)(cid:12) ≤ C ( T, ξ )(( T − t ) + ( T − t )( x − ¯ x )) n +12 . Note that if ¯ x = x , the bound reduces to C ( T, ξ )( T − t ) n +1 . The conditional expectations are approximated using the COS method, which was developed in [9] andapplied to FBSDEs with jumps in [23]. The conditional expectations arising in the equations (11)-(12)are all of the form E n [ h ( t n +1 , X n +1 )] or E n [ h ( t n +1 , X n +1 )∆ W n +1 ]. The COS formula for the first type ofconditional expectation reads E xn [ h ( t n +1 , X n +1 )] ≈ J − (cid:88) (cid:48) j =0 H j ( t n +1 )Re (cid:18) ˆΓ (cid:18) t n , x ; t n +1 , jπb − a (cid:19) exp (cid:18) ijπ − ab − a (cid:19)(cid:19) , where (cid:88) (cid:48) denotes an ordinary summation with the first term weighted by one-half, J > H j ( t n +1 ) denotes the j th Fourier-cosine coefficients of the function h ( t n +1 , x ) and ˆΓ ( t n , x ; t n +1 , ξ ) is the conditional characteristic function of the process X n +1 given X n = x .For the second type of conditional expectation, using integration by parts, we obtain E xn [ h ( t n +1 , X n +1 )∆ W n ] 12 ∆ tσ ( t n , x ) J − (cid:88) (cid:48) j =0 H j ( t n +1 )Re (cid:18) i jπb − a ˆΓ (cid:18) t n , x ; t n +1 , jπb − a (cid:19) exp (cid:18) ijπ − ab − a (cid:19)(cid:19) . See [23] for the full derivations.
Remark 5.
Note that these formulas are obtained by using an Euler approximation of the forward processand using the 2nd-order approximation of the characteristic function of the actual process. We have foundthis to be more exact than using the characteristic function of the Euler process, which is equivalent to usingjust the 0th-order approximation of the characteristic function.Finally we need to approximate the Fourier-cosine coefficients H j ( t n +1 ) of h ( t n +1 , x ) at time points t n ,where n = 0 , ..., N . The Fourier-cosine coefficient of h at time t n +1 is defined by H j ( t n +1 ) = 2 b − a (cid:90) ba h ( t n +1 , x ) cos (cid:18) jπ x − ab − a (cid:19) dx. Due to the structure of the approximated characteristic function of the local L´evy process, see (15), thecoefficients of the functions z ( t n +1 , x ) and the explicit part of y ( t n +1 , x ) can be computed using the FFTalgorithm, as we do in Appendix A, because of the matrix in (23) being of a certain form with constantdiagonals. In order to determine F j ( t n +1 ), the Fourier-Cosine coefficient of the function f ( t n +1 , x, y ( t n +1 , x ) , z ( t n +1 , x )) , due to the intricate dependence on the functions z and y we choose to approximate the integral in F j by adiscrete Fourier-Cosine transform (DCT). For the DCT we compute the integrand, and thus the functions z ( t n +1 , x ) and y ( t n +1 , x ), on an equidistant x -grid. Note that in this case we can easily approximate all Fourier-Cosine coefficients with a DCT (instead of the FFT). If we take J grid points defined by x i := a + ( i + ) b − aJ and ∆ x = b − aJ we find, using the mid-point integration rule, the approximation H j ( t n +1 ) ≈ J J − (cid:88) (cid:48) i =0 h ( t n +1 , x i ) cos (cid:18) jπ i + 12 J (cid:19) , which can be calculated using the DCT algorithm, with a computational complexity of O ( J log J ). Remark 6.
We define the truncation range [ a, b ] as follows:[ a, b ] := (cid:20) c − L (cid:113) c + √ c , c + L (cid:113) c + √ c (cid:21) , (16)where c n is the n th cumulant of log-price process log S , as proposed in [8]. The cumulants are calculatedusing the 0th-order approximation of the characteristic function. The method in Section 4 allows us to compute the XVA as in (5), consisting of CVA, DVA, MVA, KVAand FVA. In this section, we apply this method to computing Bermudan derivative values with XVA. The13esulting method – the solution of the non-linear XVA PDE through a BSDE-type method – is an efficientalternative to finite-difference methods as well as to the Monte-Carlo based method developed in [22]. Theefficiency is both due to the availability of the characteristic function in closed form through the adjointexpansion method and the fast convergence of the COS method. Furthermore, in finite difference methodscomplications may arise in the implementation of the scheme for jump diffusions. Since our proposed methodworks in the Fourier space, the jump component is easily handled by means of an additional term in thecharacteristic function and does not cause any further difficulties.For the CVA component in the XVA we develop an alternative method, which due to the ability of theFFT, results in a particularly efficient computation.
Consider an OTC derivative contract between the bank B and the counterparty C on the underlying asset S t given by (1) with γ ( t, x ) = 0 with a Bermudan-type exercise possibility: there is a finite set of so-calledexercise moments { t , ..., t M } prior to the maturity, with 0 ≤ t < t < · · · < t M = T . The payoff from thepoint-of-view of bank B is given by φ ( t m , X t m ). Denote ˆ u ( t, x ) to be the risky Bermudan option value and c ( t, x ) the continuation value. By the dynamic programming approach, the value for a Bermudan derivativewith XVA and M exercise dates t , ..., t M can be expressed by a backward recursion asˆ u ( t M , x ) = φ ( t M , x ) , and the continuation value solves the non-linear PIDE defined in (6) Lc ( t, x ) = f ( t, x, c ( t, x )) , t ∈ [ t m − , t m [ c ( t m , x ) = ˆ u ( t m , x )ˆ u ( t m − , x ) = max { Φ( t m − , x ) , c ( t m − , x ) } , m ∈ { , . . . , M } . The derivative value is set to be ˆ u ( t, x ) = c ( t, x ) for t ∈ ] t m − , t m [, and, if t >
0, also for t ∈ [0 , t [. Thepayoff function might take on various forms:1. (Portfolio) Following [22], we can consider X t to be the process of a portfolio which can take on bothpositive and negative values. Then, when exercised at time t m , bank B receives the portfolio so that φ ( t m , x ) = e x .2. (Bermudan option) In case the Bermudan contract is an option, the option value to the bank cannot have a negative value for the bank. At the same time, in case of default of the bank itself, thecounterparty loses nothing. In this case the framework simplifies to one with unilateral collateralizationand default risk and the payoff at time t m , if exercised, is given by φ ( t m , x ) = ( K − e x ) + for a put and φ ( t m , x ) = ( e x − K ) + for a call with K being the strike price.3. (Bermudan swaptions) A Bermudan swaption is an option in which the holder, bank B , has the rightto exercise and enter into an underlying swap with fixed end date t M +1 . If the swaption is exercised14t time t m the underlying swap starts with payment dates T m = { t m +1 , ..., t M +1 } . Working under theforward measure corresponding to the last reset date t M , the payoff function is given by φ ( t m , x ) = N S (cid:32) M (cid:88) k = m P ( t m , t k +1 , x ) P ( t m , t M ) ∆ t (cid:33) max( c p ( S ( t m , T m , x ) − K ) , , where N S is the notional, c p = 1 for a payer swaption and c p = − P ( t m , t k , x )is the price of a ZCB conditional on X t m = x and S ( t m , T m , x ) is the forward swap rate given by S ( t m , T m , x ) = (cid:18) − P ( t m , t m +1 , x ) P ( t m , t M , x ) (cid:19) (cid:14) (cid:32) M (cid:88) k = m P ( t m , t k +1 , x ) P ( t m , t M , x ) ∆ t (cid:33) . To solve for the continuation value we define a partition with N steps t m − = t ,m < t ,m < t ,m < ... 1. In order to compute this at each time step t n,m we thus need to evaluate c ( t n,m , x ) on the x -grid with J equidistant points using formula (17). The matrix-vector product in theformula results in a computational time of order O ( J ). Remark 7 (Convergence of the Picard iterations) . A Picard iteration is used to find the fixed-point c of c = ∆ tθ f ( t n,m , x, c ) + h ( t n,m , x ) , where f ( t, x, c ) and h ( t, x ) are respectively the implicit and explicit partsof the equation. Due to the computational domain of c ( t, x ) being bounded by [ a, b ], we can thus say that f ( t, x, c ( t, x )) is also bounded. If the driver function f ( t, x, c ) is Lipschitz continuous in c , i.e. ∃ L Lipz suchthat | f ( t, x, c ) − f ( t, x, c ) | ≤ L Lipz | c − c | , and ∆ t n is small enough such that ∆ tθ L Lipz < 1, a unique15xed-point exists and the Picard iterations converge towards that point for any initial guess. In particular,for the XVA case the non-linearity is of the form f ( t, x, c ) = − r max( c, L Lipz = 1. Thus for ∆ t sufficiently small, the Picard iteration converges to a unique fixed-point.The total algorithm for computing the value of a Bermudan contract with XVA can be summarised asin Algorithm 1 in Figure 5.1. The total computational time for the algorithm is of order O ( M · N ( J + J + P J + J log J )) , (18)consisting of the computation for M · N times the computation of the characteristic function on the x -grid (due to the availability of the analytical approximation) of O ( J ), computation of the matrix-vectormultiplications in the formulas for c ( t n,m , x ) and z ( t n,m , x ) of O ( J ), initialization of the Picard methodwith E n [ c ( t n +1 , X n +1 ] in O ( J ) operations, computation of the P Picard approximations for c ( t n,m , x ) in O ( P J ) and computing the Fourier coefficients F j ( t n ) and C j ( t n ) with the DCT in O ( J log J ) operations.1. Define the x -grid with J grid points given by x i = a + ( i + ) b − aJ for i = 0 , ..., J − c ( t N,M , x ) = ˆ u ( t M , x ) on the x -grid and compute theterminal coefficients C j ( t M ) and F j ( t M ) using the DCT.3. Recursively for the exercise dates m = M − , ..., n = N − , ..., c ( t n,m , x ) using formula (17) and use this to determine f ( t n,m , x, c ( t n,m , x )) onthe x -grid.ii. Subsequently, use these to determine F j ( t n,m ) and C j ( t n,m ) using the DCT.(b) Compute the new terminal condition c ( t N,m − , x ) = max { φ ( t ,m , x ) , c ( t ,m , x ) } (either ana-lytically or numerically) and the corresponding Fourier-cosine coefficient.4. Finally ˆ u ( t , x ) = c ( t , , x ).Figure 5.1: Algorithm 1: Bermudan derivative valuation with XVA In this section we present an efficient alternative way of calculating the CVA term in (5) in the case ofunilateral CCR using a Fourier-based method. Due to the ability of using the FFT this method is considerablyfaster for computing the CVA than the method presented in Section 5.1. We use the definition of CVA attime t given by CVA( t ) = ˆ u ( t, X t ) − u ( t, X t ) , u ( t, X t ) is as usual the default-free value of the Bermudan option ( γ ( t, x ) = 0), while ˆ u ( t, X t ) is thevalue including default ( γ ( t, x ) (cid:54) = 0). We consider the model as defined in (1). We will compute u ( t, X t ) andˆ u ( t, X t ) using the COS method and the approximation of the characteristic function (as derived in Section4.3), without default and with default, respectively. In case of a default the payoff becomes zero. Notethat the risky option value ˆ u ( t, x ) computed with the characteristic function for a defaultable underlyingcorresponds exactly to the option value in which the counterparty might default, with the probablity ofdefault, P D ( t ), defined as in (2). Thus, in this case we have unilateral CCR and ζ = τ C , the default time ofthe counterparty.Using the definition of the defaultable S t , it is well-known (see, for instance, [16, Section 2.2]) that therisky no-arbitrage value of the Bermudan option on the defaultable asset S t at time t is given byˆ u ( t, X t ) = { ζ>t } sup τ ∈{ t ,...,t M } E (cid:104) e − (cid:82) τt ( r + γ ( s,X s )) ds φ ( τ, X τ ) | X t (cid:105) . Remark 8 (Wrong-way risk) . By allowing the dependence of the default intensity on the underlying, asimplified form of wrong-way risk is already incorporated into the CVA valuation.For a Bermudan put option with strike price K , we simply have φ ( t, x ) = ( K − x ) + . By the dynamicprogramming approach, the option value can be expressed by a backward recursion asˆ u ( t M , x ) = { ζ>t M } max( φ ( t M , x ) , , and c ( t, x ) = E (cid:104) e (cid:82) tmt ( r + γ ( s,X s )) ds ˆ u ( t m , X t m ) | X t = x (cid:105) , t ∈ [ t m − , t m [ˆ u ( t m − , x ) = { ζ>t m − } max { φ ( t m − , x ) , c ( t m − , x ) } , m ∈ { , . . . , M } . (19)Thus to find the risky option price ˆ u ( t, X t ) one uses the defaultable asset with γ ( t, x ) representing the defaultintensity of the counterparty and in order to get the default-free value u ( t, X t ) one uses the default-free assetby setting γ ( t, x ) = 0. The CVA adjustment is calculated as the difference between the two. Both ˆ u ( t, x )and u ( t, x ) are calculated using the approximated characteristic function and the COS method applied tothe continuation value [2]. Due to the characteristic function being of the form (15), we are able to use theFFT in the matrix-vector multiplication when computing the continuation values of the Bermudan optionwith and without default, reducing this operation from O ( J ) to O ( J log J ). For more details, we referto Appendix A. The total complexity of the calculation of the CVA value for a Bermudan option with M exercise dates is then O ( M J log J ). Comparing this to (18), in which the most time-consuming operationswere indeed the matrix-vector products of order O ( J ) that resulted from the computation of the functionson the x -grid of size J , we conclude that the method for CVA computation is indeed significantly faster dueto the ability of using the FFT. In practice CVA is hedged and thus practitioners require efficient ways to compute the sensitivity of theCVA with respect to the underlying. The widely used bump- and revalue- method, while resulting in precise17alculations, might be slow to compute. Using the Fourier-based approach we find explicit formulas allowingfor an easy computation of the first- and second-order derivatives of the CVA with respect to the underlying.For the first-order and second-order Greeks we have∆ = e − r ( t − t ) J − (cid:88) (cid:48) j =0 Re (cid:18) e ijπ x − ab − a (cid:18) ijπb − a g dn, (cid:18) t , t , jπb − a (cid:19) + g dn, (cid:18) t , t , jπb − a (cid:19)(cid:19)(cid:19) V dj ( t ) − e − r ( t − t ) J − (cid:88) (cid:48) j =0 Re (cid:18) e ijπ x − ab − a (cid:18) ijπb − a g rn, (cid:18) t , t , jπb − a (cid:19) + g rn, (cid:18) t , t , jπb − a (cid:19)(cid:19)(cid:19) V rj ( t ) ,∂ ∆ ∂X = e − r ( t − t ) J − (cid:88) (cid:48) j =0 Re (cid:18) e ijπ x − ab − a (cid:18) − ijπb − a g dn, (cid:18) t , t , jπb − a (cid:19) − g dn, (cid:18) t , t , jπb − a (cid:19) + 2 ijπb − a g dn, (cid:18) t , t , jπb − a (cid:19) + (cid:18) ijπb − a (cid:19) g dn, (cid:18) t , t , jπb − a (cid:19) + 2 g dn, (cid:18) t , t , jπb − a (cid:19) (cid:19)(cid:19) V dj ( t ) − e − r ( t − t ) J − (cid:88) (cid:48) j =0 Re (cid:18) e ijπ x − ab − a (cid:18) − ijπb − a g rn, (cid:18) t , t , jπb − a (cid:19) − g rn, (cid:18) t , t , jπb − a (cid:19) − ijπb − a g rn, (cid:18) t , t , jπb − a (cid:19) + (cid:18) ijπb − a (cid:19) g rn, (cid:18) t , t , jπb − a (cid:19) + 2 g rn, (cid:18) t , t , jπb − a (cid:19) (cid:19)(cid:19) V j ( t ) r , where V dk and V rk are the Fourier-cosine coefficients with the defaultable and default-free characteristicfunction terms, g dn,h and g rn,h , respectively. In this section we present numerical examples to justify the accuracy of the methods in practice. We computethe XVA with the method presented in Section 5.1 and the CVA in the case of unilateral CCR with themethod from Section 5.2, which we show is more efficient for cases in which one only needs to compute theCVA. We compare the results of solving the BSDE with the COS method and the adjoint expansion of thecharacteristic function to the values obtained by using a least-squares Monte-Carlo method for computingthe conditional expected values in the BSDE as done in e.g. [1].The computer used in the experiments has an Intel Core i7 CPU with a 2.2 GHz processor. We use thesecond-order approximation of the characteristic function. We have found this to be sufficiently accurate bynumerical experiments and theoretical error estimates. The formulas for the second-order approximation aresimple, making the methods easy to implement. Here, we check the accuracy of the method from Section 5.1. We will compute the Bermudan option valuewith XVA using a simplified driver function given by f ( t, ˆ u ( t, x )) = − r max(ˆ u ( t, x ) , X t to be a portfolio process and the payoff, ifexercised at time t m , to be given by Φ( t m , x ) = x . In this case the value we can receive at every exercise date18s the value of the portfolio. Consider the model in Section 2 without default, with a local jump measureand a local volatility function with CEV-like dynamics and Gaussian jumps defined by σ ( x ) = be βx , (20) ν ( x, dq ) = λe βx √ πδ exp (cid:18) − ( q − m ) δ (cid:19) dq. (21)We assume the following parameters in equations (20)-(21), unless otherwise mentioned: b = 0 . β = − λ = 0 . δ = 0 . m = − . r = 0 . K = 1 and X = 0 (so that S = 1). In the LSM the number of timesteps is taken to be 100 and we simulate 10 paths. In the COS method we take J = 256, θ = 0 . N = 10, M = 10, making the total number of time steps N · M = 100. The truncation range is determinedas in (16) with L = 10. Due to the state-dependent coefficients in the underlying dynamics in (20)-(21) weuse the approximated characteristic function as derived in Section 4.2 with the second-order approximation,i.e. ˆΓ (2) ( t, x ; T, ξ ) and take ¯ x = x , where x = { x i } J − i =0 . Note that we thus compute the values, includingthose of the characteristic function, on the complete x -grid. In the final iteration when computing ˆ u ( t , X )we use ¯ x = X .In Table 6.1 we analyse the error in the approximation of ˆ u ( t , X ) with S = 0 . N and the number of grid points (and Fourier-cosine coefficients) J . We comparethe approximated COS value to the 95% confidence interval obtained by a LSM. Accurate results are quicklyobtained for small values of both J and N . In Figure 6.1 we plot the upper bound of the 95% confidenceinterval of the absolute error in the approximation for varying J and N . We observe approximately a linearconvergence and note that the error stops decreasing at some point for increasing values of J and N . Thiscan be due to the error being dominated by the approximated characteristic function. In particular weobserve that J = 32 and N = 10 seem to be sufficient parameters to achieve a satisfactory accuracy in theapproximation.The results for ˆ u ( t , X ) of the COS approximation method compared to a 95% confidence interval ofthe value obtained through a LSM are presented in Table 6.1. These results show that our method is ableto solve non-linear PIDEs accurately. The CPU time of the approximating method depends on the numberof time steps M · N and is approximately 5 · ( N · M ) ms. N = 1 N = 10 N = 20 N = 30 J = 8 6.4E-03 − − − − J = 16 2.3E-03 − − − − J = 32 1.7E-03 − − − − J = 64 1.4E-03 − − − − J = 128 1.7E-04 − − − − J = 256 2.1E-04 − − − − u (0 , X ) with S = 0 . J and N .19igure 6.1: Convergence of the upper bound of the 95% confidence interval of the absolute error in the COSapproximation ˆ u (0 , X ) with S = 0 . J and N .maturity T S MC value with XVA COS value with XVA0.5 0 0.03770 − − − − − − − − − − − − T = 0 . , 1) in the CEV-like modelfor the 2nd-order approximation of the characteristic function, and an LSM comparison. In this section we validate the accuracy of the method presented in Section 5.2 and compute the CVA inthe case of unilateral CCR under the model dynamics given in Section 2 with a local jump measure and alocal volatility function with CEV-like dynamics, Gaussian jumps defined by defined as in (21) and a localdefault function γ ( x ) = ce βx . We assume the same parameters as in Section 6.2, except r = 0 . 05 and wetake c = 0 . paths. In the COS method we take L = 10 and J = 100. Again, due to the state-dependent coefficients20n the underlying dynamics we use the approximated characteristic function as derived in Section 4.2 withthe second-order approximation, i.e. ˆΓ (2) ( t, x ; T, ξ ) and take ¯ x = X .The results for the CVA valuation with the FFT-based method and with LSM are presented in Table6.2. The CPU time of the LSM is at least 5 times the CPU time of the approximating method, which for M exercise dates is approximately 3 · M ms, thus more efficient than the computation of the XVA with themethod in Section 5.1. The optimal exercise boundary in Figure 6.2 shows that the exercise region becomeslarger when the probability of default increases; this is to be expected: in case of the default probabilitybeing greater, the option of exercising early is more valuable and used more often.maturity T strike K MC CVA COS CVA0.5 0.6 4 . · − − . · − . · − − · − − − − − − − − − − − T = 0 . , 1) in the CEV-like model forthe 2nd-order approximation of the characteristic function, and an LSM comparison.Figure 6.2: Optimal exercise boundary for a Bermudan put option (10 exercise dates, expiry T = 1) in theCEV-like model with varying default c = 0 , . , . 2. 21 Conclusion In this paper we considered pricing Bermudan derivatives under the presence of XVA, consisting of CVA,DVA, MVA, FVA and KVA. We derived the replicating portfolio with cashflows corresponding to the differentrates for different types of lending. This resulted in the PIDE in (6) and its corresponding BSDE (9). Wepropose to solve the BSDE using a Fourier-cosine method for the resulting conditional expectations andan adjoint expansion method for determining an approximation of the characteristic function of the localL´evy model in (1). This approach is extended to Bermudan option pricing in Section 5.1. In Section 5.2we presented an alternative for computing the CVA term in the case of unilateral collateralization (as isthe case when the derivative is an option) without the use of BSDEs. This results in an even more efficientmethod due to the ability to use the FFT. We verify the accuracy of both methods in Sections 6.1 and6.2 by comparing it to a LSM and conclude that the method from Section 5.1 is able to achieve a rapidconvergence and gives, already for small values of the discretization parameters an accurate result. Thealternative method for CVA computation from Section 5.2 is indeed more efficient than the BSDE methodfor computing just the CVA term. Acknowledgments We thank two anonymous referees for the comments and suggestions that have improved the quality of thispaper. This research is supported by the European Union in the the context of the H2020 EU Marie CurieInitial Training Network project named WAKEUPCALL. A The COS formulae Let, as usual, J denote the number of Fourier-cosine coefficients. Remembering that the expected value c ( t, x ) in (19) can be rewritten in integral form, we have c ( t, x ) = e − r ( t m − t ) (cid:90) R v ( t m , y )Γ( t, x ; t m , dy ) , t ∈ [ t m − , t m [ , where, v ( t m , y ) can be either u ( t m , y ) or ˆ u ( t m , y ). Then we use the Fourier-cosine expansion to get theapproximation:ˆ c ( t, x ) = e − r ( t m − t ) J − (cid:88) (cid:48) j =0 Re (cid:18) e − ijπ ab − a ˆΓ (cid:18) t, x ; t m , jπb − a (cid:19)(cid:19) V j ( t m ) , t ∈ [ t m − , t m [ (22) V j ( t m ) = 2 b − a (cid:90) ba cos (cid:18) jπ y − ab − a (cid:19) max { φ ( t m , y ) , c ( t m , y ) } dy, with φ ( t, x ) = ( K − e x ) + .We can recover the coefficients ( V j ( t m )) j =0 , ,...,J − from ( V j ( t m +1 )) j =0 , ,...,J − . To this end, we split theintegral in the definition of V j ( t m ) into two parts using the early-exercise point x ∗ m , which is the point where22he continuation value is equal to the payoff, i.e. c ( t m , x ∗ m ) = φ ( t m , x ∗ m ); this point can easily be found byusing the Newton method. Thus, we have V j ( t m ) = F j ( t m , x ∗ m ) + C j ( t m , x ∗ m ) , m = M − , M − , ..., , where F j ( t m , x ∗ m ) := 2 b − a (cid:90) x ∗ m a φ ( t m , y ) cos (cid:18) jπ y − ab − a (cid:19) dy,C j ( t m , x ∗ m ) := 2 b − a (cid:90) bx ∗ m c ( t m , y ) cos (cid:18) jπ y − ab − a (cid:19) dy, and V j ( t M ) = F j ( t M , log K ) . The coefficients F j ( t m , x ∗ m ) can be computed analytically using x ∗ m ≤ log K , and by inserting the ap-proximation (22) for the continuation value into the formula for C j ( t m , x ∗ m ) have the following coefficientsˆ C j for m = M − , M − , ..., C j ( t m , x ∗ m ) = 2 e − r ( t m +1 − t m ) b − a · J − (cid:88) (cid:48) k =0 V k ( t m +1 ) (cid:90) bx ∗ m Re (cid:18) e − ikπ ab − a ˆΓ (cid:18) t m , x ; t m +1 , kπb − a (cid:19)(cid:19) cos (cid:18) jπ x − ab − a (cid:19) dx. From (15) we know that the n th-order approximation of the characteristic function is of the form:ˆΓ ( n ) ( t m , x ; t m +1 , ξ ) = e iξx n (cid:88) h =0 ( x − ¯ x ) h g n,h ( t m , t m +1 , ξ ) , where the coefficients g n,h ( t, T, ξ ), with 0 ≤ k ≤ n , depend only on t, T and ξ , but not on x . Remark 9 (The defaultable and default-free characteristic functions) . To find u ( t, x ) we useˆΓ r ( t m , x ; t m +1 , ξ ) := e iξx n (cid:88) h =0 ( x − ¯ x ) h g rn,h ( t m , t m +1 , ξ ) , the characteristic function with γ ( t, x ) = 0. For ˆ u ( t, x ) we useˆΓ d ( t m , x ; t m +1 , ξ ) := e iξx n (cid:88) h =0 ( x − ¯ x ) h g dn,h ( t m , t m +1 , ξ ) , where γ ( t, x ) is chosen to be some specified function.Using (15) we can write the Fourier coefficients of the continuation value in vectorized form as: ˆC ( t m , x ∗ m ) = n (cid:88) h =0 e − r ( t m +1 − t m ) Re (cid:0) V ( t m +1 ) M h ( x ∗ m , b )Λ h (cid:1) , where V ( t m +1 ) is the vector [ V ( t m +1 ) , ..., V J − ( t m +1 )] T and M h ( x ∗ m , b )Λ h is a matrix-matrix product with M h a matrix with elements { M hk,j } J − k,j =0 defined as M hk,j ( x ∗ m , b ) := 2 b − a (cid:90) bx ∗ m e ijπ x − ab − a ( x − ¯ x ) h cos (cid:18) kπ x − ab − a (cid:19) dx, (23)23nd Λ h is a diagonal matrix with elements g n,h (cid:16) t m , t m +1 , jπb − a (cid:17) , j = 0 , . . . , J − . One can show, see [2], that the resulting matrix M h is a sum of a Hankel and Toeplitz matrix and thus theresulting matrix vector product can be calculated using a FFT. References [1] C. Bender and J. Steiner , Least-Squares Monte Carlo for Backward SDEs , in Numerical methodsin finance, Springer, 2012, pp. 257–289.[2] A. Borovykh, A. Pascucci, and C. W. Oosterlee , Pricing Bermudan options under local l´evymodels with default , Journal of Mathematical Analysis and Applications, 450 (2017), pp. 929–953.[3] C. Burgard and M. Kjaer , Partial differential equation representations of derivatives with bilateralcounterparty risk and funding costs , Journal of Credit Risk, 7 (2011), pp. 75–93.[4] A. Capponi, S. Pagliarani, and T. Vargiolu , Pricing vulnerable claims in a L´evy-driven model ,Finance Stoch., 18 (2014), pp. 755–789.[5] P. Carr and V. Linetsky , A jump to default extended CEV model: an application of Bessel processes ,Finance Stoch., 10 (2006), pp. 303–330.[6] C. de Graaf, Q. Feng, D. Kandhai, and C. Oosterlee , Efficient computation of exposure profilesfor counterparty credit risk , International Journal of Theoretical and Applied Finance, 4 (2014).[7] (cid:32)L. Delong , Backward stochastic differential equations with jumps and their actuarial and financialapplications: BSDEs with jumps , Springer Science & Business Media, 2013.[8] F. Fang and C. W. Oosterlee , A novel pricing method for European options based on Fourier-cosineseries expansions , SIAM J. Sci. Comput., 31 (2008/09), pp. 826–848.[9] , Pricing early-exercise and discrete Barrier options by Fourier-cosine series expansions , Numer.Math., 114 (2009), pp. 27–62.[10] Q. Feng and C. Oosterlee , Monte Carlo calculation of exposure profiles and Greeks for Bermudanand Barrier options under the Heston Hull-White model , submitted, (2014).[11] A. Green, C. Kenyon, and C. Dennis , KVA: Capital Valuation Adjustment , Risk, 12 (2014).[12] J. Hull and A. White , The FVA debate , Risk, 7 (2012), pp. 83–85.[13] , XVAs: A gap between theory and practice , Risk, (2016).[14] S. Jain and C. W. Oosterlee , The stochastic grid bundling method: Efficient pricing of Bermudanoptions and their Greeks , Appl. Math. Comput., 269 (2015), pp. 412–432.2415] A. Lesniewski and A. Richter , Managing counterparty credit risk via BSDEs , submitted, (2016).[16] V. Linetsky , Pricing equity derivatives subject to bankruptcy , Math. Finance, 16 (2006), pp. 255–282.[17] F. Longstaff and E. Schwartz , Valuing American options by simulation: A simple least-squaresapproach , Rev. Financ. Stud., 14 (2001), pp. 113–147.[18] M. Lorig, S. Pagliarani, and A. Pascucci , A family of density expansions for L´evy-type processes ,Ann. Appl. Probab., 25 (2015), pp. 235–267.[19] S. Pagliarani, A. Pascucci, and C. Riga , Adjoint expansions in local L´evy models , SIAM J.Financial Math., 4 (2013), pp. 265–296.[20] A. Pascucci , PDE and martingale methods in option pricing , vol. 2 of Bocconi & Springer Series,Springer, Milan; Bocconi University Press, Milan, 2011.[21] V. Piterbarg , Funding beyond discounting: Collateral agreements and derivatives pricing , Risk, 2(2010), pp. 97 –102.[22] , A non-linear PDE for XVA by forward Monte Carlo , Risk, 10 (2015).[23] M. J. Ruijter and C. W. Oosterlee , A Fourier-cosine method for an efficient computation ofsolutions to BSDEs , SIAM J. Sci. Comput., 37 (2015), pp. 859–889.[24] W. Zhao, Y. Li, and G. Zhang , A generalized θ -scheme for solving backward stochastic differentialequations.-scheme for solving backward stochastic differentialequations.