Efficient Computation of Various Valuation Adjustments Under Local L\'evy Models

TL;DR

This paper introduces a Fourier-based numerical method to efficiently price Bermudan derivatives with XVAs under local Le9vy models, accommodating local volatility and jump measures despite the lack of characteristic functions.

Contribution

It develops an innovative Fourier-based approach for solving FBSDEs in local Le9vy models, enabling accurate valuation of derivatives with XVAs.

Findings

Effective pricing of Bermudan options with XVAs under complex models

Approximation method compensates for unavailable characteristic functions

Framework accommodates local volatility and jump dynamics

Abstract

Various valuation adjustments, or XVAs, can be written in terms of non-linear PIDEs equivalent to FBSDEs. In this paper we develop a Fourier-based method for solving FBSDEs in order to efficiently and accurately price Bermudan derivatives, including options and swaptions, with XVA under the flexible dynamics of a local L\'evy model: this framework includes a local volatility function and a local jump measure. Due to the unavailability of the characteristic function for such processes, we use an asymptotic approximation based on the adjoint formulation of the problem.

Tables

Table 1. Table 3.1: Definitions of the rates used throughout the paper.

Rate	Definition	Rate	Definition
$r$	the risk-free rate	$r_R$	the rate received on funding secured by the underlying asset
$r_D$	the dividend rate in case the stock pays dividends	$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 counterparty $C$
${\lambda }_B$	${\lambda }_B:=r_B-r$	${\lambda }_C$	${\lambda }_C:=r_C-r$
${\lambda }_F$	${\lambda }_F:=r_F-r$	$R_B$	the recovery rate of the bank
$R_C$	the recovery rate of the counterparty

Table 2. Table 6.1: The 95% confidence interval of the absolute error in the COS approximation of u ^ ​ ( 0 , X 0 ) ^ 𝑢 0 subscript 𝑋 0 \hat{u}(0,X_{0}) with S 0 = 0.4 subscript 𝑆 0 0.4 S_{0}=0.4 compared to a LSM for varying parameters J 𝐽 J and N 𝑁 N .

	$N=1$	$N=10$	$N=20$	$N=30$
$J=8$	6.4E-03$-$6.9E-03	4.3E-03$-$4.8E-03	4.9E-03$-$5.3E-03	5.3E-03$-$5.8E-03
$J=16$	2.3E-03$-$2.7E-03	8.8E-04$-$1.3E-03	6.2E-04$-$1.1E-03	5.4E-04$-$9.2E-04
$J=32$	1.7E-03$-$2.0E-03	4.2E-04$-$8.3E-04	2.4E-04$-$6.3E-04	1.6E04$-$5.8E-04
$J=64$	1.4E-03$-$1.9E-03	2.2E-04$-$6.5E-04	1.6E-04$-$2.3E-04	1.2E-04$-$2.9E-04
$J=128$	1.7E-04$-$6.0E-04	2.1E-04$-$6.6E-04	2.3E-04$-$6.5E-04	1.9E-04$-$6.1E-04
$J=256$	2.1E-04$-$6.6E-04	3.7E-04$-$7.7E-04	1.5E-04$-$5.7E-04	1.2E-04$-$3.1E-04

Table 3. Table 6.2: A Bermudan put option with XVA (10 exercise dates, expiry T = 0.5 , 1 𝑇 0.5 1 T=0.5,1 ) in the CEV-like model for the 2nd-order approximation of the characteristic function, and an LSM comparison.

maturity $T$	$S_0$	MC value with XVA	COS value with XVA
0.5	0	0.03770$-$0.03838	0.03809
	0.2	0.2326$-$0.2330	0.2320
	0.4	0.4251$-$0.4254	0.4243
	0.6	0.6169$-$0.6171	0.6158
	0.8	0.8077$-$0.8079	0.8069
	1	1.000$-$1.000	1.0000
1	0	0.07374$-$0.07453	0.07228
	0.2	0.2611$-$0.2617	0.2606
	0.4	0.4461$-$0.4465	0.4454
	0.6	0.6288$-$0.6291	0.6288
	0.8	0.8126$-$0.8129	0.8113
	1	1.001$-$1.001	1.000

Table 4. Table 6.3: CVA for a Bermudan put option (10 exercise dates, expiry T = 0.5 , 1 𝑇 0.5 1 T=0.5,1 ) in the CEV-like model for the 2nd-order approximation of the characteristic function, and an LSM comparison.

maturity $T$	strike $K$	MC CVA	COS CVA
0.5	0.6	$4.200\cdot {10}^{-4}-4.807\cdot {10}^{-4}$	$1.113\cdot {10}^{-4}$
	0.8	0.001525$-$0.001609	9.869$\cdot {10}^{-4}$
	1	0.01254$-$0.01273	0.01138
	1.2	0.005908$-$0.005931	0.005937
	1.4	0.006657$-$0.06758	0.006898
	1.6	0.007795$-$0.008008	0.007883
1	0.6	8.673E-04$-$9.574E-04	4.463E-04
	0.8	0.005817$-$0.006040	0.003535
	1	0.02023$-$0.02054	0.01882
	1.2	0.01221$-$0.01222	0.1272
	1.4	0.01378$-$0.01391	0.01360
	1.6	0.01532$-$0.01502	0.01554

Equations

$$S_t$$

$$d{X}_t$$

$$d{\tilde{N}}_t(t,X_{t-},dq)$$

$$\zeta$$

$$\nu (t,X_{t-},dq)=a(t,X_{t-})\nu (dq).$$

$$\text{PD}(t):=\mathbb{P}(\zeta \le t)=1-\mathbb{E}\left[e^{-{\int }_0^t\gamma (s,X_s)ds}\right].$$

$${\int }_{\mathbb{R}}{e}^{|q|}a(t,x)\nu (dq)<\infty .$$

$$\mu (t,x)=\gamma (t,x)+r-\frac{{\sigma }^2(t,x)}{2}-a(t,x){\int }_{\mathbb{R}}\nu (dq)(e^q-1-q),$$

$$Lu(t,x)=$$

$$+a(t,x){\int }_{\mathbb{R}}\nu (dq)(u(t,x+q)-u(t,x)-q{\partial }_{x}u(t,x)).$$

$$\text{TVA}:=\hat{u}(t,x)-u(t,x)=\text{CVA}+\text{DVA}+\text{KVA}+\text{MVA}+\text{FVA}.$$

$$Lu(t,x)$$

$$u(T,x)$$

$${\theta }_t^B$$

$${\theta }_t^C$$

$${\theta }_{\tau }=1_{{\tau }^C<{\tau }^B}{\theta }_{\tau }^C+1_{{\tau }^B<{\tau }^C}{\theta }_{\tau }^B,$$

$$d{P}_t^B$$

$$d{P}_t^C$$

$$\Pi (t)=-\hat{u}(t,x)+{\alpha }_B(t)P_t^B+{\alpha }_C(t)P_t^C+g(t).$$

$$dg(t)=$$

$$+r(-\hat{u}(t,x)-I^V(t)+I_{TC}-{\alpha }_B(t)P_t^B)+{\lambda }_F(-\hat{u}(t,x)-I^V(t)+I_{TC}-{\alpha }_B(t)P_t^B{)}^{-}]dt.$$

$$dg(t)=-r\hat{u}(t,x)dt.$$

$$d\Pi (t)=$$

$$d\hat{u}(t,x)=$$

$$-({\theta }^B-\hat{u}(t,x))d{\mathcal{J}}_t^B-({\theta }^C-\hat{u}(t,x))d{\mathcal{J}}_t^C,$$

$$d\Pi =$$

$$+({\theta }^B-\hat{u}(t,x))d{\mathcal{J}}_t^B+({\theta }^C-\hat{u}(t,x))d{\mathcal{J}}_t^C-{\alpha }^B(t)P_{t-}^{B}d{\mathcal{J}}_t^B-{\alpha }^C(t)P_{t-}^{C}d{\mathcal{J}}_t^C$$

$$+[{\alpha }^B(t){\lambda }_B{P}_t^B+{\alpha }^C(t){\lambda }_C{P}_t^C+(r_{TC}+r)I^{TC}-r_{FC}{I}^{FC}-(r_I+r)I^V(t)$$

$$-r_{K}K(t)+r\hat{u}(t,x)+{\lambda }_F(-\hat{u}(t,x)-I^V(t)+I^{TC}-{\alpha }^B(t)P_t^B{)}^{-}]dt.$$

$${\alpha }_B=-\frac{{\theta }^B-\hat{u}(t,x)}{P_B},{\alpha }_C=-\frac{{\theta }^C-\hat{u}(t,x)}{P_C},$$

$$d\Pi =$$

$$+[-({\theta }^B-\hat{u}(t,x)){\lambda }_B-({\theta }^C-\hat{u}(t,x)){\lambda }_C+(r_{TC}+r)I^{TC}-r_{FC}{I}^{FC}-(r_I+r)I^V(t)$$

$$-r_{K}K(t)+r\hat{u}(t,x)+{\lambda }_F({\theta }^B-I^V(t)+I^{TC}{)}^{-}]dt.$$

$$L\hat{u}(t,x)=$$

$$f(t,x,\hat{u}(t,x))=$$

$$-(r_I+r)I^V(t)-r_{K}K(t)+r\hat{u}(t,x)+{\lambda }_F({\theta }^B-I^V(t)+I^{TC}{)}^{-}.$$

$$\displaystyle Y_{t}=\$$

$$-{\int }_t^T{\int }_{\mathbb{R}}{V}_s(q)d{\tilde{N}}_s(s,X_s,q),$$

Full text

Keywords: Fast Fourier Transform, CVA, XVA, BSDE, characteristic function

1 Introduction

After the financial crisis in 2007, it was recognized that Counterparty Credit Risk (CCR) poses a substantial risk 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 counterparty default event for over-the-counter (OTC) uncollateralized derivatives. The CVA is the expected loss arising from a default by the counterparty and can be defined as the difference between the risky value and the current risk-free value of a derivatives contract. CVA is calculated and hedged in the same way as derivatives by many banks, therefore having efficient ways of calculating the value and the Greeks of these adjustments is important.

One common way of pricing CVA is to use the concept of expected exposure, defined as the mean of the exposure distribution at a future date. Calculating these exposures typically involve computationally time-consuming Monte Carlo procedures, like nested Monte Carlo schemes or the more efficient least squares Monte Carlo method (LSM)[17]. Recently the Stochastic Grid Bundling method (SGBM)[14] was introduced as an improvement of the standard LSM. This method was extended to pricing CVA for Bermudan options in [10]. Another recently introduced alternative is the so-called finite-differences Monte Carlo method (FDMC) [6]. The FDMC method uses the scenario generation from the Monte Carlo method combined with finite-difference option valuation.

Besides CVA, many other valuation adjustments, collectively called XVA, have been introduced in option pricing in the recent years, causing a change in the way derivatives contracts are priced. For instance, a companies own credit risk is taken into account with a debt value adjustment (DVA). The DVA is the expected gain that will be experienced by the bank in the event that the bank defaults on its portfolio of derivatives with a counterparty. To reduce the credit risk in a derivatives contract, the parties can include a credit support annex (CSA), requiring one or both of the parties to post collateral. Valuation of derivatives under CSA was first done in [21]. A margin valuation adjustment (MVA) arises when the parties are required to post an initial margin. In this case the cost of posting the initial margin to the counterparty over the length of the contract is known as MVA. Funding value adjustments (FVA) can be interpreted as a funding cost or benefit associated to the hedge of market risk of an uncollateralized transaction through a collateralized market. While there is still a debate going on about whether to include or exclude this adjustment, see [13] and [12] for an in-depth overview of the arguments, most dealers now seem to indeed take into account the FVA. The capital value adjustment (KVA) refers to the cost of funding the additional capital that is required for derivative trades. This capital acts as a buffer against unexpected losses and thus, as argued in [11], has to be included in derivative pricing.

For pricing in the presence of XVA, one needs to redefine the pricing partial differential equation (PDE) by constructing a hedging portfolio with cashflows that are consistent with the additional funding requirements. This has been done for unilateral CCR in [21], bilateral CCR and XVA in [3] and extended to stochastic rates in [15]. This results in a non-linear option valuation PDE.

Non-linear PDEs can be solved by e.g. finite-difference methods or the LSM for solving the corresponding backward stochastic differential equation (BSDE). In [22] an efficient forward simulation algorithm that gives the solution of the non-linear PDE as the optimum over solutions of related but linear PDEs is introduced, with the computational cost being of the same order as one forward Monte Carlo simulation. The downside of these numerical methods is the computational time that is required to reach an accurate solution. An efficient alternative might be to use Fourier methods for solving the (non-)linear PDE or related BSDE, such as the COS method, as was introduced in [8], extended to Bermudan options in [9] and to BSDEs in [23]. In certain cases the efficiency of these methods is further increased due to the ability to additionally use the fast Fourier transform (FFT).

In this paper we consider an exponential Lévy-type model with a state-dependent jump measure and propose an efficient Fourier-based method to solve for Bermudan derivatives, including options and swaptions, with XVA. We derive, in the presence of state-dependent jumps, a non-linear partial integro-differential equation (PIDE) and its corresponding BSDE for an OTC derivative between a bank $B$ and its counterparty $C$ in the presence of CCR, bilateral collateralization, MVA, FVA and KVA, by setting up a hedging portfolio in which we focus on hedging the default risks and take into account the different rates associated with different types of lending. We extend the Fourier-based method known as the BCOS method, developed in [23], to solve the BSDE under Lévy models with non-constant coefficients. As this method requires the knowledge of the characteristic function of the forward process, which, in the case of the Lévy process with variable coefficients, is not known, we will use an approximation of the characteristic function obtained by the adjoint expansion method developed in [19], [18] and extended to the defaultable Lévy process with a state-dependent jump measure in [2]. Compared to other state-of-the-art methods for calculating XVAs, like Monte Carlo methods and PDE solvers, our method is more efficient and/or flexible. The efficiency is both due to the availability of the characteristic function in closed form through the adjoint expansion method and the fast convergence of the COS method. Furthermore we propose an alternative Fourier-based method for explicitly pricing the CVA term in case of unilateral CCR for Bermudan derivatives under the local Lévy model. The advantage of this method is that is allows us to use the FFT, resulting in a fast and efficient calculation. 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évy models with non-constant coefficients. In Section 3 we derive the non-linear PIDE and corresponding BSDE for pricing contracts under XVA. In Section 4 we propose the Fourier-based method for solving this BSDE and in Section 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 examples validating the accuracy and efficiency of the proposed methods.

2 The model

We consider a defaultable asset $S_{t}$ whose risk-neutral dynamics are given by

[TABLE]

where $d\tilde{N}_{t}(t,X_{t-},dq)$ is a compensated random measure with state-dependent Lévy measure

[TABLE]

The default time $\zeta$ of $S_{t}$ is defined in a canonical way as the first arrival time of a doubly stochastic Poisson process with local intensity function $\gamma(t,x)\geq 0$, and $\varepsilon\sim\mathrm{Exp}(1)$ and is independent of $X_{t}$. This way of modeling default is also considered in a diffusive setting in [5] and for exponential Lévy models in [4]. Thus, our model includes a local volatility function, a local jump measure, and a default probability which is dependent on the underlying. We define the filtration at time $t$ of the market observer to be $\mathcal{G}_{t}=\mathcal{F}^{X}_{t}\vee\mathcal{F}^{D}_{t}$, where $\mathcal{F}^{X}_{t}$ is the filtration generated by $X$ upto time $t$ and $\mathcal{F}_{t}^{D}:=\sigma(\{\zeta\leq u\},u\leq t)$, for $t\geq 0$, is the filtration of the default. Using this definition of default, the probability of default is

[TABLE]

We assume furthermore

[TABLE]

Imposing that the discounted asset price $\tilde{S}_{t}:=e^{-rt}S_{t}$ is a $\mathcal{G}$-martingale under the risk-neutral measure, we get the following restriction on the drift coefficient:

[TABLE]

with $r$ being the risk-free (collateralized) rate. In the whole of the paper we assume deterministic, constant interest rates, while the derivations can easily be extended to time-dependent rates. The integro-differential operator of the process is given by (see e.g. [20])

[TABLE]

3 XVA computation

Consider a bank $B$ and its counterparty $C$, both of them might default. Assume they enter into a contract paying $\Phi(S_{t})$ at maturity. Let $\phi(x)=\Phi(e^{x})$, and assume the risk-neutral dynamics of the underlying as in (3) with the drift given by (7). Define $\hat{u}(t,x)$ to be the value to the bank of the (default risky) portfolio with valuation adjustments referred to as XVA and $u(t,x)$ to be the risk-free value. Note that the difference between these two values is called the total valuation adjustment and in our setting this consists of

[TABLE]

The risk-free value $u(t,x)$ solves a linear PIDE:

[TABLE]

where $L$ is given in (8). Assuming the dynamics in (3), this linear PIDE can be solved with the methods presented in [2].

3.1 Derivative pricing under CCR and bilateral CSA agreements

In [3], the authors derive an extension to the Black-Scholes PDE in the presence of a bilateral counterparty risk in a jump-to-default model with the underlying being a diffusion, using replication arguments that include the funding costs. In [15] this derivation is extended to a multivariate diffusion setting with stochastic rates in 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 are obliged to hold regulatory capital, the cost of which is the KVA and face the costs of funding uncollateralized positions through collateralized markets, known as FVA. Both [3] and [15] extend the approach of [21], in which unilateral collateralization was considered. We extend their approach to derive the value of $\hat{u}(t,x)$ when the underlying follows the jump-diffusion defined in (3). We assume a one-dimensional underlying diffusion and consider all rates to be deterministic and, for ease of notation, constant. We specify different rates, defined in Table 3.1, for different types of lending.

Assume that the parties $B$ and $C$ enter into a derivative contract on the spot asset that pays the bank $B$ the amount $\phi(X_{T})$ at maturity $T$. The value of this derivative to the bank at time t is denoted by $\hat{u}(t,x,\mathcal{J}^{B},\mathcal{J}^{C})$ and depends on the value of the underlying $X$ and the default states $\mathcal{J}^{B}$ and $\mathcal{J}^{C}$ of the bank $B$ and counterparty $C$, respectively. Define $I^{TC}$ to be the initial margin posted by the bank to the counterparty, $I^{FC}$ the initial margin posted by the counterparty to the bank and $I^{V}(t)$ to be the variation margin on which a rate $r_{I}$ is paid or received. The initial margin is constant throughout the duration of the contract. 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 counterparty or the bank, the value of the derivative to the bank $\hat{u}(t,x)$ is determined with a mark-to-market rule $M$, which may be equal to either the derivative value $\hat{u}(t,x,0,0)$ prior to default or the risk-free derivative value $u(t,x)$, depending on the specifications in the ISDA master agreement. Denote by $\tau^{B}$ and $\tau^{C}$ the random default times of the bank and the counterparty respectively. We will use the notation $x^{+}=\max(x,0)$ and $x^{-}=\min(x,0)$. In a situation in which the counterparty defaults, the bank is already in the possession of $I^{V}+I^{FC}$. If the outstanding value $M-(I^{V}+I^{FC})$ is negative, the bank has to pay the full amount $(M-I^{V}-I^{FC})^{-}$, while if the contract has a positive value to the bank, it will recover only $R_{C}(M-I^{V}-I^{FC})^{+}$. Using a similar argument in case the bank defaults, we find the following boundary conditions:

[TABLE]

so that the portfolio value at default is given by

[TABLE]

with $\tau=\min(\tau^{B},\tau^{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 zero otherwise. Assume the dynamics for $P^{B}_{t}$ and $P_{t}^{C}$ to be given by $P_{t}^{B}={\mathds{1}}_{\{\tau^{B}>t\}}e^{r_{B}t}$ and $P_{t}^{C}={\mathds{1}}_{\{\tau^{C}>t\}}e^{r_{C}t}$, so that

[TABLE]

with $\mathcal{J}_{t}^{B}={\mathds{1}}_{\tau^{B}\leq t}$ and $\mathcal{J}_{t}^{C}={\mathds{1}}_{\tau^{C}\leq t}$, where the default times $\tau^{B}$ and $\tau^{C}$ are defined in a canonical way as the first arrival time of a doubly stochastic Poisson process with intensity functions $\gamma^{B}$ and $\gamma^{C}$, respectively (see also the definition of the defaultable asset in (3)). We define the market interest rates for $B$ and $C$ to be $r_{B}=r+\gamma^{B}$ and $r_{C}=r+\gamma^{C}$, so that by the usual arguments (see, for instance, [16, Section 2.2]) the discounted bonds $e^{-rt}P_{t}^{B}$ and $e^{-rt}P_{t}^{C}$ are martingales under the risk-neutral measure.

We construct a hedging portfolio consisting of the shorted derivative, $\alpha_{C}$ units of $P^{C}$, $\alpha_{B}$ units of $P^{B}$ and $g$ units of cash:

[TABLE]

In other words, since we assume both the underlying asset process and the tradeable bonds $P_{B}$ and $P_{C}$ to be risk-neutral, we focus on hedging the risk arising from the defaults of both $B$ and $C$ by means of the default-risky bonds.

If the value of the derivative is positive to $B$, it will incur a cost at the counterparties’ default. To hedge this, $B$ shorts $P^{C}$, i.e. $\alpha_{C}\leq 0$. If we assume $B$ can borrow the bond close to the risk-free rate $r$ (i.e. no haircut) through a repurchase agreement, it will incur financing costs of $r\alpha_{C}(t)P_{t}^{C}dt$. The cashflows from the collateralization follow from the rate $r_{TC}$ received and $r_{FC}$ 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}<0$, and the bank pays collateral respectively. From holding the regulatory capital we incur a cost of $r_{K}K(t)$. Finally, the rates $r$ and $r_{F}$ are respectively received or paid on the surplus cash in the account. This cash consists of the gap between the shorted derivative value and the collateral and the cost of buying $\alpha_{B}$ bonds $P^{B}$ in order for $B$ to hedge its own default, i.e. $-\hat{u}(t,x)-I^{V}(t)+I^{TC}-\alpha_{B}(t)P_{t}^{B}$. Thus, the total change in the cash account is given by

[TABLE]

Note that this is in contrast with the change in cash in a portfolio without the XVA arising from the different types of funding, i.e. where we assume the cash in the portfolio simply earns the risk-free rate

[TABLE]

Assuming the portfolio is self-financing we have

[TABLE]

Applying Itô’s Lemma to $\hat{u}(t,x)$ gives us:

[TABLE]

with the operator $L$ as in (8). Thus, we find,

[TABLE]

By choosing

[TABLE]

we hedge the jump-to-default risk in the hedging portfolio, i.e.,

[TABLE]

Then, using the fact that the portfolio has to satisfy the martingale condition in the risk-neutral world, i.e. $\mathbb{E}[d\Pi]=0$, we find the non-linear pricing PIDE to be

[TABLE]

where we have defined

[TABLE]

3.2 BSDE representation

In this section we will cast the PIDE in (30) in the form of a Backward Stochastic Differential Equation. In the methods where we make use of BSDEs we assume $\gamma(t,x)=0$. We begin by recalling the non-linear Feynman-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 (3). We assume $\mu$, $\sigma$ and $a$ to be Lipschitz continuous in $x$ and additionally $|a(t,x)|\leq K$. Consider the BSDE

[TABLE]

where the generator $f$ is continuous and satisfies the Lipschitz condition in the space variables, $\delta$ is a measurable, bounded function and the terminal condition $\phi(x)$ is measurable and Lipschitz continuous. Consider the non-linear PIDE

[TABLE]

If the PIDE in (35) has a solution $u(t,x)\in C^{1,2}$, the FBSDE in (34) has a unique solution $(Y_{t},Z_{t},V_{t}(q))$ that can be represented as

[TABLE]

for all $s\in[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 (30) reads

[TABLE]

where we have defined the driver function to be

[TABLE]

3.3 A simplified driver function

Following [11], one can derive that the KVA is a function of trade properties (i.e. maturity, strike) and/or the exposure at default, which in turn is a function of the portfolio value, so that the cost of holding the capital can be rewritten as $r_{K}K(t)=r_{K}c_{1}\hat{u}(t,x),$ with $c_{1}$ being a function of the trade properties. The collateral is paid when the portfolio has a negative value, and received when the collateral has a positive value. Assuming the collateral is a multiple of the portfolio value we have $I^{V}(t)=c_{2}\hat{u}(t,x)$, where $c_{2}$ is some 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 function reduces to $f(t,\hat{u}(t,x))=r_{u}(t)\max(\hat{u}(t,x),0)$, for a function $r_{u}$ here left unspecified. In this case the BSDE is similar to the one considered in [22].

4 Solving FBSDEs

In this section we extend the BCOS method from [23] to solving FBSDEs under local Lévy models with variable coefficients and jumps (without default, i.e. $\gamma(t,x)=0$). The conditional expectations resulting from the discretization of the FBSDE are approximated using the COS method. This requires the characteristic function, which we approximate using the Adjoint Expansion Method of [19] and [2].

4.1 Discretization of the BSDE

Consider the forward process $X_{t}$ as in (3) and the BSDE $Y_{t}$ as in (39) with a more general driver function $f(t,x,y,z)$. Define a partition $0=t_{0}<t_{1}<...<t_{N}=T$ of $[0,T]$ with a fixed time step $\Delta t=t_{n+1}-t_{n}$, for $n=N-1,...0$. Rewriting the set of FBSDEs we find,

[TABLE]

One can obtain an approximation of the process $Y_{t}$ by taking conditional expectations with respect to the underlying filtration $\mathcal{G}_{n}$, using the independence of $W_{t}$ and $\tilde{N}_{t}(t,X_{t-},dq)$ and by approximating the integrals that appear with a theta-method, as first done in [24] and extended to BSDEs with jumps in [23]:

[TABLE]

Let $\Delta W_{s}:=W_{s}-W_{n}$ for $t_{n}\leq s\leq t_{n+1}$. Multiplying both sides of equation (43) by $\Delta W_{n+1}$, taking conditional expectations and applying the theta-method gives

[TABLE]

Since in our scheme the terminal values are functions of time $t$ and the Markov process $X$, it is easily seen that there exist deterministic functions $y(t_{n},x)$ and $z(t_{n},x)$ so that

[TABLE]

The functions $y(t_{n},x)$ and $z(t_{n},x)$ are obtained in a backward manner using the following scheme

[TABLE]

where we have simplified notations with

[TABLE]

In the case $\theta_{1}>0$ we obtain an implicit dependence on $y(t_{n},x)$ in (50) and we use $P$ Picard iterations starting with initial guess $\mathbb{E}_{n}[y(t_{n+1},X_{n+1})]$ to determine $y(t_{n},x)$.

4.2 The characteristic function

Is it well-known (see, for instance, [16, Section 2.2]) that the risk-free pre-default price $u(t,x)$ of a European option on the defaultable asset $S_{t}$ with maturity $T$ and payoff $\phi(X_{T})$ is given by

[TABLE]

in the measure corresponding to the dynamics in (3). Thus, in order to compute the price of an option, we must evaluate functions of the form

[TABLE]

Under standard assumptions, by the Feynman-Kac theorem, $v$ can be expressed as the classical solution of the following Cauchy problem

[TABLE]

with $L$ as in (8).

The function $v$ in (55) can be represented as an integral with respect to the transition distribution of the defaultable log-price process $\log S_{t}$:

[TABLE]

where $\Gamma(t,x;T,dy)$ is the Green’s function of the PIDE in (56) and we say that its Fourier transform

[TABLE]

is the characteristic function of $\log S$. Following [19] and [2] we expand the state-dependent coefficients

[TABLE]

around some point $\bar{x}$. The coefficients $s(t,x)$, $\gamma(t,x)$ and $a(t,x)$ are assumed to be continuously differentiable with respect to $x$ up to order $n\in\mathbb{N}$.

Introduce the $n$th-order approximation of $L$ in (8):

[TABLE]

where

[TABLE]

and

[TABLE]

The basepoint $\bar{x}$ is a constant parameter which can be chosen freely. In general the simplest choice is $\bar{x}=x$ (the value of the underlying at initial time $t$).

Assume for a moment that $L_{0}$ has a fundamental solution $G^{0}(t,x;T,y)$ that is defined as the solution of the Cauchy problem

[TABLE]

In this case we define the $n$th-order approximation of $\Gamma$ as

[TABLE]

where, for any $k\geq 1$ and $(T,y)$, $G^{k}(\cdot,\cdot;T,y)$ is defined recursively through the following Cauchy problem

[TABLE]

Notice that

[TABLE]

Correspondingly, the $n$th-order approximation of the characteristic function $\hat{\Gamma}$ is defined to be

[TABLE]

Now, by transforming the simplified Cauchy problems into adjoint problems and solving these in the Fourier space we find

[TABLE]

with

[TABLE]

[TABLE]

where $\bar{\nu}(dq)=\nu(-dq)$.

Remark 3**.**

After some algebraic manipulations it can be shown, see [2], that the characteristic function approximation of order $n$ is a function of the form

[TABLE]

where the coefficients $g_{n,k}$, with $0\leq k\leq n$, depend only on $t,T$ and $\xi$, but not on $x$. The approximation formula can thus always be split into a sum of products of functions depending only on $\xi$ and functions that are linear combinations of $(x-\bar{x})^{m}e^{i\xi x}$, $m\in\mathbb{N}_{0}$.

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,1$ and assume the coefficients $s(t,x)$, $\gamma(t,x)$ and $a(t,x)$ are continuously differentiable with bounded derivatives up to order $n$. For the $n$th-order approximation $\Gamma^{(n)}(t,x;T,\xi)$, for any $\bar{x}\in\mathbb{R}$,

[TABLE]

Note that if $\bar{x}=x$, the bound reduces to $C(T,\xi)(T-t)^{n+1}$.

4.3 The COS formulae

The conditional expectations are approximated using the COS method, which was developed in [9] and applied to FBSDEs with jumps in [23]. The conditional expectations arising in the equations (50)-(51) are all of the form $\mathbb{E}_{n}[h(t_{n+1},X_{n+1})]$ or $\mathbb{E}_{n}[h(t_{n+1},X_{n+1})\Delta W_{n+1}]$. The COS formula for the first type of conditional expectation reads

[TABLE]

where $\sideset{}{{}^{\prime}}{\sum}$ denotes an ordinary summation with the first term weighted by one-half, $J>0$ is the number of Fourier-cosine coefficients we use, $H_{j}(t_{n+1})$ denotes the $j$th Fourier-cosine coefficients of the function $h(t_{n+1},x)$ and $\hat{\Gamma}\left(t_{n},x;t_{n+1},\xi\right)$ 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

[TABLE]

See [23] for the full derivations.

Remark 5**.**

Note that these formulas are obtained by using an Euler approximation of the forward process and using the 2nd-order approximation of the characteristic function of the actual process. We have found this to be more exact than using the characteristic function of the Euler process, which is equivalent to using just 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

[TABLE]

Due to the structure of the approximated characteristic function of the local Lévy process, see (71), the coefficients of the functions $z(t_{n+1},x)$ and the explicit part of $y(t_{n+1},x)$ can be computed using the FFT algorithm, as we do in Appendix A, because of the matrix in (110) being of a certain form with constant diagonals. In order to determine $F_{j}(t_{n+1})$, the Fourier-Cosine coefficient of the function

[TABLE]

due to the intricate dependence on the functions $z$ and $y$ we choose to approximate the integral in $F_{j}$ by a discrete 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+\frac{1}{2})\frac{b-a}{J}$ and $\Delta x=\frac{b-a}{J}$ we find, using the mid-point integration rule, the approximation

[TABLE]

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:

[TABLE]

where $c_{n}$ is the $n$th cumulant of log-price process $\log S$, as proposed in [8]. The cumulants are calculated using the 0th-order approximation of the characteristic function.

5 XVA computation for Bermudan derivatives

The method in Section 4 allows us to compute the XVA as in (9), consisting of CVA, DVA, MVA, KVA and FVA. In this section, we apply this method to computing Bermudan derivative values with XVA. The resulting method – the solution of the non-linear XVA PDE through a BSDE-type method – is an efficient alternative to finite-difference methods as well as to the Monte-Carlo based method developed in [22]. The efficiency is both due to the availability of the characteristic function in closed form through the adjoint expansion method and the fast convergence of the COS method. Furthermore, in finite difference methods complications may arise in the implementation of the scheme for jump diffusions. Since our proposed method works in the Fourier space, the jump component is easily handled by means of an additional term in the characteristic 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 the FFT, results in a particularly efficient computation.

5.1 XVA computation

Consider an OTC derivative contract between the bank $B$ and the counterparty $C$ on the underlying asset $S_{t}$ given by (3) with $\gamma(t,x)=0$ with a Bermudan-type exercise possibility: there is a finite set of so-called exercise moments $\{t_{1},...,t_{M}\}$ prior to the maturity, with $0\leq t_{1}<t_{2}<\cdots<t_{M}=T$. The payoff from the point-of-view of bank $B$ is given by $\phi(t_{m},X_{t_{m}})$. Denote $\hat{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 derivative with XVA and $M$ exercise dates $t_{1},...,t_{M}$ can be expressed by a backward recursion as

[TABLE]

and the continuation value solves the non-linear PIDE defined in (30)

[TABLE]

The derivative value is set to be $\hat{u}(t,x)=c(t,x)$ for $t\in]t_{m-1},t_{m}[$, and, if $t_{1}>0$, also for $t\in[0,t_{1}[$. The payoff function might take on various forms:

(Portfolio) Following [22], we can consider $X_{t}$ to be the process of a portfolio which can take on both positive and negative values. Then, when exercised at time $t_{m}$, bank $B$ receives the portfolio so that $\phi(t_{m},x)=e^{x}$. 2. 2.

(Bermudan option) In case the Bermudan contract is an option, the option value to the bank can not have a negative value for the bank. At the same time, in case of default of the bank itself, the counterparty loses nothing. In this case the framework simplifies to one with unilateral collateralization and default risk and the payoff at time $t_{m}$, if exercised, is given by $\phi(t_{m},x)=(K-e^{x})^{+}$ for a put and $\phi(t_{m},x)=(e^{x}-K)^{+}$ for a call with $K$ being the strike price. 3. 3.

(Bermudan swaptions) A Bermudan swaption is an option in which the holder, bank $B$, has the right to exercise and enter into an underlying swap with fixed end date $t_{M+1}$. If the swaption is exercised at time $t_{m}$ the underlying swap starts with payment dates $\mathcal{T}_{m}=\{t_{m+1},...,t_{M+1}\}$. Working under the forward measure corresponding to the last reset date $t_{M}$, the payoff function is given by

[TABLE]

where $N^{S}$ is the notional, $c_{p}=1$ for a payer swaption and $c_{p}=-1$ for a receiver swaption, $P(t_{m},t_{k},x)$ is the price of a ZCB conditional on $X_{t_{m}}=x$ and $S(t_{m},\mathcal{T}_{m},x)$ is the forward swap rate given by

[TABLE]

To solve for the continuation value we define a partition with $N$ steps $t_{m-1}=t_{0,m}<t_{1,m}<t_{2,m}<...<t_{n,m}<...<t_{N,m}=t_{m}$ between two exercise dates $t_{m-1}$ and $t_{m}$, with fixed time step $\Delta t_{n}:=t_{n+1,m}-t_{n,m}$. Applying the method developed in Section 4, we find the following time iteration for the continuation value:

[TABLE]

where we have defined

[TABLE]

and the Fourier-cosine coefficients are given by

[TABLE]

In order to determine the function $c(t_{n},x)$, we will perform $P$ Picard iterations. To evaluate the coefficients with a DCT we need to compute the integrands $c(t_{n+1,m},x)$ and $f(t_{n+1,m},x,c(t_{n+1,m},x))$ on the equidistant $x$-grid with $x_{i}$, for $i=0,...,J-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 (86). The matrix-vector product in the formula results in a computational time of order $O(J^{2})$.

Remark 7** (Convergence of the Picard iterations).**

A Picard iteration is used to find the fixed-point $c$ of $c=\Delta t\theta_{1}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 parts of 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. $\exists$ $L^{Lipz}$ such that $|f(t,x,c_{1})-f(t,x,c_{2})|\leq L^{Lipz}|c_{1}-c_{2}|$, and $\Delta t_{n}$ is small enough such that $\Delta t\theta_{1}L^{Lipz}<1$, a unique fixed-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,0)$, and this is Lipschitz continuous with $L^{Lipz}=1$. Thus for $\Delta 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 as in Algorithm 1 in Figure 5.1. The total computational time for the algorithm is of order

[TABLE]

consisting of the computation for $M\cdot 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-vector multiplications in the formulas for $c(t_{n,m},x)$ and $z(t_{n,m},x)$ of $O(J^{2})$, initialization of the Picard method with $\mathbb{E}_{n}[c(t_{n+1},X_{n+1}]$ in $O(J^{2})$ operations, computation of the $P$ Picard approximations for $c(t_{n,m},x)$ in $O(PJ)$ and computing the Fourier coefficients $F_{j}(t_{n})$ and $C_{j}(t_{n})$ with the DCT in $O(J\log_{2}J)$ operations.

5.2 An alternative for CVA computation

In this section we present an efficient alternative way of calculating the CVA term in (9) in the case of unilateral CCR using a Fourier-based method. Due to the ability of using the FFT this method is considerably faster for computing the CVA than the method presented in Section 5.1. We use the definition of CVA at time $t$ given by

[TABLE]

where $u(t,X_{t})$ is as usual the default-free value of the Bermudan option ($\gamma(t,x)=0$), while $\hat{u}(t,X_{t})$ is the value including default ($\gamma(t,x)\neq 0$). We consider the model as defined in (3). We will compute $u(t,X_{t})$ and $\hat{u}(t,X_{t})$ using the COS method and the approximation of the characteristic function (as derived in Section 4.3), without default and with default, respectively. In case of a default the payoff becomes zero. Note that the risky option value $\hat{u}(t,x)$ computed with the characteristic function for a defaultable underlying corresponds exactly to the option value in which the counterparty might default, with the probablity of default, $PD(t)$, defined as in (5). Thus, in this case we have unilateral CCR and $\zeta=\tau_{C}$, the default time of the counterparty.

Using the definition of the defaultable $S_{t}$, it is well-known (see, for instance, [16, Section 2.2]) that the risky no-arbitrage value of the Bermudan option on the defaultable asset $S_{t}$ at time $t$ is given by

[TABLE]

Remark 8** (Wrong-way risk).**

By allowing the dependence of the default intensity on the underlying, a simplified form of wrong-way risk is already incorporated into the CVA valuation.

For a Bermudan put option with strike price $K$, we simply have $\phi(t,x)=\left(K-x\right)^{+}$. By the dynamic programming approach, the option value can be expressed by a backward recursion as

[TABLE]

and

[TABLE]

Thus to find the risky option price $\hat{u}(t,X_{t})$ one uses the defaultable asset with $\gamma(t,x)$ representing the default intensity of the counterparty and in order to get the default-free value $u(t,X_{t})$ one uses the default-free asset by setting $\gamma(t,x)=0$. The CVA adjustment is calculated as the difference between the two. Both $\hat{u}(t,x)$ and $u(t,x)$ are calculated using the approximated characteristic function and the COS method applied to the continuation value [2]. Due to the characteristic function being of the form (71), we are able to use the FFT in the matrix-vector multiplication when computing the continuation values of the Bermudan option with and without default, reducing this operation from $O(J^{2})$ to $O(J\log_{2}J)$. For more details, we refer to Appendix A. The total complexity of the calculation of the CVA value for a Bermudan option with $M$ exercise dates is then $O(MJ\log_{2}J)$. Comparing this to (90), in which the most time-consuming operations were indeed the matrix-vector products of order $O(J^{2})$ that resulted from the computation of the functions on the $x$-grid of size $J$, we conclude that the method for CVA computation is indeed significantly faster due to the ability of using the FFT.

5.2.1 Hedging CVA

In practice CVA is hedged and thus practitioners require efficient ways to compute the sensitivity of the CVA with respect to the underlying. The widely used bump- and revalue- method, while resulting in precise calculations, might be slow to compute. Using the Fourier-based approach we find explicit formulas allowing for 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

[TABLE]

[TABLE]

where $V_{k}^{d}$ and $V_{k}^{r}$ are the Fourier-cosine coefficients with the defaultable and default-free characteristic function terms, $g_{n,h}^{d}$ and $g_{n,h}^{r}$, respectively.

6 Numerical experiments

In this section we present numerical examples to justify the accuracy of the methods in practice. We compute the XVA with the method presented in Section 5.1 and the CVA in the case of unilateral CCR with the method from Section 5.2, which we show is more efficient for cases in which one only needs to compute the CVA. We compare the results of solving the BSDE with the COS method and the adjoint expansion of the characteristic function to the values obtained by using a least-squares Monte-Carlo method for computing the 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 the second-order approximation of the characteristic function. We have found this to be sufficiently accurate by numerical experiments and theoretical error estimates. The formulas for the second-order approximation are simple, making the methods easy to implement.

6.1 A numerical example for XVA

Here, we check the accuracy of the method from Section 5.1. We will compute the Bermudan option value with XVA using a simplified driver function given by $f(t,\hat{u}(t,x))=-r\max(\hat{u}(t,x),0)$. Our method is easily extendible to the driver function in Section 3.2. Consider $X_{t}$ to be a portfolio process and the payoff, if exercised at time $t_{m}$, to be given by $\Phi(t_{m},x)=x$. In this case the value we can receive at every exercise date is the value of the portfolio. Consider the model in Section 2 without default, with a local jump measure and a local volatility function with CEV-like dynamics and Gaussian jumps defined by

[TABLE]

We assume the following parameters in equations (101)-(102), unless otherwise mentioned: $b=0.15$, $\beta=-2$, $\lambda=0.2$, $\delta=0.2$, $m=-0.2$, $r=0.1$, $K=1$ and $X_{0}=0$ (so that $S_{0}=1$). In the LSM the number of time steps is taken to be 100 and we simulate $10^{5}$ paths. In the COS method we take $J=256$, $\theta_{1}=0.5$ and $N=10$, $M=10$, making the total number of time steps $N\cdot M=100$. The truncation range is determined as in (78) with $L=10$. Due to the state-dependent coefficients in the underlying dynamics in (101)-(102) we use the approximated characteristic function as derived in Section 4.2 with the second-order approximation, i.e. $\hat{\Gamma}^{(2)}(t,x;T,\xi)$ and take $\bar{x}=x$, where $x=\{x_{i}\}_{i=0}^{J-1}$. Note that we thus compute the values, including those of the characteristic function, on the complete $x$-grid. In the final iteration when computing $\hat{u}(t_{0},X_{0})$ we use $\bar{x}=X_{0}$.

In Table 6.1 we analyse the error in the approximation of $\hat{u}(t_{0},X_{0})$ with $S_{0}=0.4$ for different values of the discretization parameter $N$ and the number of grid points (and Fourier-cosine coefficients) $J$. We compare the approximated COS value to the 95% confidence interval obtained by a LSM. Accurate results are quickly obtained for small values of both $J$ and $N$. In Figure 6.1 we plot the upper bound of the 95% confidence interval of the absolute error in the approximation for varying $J$ and $N$. We observe approximately a linear convergence and note that the error stops decreasing at some point for increasing values of $J$ and $N$. This can be due to the error being dominated by the approximated characteristic function. In particular we observe that $J=32$ and $N=10$ seem to be sufficient parameters to achieve a satisfactory accuracy in the approximation.

The results for $\hat{u}(t_{0},X_{0})$ of the COS approximation method compared to a 95% confidence interval of the value obtained through a LSM are presented in Table 6.2. These results show that our method is able to solve non-linear PIDEs accurately. The CPU time of the approximating method depends on the number of time steps $M\cdot N$ and is approximately $5\cdot(N\cdot M)$ ms.

6.2 A numerical example for CVA

In this section we validate the accuracy of the method presented in Section 5.2 and compute the CVA in the case of unilateral CCR under the model dynamics given in Section 2 with a local jump measure and a local volatility function with CEV-like dynamics, Gaussian jumps defined by defined as in (102) and a local default function $\gamma(x)=ce^{\beta x}$. We assume the same parameters as in Section 6.2, except $r=0.05$ and we take $c=0.1$ in the default function. In the LSM the number of time steps is taken to be 100 and we simulate $10^{5}$ paths. In the COS method we take $L=10$ and $J=100$. Again, due to the state-dependent coefficients in the underlying dynamics we use the approximated characteristic function as derived in Section 4.2 with the second-order approximation, i.e. $\hat{\Gamma}^{(2)}(t,x;T,\xi)$ and take $\bar{x}=X_{0}$.

The results for the CVA valuation with the FFT-based method and with LSM are presented in Table 6.3. 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\cdot M$ ms, thus more efficient than the computation of the XVA with the method in Section 5.1. The optimal exercise boundary in Figure 6.2 shows that the exercise region becomes larger when the probability of default increases; this is to be expected: in case of the default probability being greater, the option of exercising early is more valuable and used more often.

7 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 different rates for different types of lending. This resulted in the PIDE in (30) and its corresponding BSDE (39). We propose to solve the BSDE using a Fourier-cosine method for the resulting conditional expectations and an adjoint expansion method for determining an approximation of the characteristic function of the local Lévy model in (3). This approach is extended to Bermudan option pricing in Section 5.1. In Section 5.2 we presented an alternative for computing the CVA term in the case of unilateral collateralization (as is the case when the derivative is an option) without the use of BSDEs. This results in an even more efficient method due to the ability to use the FFT. We verify the accuracy of both methods in Sections 6.1 and 6.2 by comparing it to a LSM and conclude that the method from Section 5.1 is able to achieve a rapid convergence and gives, already for small values of the discretization parameters an accurate result. The alternative method for CVA computation from Section 5.2 is indeed more efficient than the BSDE method for computing just the CVA term.

Acknowledgments

We thank two anonymous referees for the comments and suggestions that have improved the quality of this paper. This research is supported by the European Union in the the context of the H2020 EU Marie Curie Initial Training Network project named WAKEUPCALL.

Appendix A The COS formulae

Let, as usual, $J$ denote the number of Fourier-cosine coefficients. Remembering that the expected value $c(t,x)$ in (94) can be rewritten in integral form, we have

[TABLE]

where, $v(t_{m},y)$ can be either $u(t_{m},y)$ or $\hat{u}(t_{m},y)$. Then we use the Fourier-cosine expansion to get the approximation:

[TABLE]

with $\phi(t,x)=\left(K-e^{x}\right)^{+}$.

We can recover the coefficients $\left(V_{j}(t_{m})\right)_{j=0,1,...,J-1}$ from $\left(V_{j}(t_{m+1})\right)_{j=0,1,...,J-1}$. To this end, we split the integral in the definition of $V_{j}(t_{m})$ into two parts using the early-exercise point $x_{m}^{*}$, which is the point where the continuation value is equal to the payoff, i.e. $c(t_{m},x_{m}^{*})=\phi(t_{m},x_{m}^{*})$; this point can easily be found by using the Newton method. Thus, we have

[TABLE]

where

[TABLE]

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}^{*}\leq\log K$, and by inserting the approximation (104) for the continuation value into the formula for $C_{j}(t_{m},x_{m}^{*})$ have the following coefficients $\hat{C}_{j}$ for $m=M-1,M-2,...,1$:

[TABLE]

From (71) we know that the $n$th-order approximation of the characteristic function is of the form:

[TABLE]

where the coefficients $g_{n,h}(t,T,\xi)$, with $0\leq k\leq n$, depend only on $t,T$ and $\xi$, but not on $x$.

Remark 9** (The defaultable and default-free characteristic functions).**

To find $u(t,x)$ we use

[TABLE]

the characteristic function with $\gamma(t,x)=0$. For $\hat{u}(t,x)$ we use

[TABLE]

where $\gamma(t,x)$ is chosen to be some specified function.

Using (71) we can write the Fourier coefficients of the continuation value in vectorized form as:

[TABLE]

where $\mathbb{V}(t_{m+1})$ is the vector $[V_{0}(t_{m+1}),...,V_{J-1}(t_{m+1})]^{T}$ and $\mathcal{M}^{h}(x_{m}^{*},b)\Lambda^{h}$ is a matrix-matrix product with $\mathcal{M}^{h}$ a matrix with elements $\{M_{k,j}^{h}\}_{k,j=0}^{J-1}$ defined as

[TABLE]

and $\Lambda^{h}$ is a diagonal matrix with elements

[TABLE]

One can show, see [2], that the resulting matrix $\mathcal{M}^{h}$ is a sum of a Hankel and Toeplitz matrix and thus the resulting matrix vector product can be calculated using a FFT.