# CVA in fractional and rough volatility models

E. Alòs111Dept. of Economics and Business, University Pompeu Fabra, Barcelona,
<EMAIL_ADDRESS>F. Antonelli222DISIM, University of L’Aquila,
<EMAIL_ADDRESS>A. Ramponi333Dept. Economics and Finance,
University of Roma - Tor Vergata<EMAIL_ADDRESS>S.
Scarlatti444Dept. Economics and Finance, University of Roma - Tor Vergata,
<EMAIL_ADDRESS>
###### Abstract

In this work we present a general representation formula for the price of a
vulnerable European option, and the related CVA in stochastic (either rough or
not) volatility models for the underlying’s price, when admitting correlation
with the default event. We specialize it for some volatility models and we
provide price approximations, based on the representation formula. We study
numerically their accuracy, comparing the results with Monte Carlo
simulations, and we run a theoretical study of the error. We also introduce a
seminal study of roughness influence on the claim’s price.

Keywords: Credit Value Adjustment; Vulnerable Options; Rough volatility
models; Intensity approach.

##  1 Introduction

In recent years, the interest in including correctly the effects due to
default risks in the derivatives’ evaluation has grown immensely, generating a
whole new field in mathematical finance. As it is usually the case, the
seminal papers work on simple market models (see [26], [24] [16], [29]) built
the way towards the extension to more complex and realistic ones. We refer the
reader to [10], [11], [12], [13], and [14] for a sample of the interest and
width of the topic.

Surely, considering stochastic volatility for derivatives subject to credit
risks is an important modeling feature, even more so under volatility
roughness that better reflects the market behavior, as much of the recent
literature shows, and we refer the reader to [6], [3], [9], [20] just to quote
some.

Inevitably, including credit risk features and stochastic volatility in the
same model makes the picture more complex and harder to manage. Often, this
complexity leaves Monte Carlo simulations as the only option to evaluate
derivatives’ prices in such a framework, paying a remarkable price in terms of
computational times.

Therefore, extensive work was done to develop representation formulas that
might lead to handier expressions and to finding alternative computational
methods, as in [1], and more recently in [8] and [7].

Lately, see for instance [21], deep learning techniques have been applied to
evaluate a derivative depending on multiple risks, whose price is corrected by
a family of adjustments that go under the acronym of XVA. This approach is
based on the price characterization as solution of a Backward Stochastic
Differential Equation and the application of machine learning to those
equations. Nevertheless, those algorithms though highly efficient in the
predictive phase, remain computationally costly in the learning phase still
based on Monte Carlo simulations. Moreover, the network architecture (e.g. the
number of hidden layers and the number of neurons per layer), as well as the
choice of the hyperparameters which characterize the implementation of the
learning algorithm (the activation function, the minimization procedure, the
mini-batch size, etc.) certainly require very careful programming work and
skills.

Here, we focus our attention on vulnerable European options, that are options
subject to some default event concerning the solvability of the issuer, and we
provide a general representation formula for the price and the required Credit
Value Adjustment (CVA) for stochastic (either rough or not) volatility market
models, when correlations among the driving processes are present.

The underlying idea is inherited from the papers by Alòs et al. ([1], [5]),
where the authors remark that those prices have to depend upon the integrated
mean variance (or the zero-strike variance swap process) and the Black &
Scholes pricing function, deriving a price representation by applying either
anticipative or non-anticipative Stochastic Calculus. To extend that technique
to a model including also the adjustment coming from the default
probabilities, we use the so called “intensity approach” ([11], [12] just to
quote some), that we assume to be a diffusion. Correlations among all the
processes, underlying, volatility and intensity, are admitted, making much
harder to exploit the properties of the single processes (such as affinity).
In particular, the intensity-underlying correlation is related to the so-
called right/wrong-way risk.

Exploiting the fact that the vulnerable option final value is the same as in
the default-free case, the explicit knowledge of the Black & Scholes pricing
function, and anticipating Stochastic Calculus, we are able to provide a
general and rather manageable representation formula for the CVA (and
consequently the adjusted price) in terms of the correlations between the
asset price and, respectively, the stochastic volatility and the intensity
process. This representation turns out to be quite general, including at once
both classical and rough stochastic volatility models.

Further, on the basis of that representation, we develop a manageable
approximation formula, which never employs the computationally costly Monte
Carlo simulations. To achieve this goal, we first notice that the zero-strike
variance swap process is the optional projection of the integrated mean
variance onto the market filtration, which allows developing the main term in
a non-anticipating fashion by classical Itô formula. Finally, freezing some
terms at the initial time, and using the well-known Clark-Ocone-Haussman
formula, we are able to write an approximation that makes the covariation
processes explicit.

When choosing the intensity in the class of the affine models and the
volatility either a diffusion or a rough one, the formula becomes quite easily
computable, as we show by selecting a CIR default intensity and a volatility
following a Heston, SABR or Rough Bergomi model. In those specialized models,
the approximation takes the shape of a “first-order” expansion with respect to
the correlation parameters among the Brownian motions driving the processes
dynamics.

Finally, we briefly discuss a theoretical error estimate procedure and we test
our method accuracy and efficiency by running a numerical study of the
rBergomi price model coupled with a CIR default intensity model, taking the
corresponding Monte Carlo price as benchmark. The relative errors turn out to
be rather small, while there is a very remarkable gain in computational times
dropping from hours down to few tens of seconds. We test our results by
varying the values of the correlation parameters and the maturity, moreover we
perform a brief sensitivity analysis with respect to the Hurst parameter, to
enucleate the roughness influence on the claim’s price.

The paper is structured as follows. In the next two sections, we first
introduce the theoretical framework for CVA evaluation and then we provide our
general representation and approximation formulas for a market model
characterized by a stochastic (either rough or not) volatility. Section 4
specializes the approximation formula to the previously mentioned models,
while section 5 is dedicated to a brief discussion of the error estimates. In
the final section, restricting to the more interesting rBergomi model, we
present the numerical analysis. We refer the reader to the appendix for a very
short primer on Malliavin Calculus.

##  2 CVA Evaluation of Defaultable European Claims

In this section, we briefly describe the general framework concerning
defaultable derivatives. Let $[0,T]$ be a finite time interval and
$(\Omega,\mathcal{F},P)$ a complete probability space endowed with a
filtration $\\{{\mathcal{F}}_{t}\\}_{t\in[0,T]}$, augmented with the $P-$null
sets and made right continuous. We also assume that all the processes have a
cádlág version.

The market is described by an interest rate $r_{t}$, and a process $X_{t}$,
representing an asset’s log-price, which may depend on multiple stochastic
factors. The filtration $\\{{\mathcal{F}}_{t}\\}_{\\{t\in[0,T]\\}}$ is rich
enough to support all those processes, and we assume to be in absence of
arbitrage, with $P$ a given risk neutral measure, selected by some criterion.
We denote the discount factor by $B(t,s)=\mathrm{e}^{-\int_{t}^{s}r_{u}du}$.

In this market, a defaultable European contingent claim paying $f(X_{T})$ at
maturity is traded, where $f$ is some function to be specified. During its
lifetime the claim is subject to default and we denote by $\tau$ (not
necessarily a stopping time w.r.t. the filtration ${\mathcal{F}}_{t}$) the
default time. At default, the contingent’s value may be partially recovered by
the creditor and $Z_{t}$ denotes an ${\mathcal{F}}_{t}-$measurable bounded
recovery process.

To properly evaluate this type of derivative, we need to include the
information generated by the default time. We denote by ${\mathcal{G}}_{t}$
the progressively enlarged filtration, that makes $\tau$ a
${\mathcal{G}}_{t}-$stopping time, that is
${\mathcal{G}}_{t}={\mathcal{F}}_{t}\lor\sigma(\\{\tau\leq t\\})$. Hence,
denoting by $H_{t}=\mathbf{1}_{\\{\tau\leq t\\}}$ and by ${\mathcal{H}}_{t}$
its natural filtration, we choose
${\mathcal{G}}_{t}={\mathcal{F}}_{t}\lor{\mathcal{H}}_{t}$.

We make the fundamental assumption, known as the H-hypothesis (see e.g. [19]
and [18] and the references therein), that

(H) Every ${\mathcal{F}}_{t}-$martingale remains a
${\mathcal{G}}_{t}-$martingale.

Under this assumption $B(t,s)S_{s},\,s\geq t$ remains a
${\mathcal{G}}_{s}-$martingale under the unique extension of the risk neutral
probability to the filtration ${\mathcal{G}}_{s}$. (To keep notation light, we
do not indicate the probability we use for the expectations, assuming to be
the one corresponding to the filtration in use).

In this setting, the price a defaultable claim, with final value $f(X_{T})$,
default time $\tau$ and recovery process $\\{Z_{t}\\}_{t}$ is given by

$c^{d}(t,T)=\mathbf{E}[B(t,T)f(X_{T})1_{\\{\tau>T\\}}+B(t,\tau)Z_{\tau}1_{\\{t<\tau\leq
T\\}}|{\mathcal{G}}_{t}],\quad t\in[0,T],$ (1)

while the corresponding default free value is
$c(t,T)=\mathbf{E}[B(t,T)f(X_{T})|{\mathcal{F}}_{t}].$

In many situations, investors do not know the default time and they may
observe only whether it happened or not. The actual observable quantity is the
asset price, therefore it is interesting to write the pricing formula (1) in
terms of ${\mathcal{F}}_{t}$, rather than in terms of ${\mathcal{G}}_{t}$. For
that we have the following key Lemma, see [12] or [10].

###### Lemma 2.1

For any integrable $\mathcal{G}-$measurable r.v. $Y$, the following equality
holds

$\mathbf{E}\Big{[}\mathbf{1}_{\\{\tau>t\\}}Y|{\mathcal{G}}_{t}\Big{]}=P(\tau>t|{\mathcal{G}}_{t})\frac{\mathbf{E}\Big{[}\mathbf{1}_{\\{\tau>t\\}}Y|{\mathcal{F}}_{t}\Big{]}}{P(\tau>t|{\mathcal{F}}_{t})}.$
(2)

Applying this lemma to the first and the second term of (1) and recalling that
$1-H_{t}=\mathbf{1}_{\\{\tau>t\\}}$ is ${\mathcal{G}}_{t}-$measurable, we
obtain

$\displaystyle\mathbf{E}[B(t,T)f(X_{T})1_{\\{\tau>T\\}}|{\mathcal{G}}_{t}]$
$\displaystyle=$
$\displaystyle\mathbf{1}_{\\{\tau>t\\}}\frac{\mathbf{E}[B(t,T)f(X_{T})1_{\\{\tau>T\\}}|{\mathcal{F}}_{t}]}{P(\tau>t|{\mathcal{F}}_{t})},$
(3) $\displaystyle\mathbf{E}[B(t,\tau)Z_{\tau}1_{\\{t<\tau\leq
T\\}}|{\mathcal{G}}_{t}]$ $\displaystyle=$
$\displaystyle\mathbf{1}_{\\{\tau>t\\}}\frac{\mathbf{E}[B(t,\tau)Z_{\tau}1_{\\{t<\tau\leq
T\\}}|{\mathcal{F}}_{t}]}{P(\tau>t|{\mathcal{F}}_{t})},$ (4)

which may be made more explicit by following the hazard process approach.

We denote the conditional distribution of the default time $\tau$ given
${\mathcal{F}}_{t}$ by

$\displaystyle F_{t}=P(\tau\leq t|{\mathcal{F}}_{t}),\qquad\forall\,t\geq 0,$
(5)

whence, for $u\geq t$, $P(\tau\leq u|{\mathcal{F}}_{t})=\mathbf{E}(P(\tau\leq
u|{\mathcal{F}}_{u})|{\mathcal{F}}_{t})=\mathbf{E}(F_{u}|{\mathcal{F}}_{t})$.
We also assume that $F_{t}(\omega)<1$ for all $t>0$ to well define the so-
called hazard process

$\Gamma_{t}:=-\ln(1-F_{t})\quad\Rightarrow\quad
F_{t}=1-\mathrm{e}^{-\Gamma_{t}}\quad\forall\,t>0,\qquad\Gamma_{0}=0.$ (6)

With this notation, by an extension of Proposition 5.1.1 of [11], we rewrite
(3) and (4) as

$\displaystyle\mathbf{E}[B(t,T)f(X_{T})1_{\\{\tau>T\\}}|{\mathcal{G}}_{t}]=$
$\displaystyle\mathbf{1}_{\\{\tau>t\\}}\mathbf{E}[B(t,T)f(X_{T})\mathrm{e}^{-(\Gamma_{T}-\Gamma_{t})}|{\mathcal{F}}_{t}],$
$\displaystyle\mathbf{E}[B(t,\tau)Z_{\tau}1_{\\{t<\tau\leq
T\\}}|{\mathcal{G}}_{t}]=$
$\displaystyle\mathbf{1}_{\\{\tau>t\\}}\mathbf{E}\Big{[}\int_{t}^{T}B(t,s)Z_{s}\mathrm{e}^{-(\Gamma_{s}-\Gamma_{t})}d\Gamma_{s}|{\mathcal{F}}_{t}\Big{]}.$

Assuming that the hazard process is differentiable with derivative
$\lambda_{t}$, called the intensity process, we arrive at the pricing formula

$\\!c^{d}(t,T)\\!=\\!\mathbf{1}_{\\{\tau>t\\}}\\!\left[\mathbf{E}\Big{(}\mathrm{e}^{-\\!\int_{t}^{T}(r_{s}\\!+\lambda_{s})ds}f(X_{T})\\!+\\!\\!\int_{t}^{T}\\!\\!\\!\\!Z_{s}\lambda_{s}\mathrm{e}^{-\\!\int_{t}^{s}(r_{u}\\!+\lambda_{u})du}ds|{\mathcal{F}}_{t}\Big{)}\right]\\!,$
(7)

recovering Lando’s formulas (3.1) and (3.3) in [26].

This formula can be specialized further, assuming fractional recovery ([16]),
$Z_{t}=Rc(t,T)$ for some $0\leq R<1$, and using the Optional Projection
Theorem (see [27] Theorem 4.16) to obtain

$\displaystyle
c^{d}(t,T)=\mathbf{1}_{\\{\tau>t\\}}\Big{[}R\mathbf{E}\Big{(}\mathrm{e}^{-\int_{t}^{T}r_{u}du}f(X_{T})|{\mathcal{F}}_{t}\Big{)}+(1-R)\mathbf{E}\Big{(}\mathrm{e}^{-\int_{t}^{T}(r_{u}+\lambda_{u})du}f(X_{T})|{\mathcal{F}}_{t}\Big{)}\Big{]},$
(8)

which was used also by Fard in [17] with $f(x)=(\mathrm{e}^{x}-K)^{+}$, and
that can be interpreted as a convex combination of the default free price and
the price with default.

As a consequence we have an expression also for the unilateral CVA, defined as
difference between the default free price and the adjusted price

$CVA(t):=\mathbf{1}_{\\{\tau>t\\}}[c(t,T)-c^{d}(t,T)]=\mathbf{1}_{\\{\tau>t\\}}(1-R)\mathbf{E}\Big{[}\mathrm{e}^{-\int_{t}^{T}r_{u}du}f(X_{T})(1-\mathrm{e}^{-\int_{t}^{T}\lambda_{u}du})|{\mathcal{F}}_{t}\Big{]}.$
(9)

Of course, the computability of these expectations will depend heavily on the
modeling choices one makes for $X$ and $\lambda$.

##  3 A representation formula for the CVA in stochastic volatility models

In this section, we consider a family of stochastic volatility models. For the
sake of simplicity, we assume zero fractional recovery ($R=0$), and the risk-
free spot rate, $r$, to be a deterministic function of time. These are not
restrictive assumptions, since the following discussion can be easily extended
to consider $0<R<1$, and increasing the dimensionality of the problem, also a
stochastic interest rate might be included. We assume that the probability
space contains at least a 3-dimensional standard ${\mathcal{F}}_{t}-$adapted
Brownian motion $\mathbf{B}_{s}=(B^{1}_{s},B^{2}_{s},B^{3}_{s})$, representing
the randomness sources of the market.

Indeed, our market model is described by an asset log-price process, $X_{t}$,
whose dynamics under a given risk neutral measure is

$dX_{t}=(r_{t}-\frac{1}{2}\sigma_{t}^{2})dt+\sigma_{t}dB_{t},\qquad
B_{t}=\sqrt{1-\eta^{2}}B^{1}_{t}+\eta B^{2}_{t},$ (10)

where the volatility process, $\sigma_{t}$, is assumed to be square-integrable
and adapted to the filtration generated by $\\{B^{2}_{t}\\}_{t\geq 0}$. The
parameter $\eta$ describes the correlation between the volatility and the log-
price.

We remark that the above model includes the classical volatility models, where
$\sigma$ is represented by a diffusion, but also those where the volatility’s
dynamics is determined by a fractional Brownian motion, in order to better
describe the long-term and short-term behavior of the implied volatility (see
[15], and [4]).

From now on, we take the shorter notation $\mathbf{E}_{t}[\cdot]$ to denote
the conditional expectation with respect to the filtration
${\mathcal{F}}_{t}$, and we consider the CVA problem given by formula (9), for
$R=0$ and $f(x)=(\mathrm{e}^{x}-\mathrm{e}^{\kappa})^{+}$, for some
$\kappa\in\mathbb{R}$. Thus we need to evaluate on $\\{\tau>t\\}$

$\mathrm{e}^{-\int_{t}^{T}r_{u}du}\mathbf{E}_{t}\Big{[}(1-\mathrm{e}^{-\int_{t}^{T}\lambda_{u}du})(\mathrm{e}^{X_{T}}-\mathrm{e}^{\kappa})^{+}\Big{]},$
(11)

which has to depend (as in the default free case) on the integrated mean
variance

$v_{t}=\frac{1}{T-t}\int_{t}^{T}\sigma_{u}^{2}du.$

For $0\leq t\leq s\leq T$, we denote by $\displaystyle
N^{t}_{T}=\mathrm{e}^{-\int_{t}^{T}\lambda_{u}du}$ and by

${\mathcal{B}}(s,x,\zeta)=\mathrm{e}^{x}N(d_{+}(s,x,\zeta))-\mathrm{e}^{\kappa}N(d_{-}(s,x,\zeta)),\,\,\text{with
}\,d_{\pm}(s,x,\zeta)\equiv
d_{\pm}=\frac{x-\kappa\pm\frac{\zeta^{2}}{2}(T-s)}{\zeta\sqrt{T-s}},$

the Black & Scholes European call pricing function, where $N$ is the standard
normal distribution function. We recall that ${\mathcal{B}}$ verifies

$\displaystyle\begin{cases}&\partial_{t}{\mathcal{B}}(s,x,\zeta)+\frac{\zeta^{2}}{2}(\partial_{xx}-\partial_{x}){\mathcal{B}}(s,x,\zeta)+r\partial_{x}{\mathcal{B}}(s,x,\zeta)-r{\mathcal{B}}(s,x,\zeta)=0\\\
&{\mathcal{B}}(T,x,\zeta)=(\mathrm{e}^{x}-\mathrm{e}^{\kappa})^{+}\end{cases},$
(12)
$\displaystyle\partial_{\zeta}{\mathcal{B}}(s,x,\zeta)=\zeta(T-s)(\partial_{xx}-\partial_{x}){\mathcal{B}}(s,x,\zeta),$
(13)
$\displaystyle\partial_{x}{\mathcal{B}}(s,x,\zeta)=\mathrm{e}^{x}N(d_{+}),\quad(\partial_{xx}-\partial_{x}){\mathcal{B}}(s,x,\zeta)=\frac{\mathrm{e}^{x}}{\zeta\sqrt{T-s}}N^{\prime}(d_{+}),$
(14)
$\displaystyle(\partial_{xxx}-\partial_{xx}){\mathcal{B}}(s,x,\zeta)=\frac{\mathrm{e}^{x}}{\sigma\sqrt{T-s}}N^{\prime}(d_{+})\left(1-\frac{d_{+}}{\zeta\sqrt{T-s}}\right),$
(15)
$\displaystyle(\partial_{xx}-\partial_{x})^{2}{\mathcal{B}}(s,x,\zeta)=\frac{\mathrm{e}^{x}N^{\prime}(d_{+})}{\zeta^{3}(T-s)^{\frac{3}{2}}}\left[d_{+}^{2}-d_{+}\zeta\sqrt{T-s}-1\right].$
(16)

Since the final condition
${\mathcal{B}}(T,x,\zeta)=(\mathrm{e}^{x}-\mathrm{e}^{\kappa})^{+}$ does not
depend on $\zeta$, we may rewrite the risk-neutral expectation (11) as

$\mathbf{E}_{t}\left[\mathrm{e}^{-\int_{t}^{T}r_{u}du}(1-N^{t}_{T}){\mathcal{B}}(T,X_{T},v_{T})\right]\mathbf{1}_{\\{\tau>t\\}}.$
(17)

###### Proposition 3.1

Assuming the market model (10), we have

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[\mathrm{e}^{-\int_{t}^{T}r_{u}du}(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]$
(18) $\displaystyle+$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\\!\\!\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]$ $\displaystyle-$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right],$

where by $D_{s}^{B^{i}}$ we denoted the Malliavin derivative with respect to
the $i-$th component of the Brownian motion $\mathbf{B}$.

Proof:  We denote $\displaystyle R_{s}^{-1}=\mathrm{e}^{-\int_{0}^{s}r_{u}du}$
and $\displaystyle Y_{t}=\int_{t}^{T}\sigma^{2}_{u}du$. We are going to apply
the anticipating Itô formula (see [5]) to the process $Y_{t}$ and the function

$F(s,x,y)=R_{s}^{-1}{\mathcal{B}}\left(s,x,\sqrt{\frac{y}{T-s}}\right).$

Exploiting relations (12) and (13), we have

$\displaystyle\partial_{s}F(s,x,y)=$ $\displaystyle
R_{s}^{-1}\Big{(}\partial_{s}{\mathcal{B}}+\frac{y}{2(T-s)}(\partial_{xx}-\partial_{x}){\mathcal{B}}-r_{s}{\mathcal{B}}\Big{)}{\mathcal{B}}\left(s,x,\sqrt{\frac{y}{T-s}}\right),$
$\displaystyle\partial_{x}F(s,x,y)=$ $\displaystyle
R_{s}^{-1}\partial_{x}{\mathcal{B}}\left(s,x,\sqrt{\frac{y}{T-s}}\right),\quad\partial_{xx}F(s,x,y)=R_{s}^{-1}{\mathcal{B}}\left(s,x,\sqrt{\frac{y}{T-s}}\right),$
$\displaystyle\partial_{y}F(s,x,y)=$
$\displaystyle\frac{1}{2}R_{s}^{-1}(\partial_{xx}-\partial_{x}){\mathcal{B}}\left(s,x,\sqrt{\frac{y}{T-s}}\right),$
$\displaystyle\partial_{xy}F(s,x,y)=$
$\displaystyle\frac{R_{s}^{-1}}{2}(\partial_{xxx}-\partial_{xx}){\mathcal{B}}\left(s,x,\sqrt{\frac{y}{T-s}}\right).$

Keeping in mind that $v_{t}=\sqrt{\frac{Y_{t}}{T-t}}$, we obtain

$\displaystyle CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[R_{[t,T]}^{-1}(1-N^{t}_{T}){\mathcal{B}}(T,X_{T},v_{T})\right]=\mathbf{E}_{t}\left[R_{[t,T]}^{-1}(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]$
$\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T})\int_{t}^{T}\frac{\eta}{2}R_{[t,s]}^{-1}\sigma_{s}(D^{-}Y)_{s}(\partial_{xxx}-\partial_{xx}){\mathcal{B}}(s,X_{s},v_{s})ds\right]$
$\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T})\int_{t}^{T}\sqrt{1-\eta^{2}}R_{[t,s]}^{-1}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})dB^{1}_{s}\right]$
$\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T})\int_{t}^{T}\eta
R_{[t,s]}^{-1}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})dB^{2}_{s}\right],$

where $R_{[t,s]}^{-1}=R_{s}^{-1}R_{t}$, and
$\displaystyle(D^{-}Y)_{s}=D^{B^{2}}_{s}\left(\int_{s}^{T}\sigma^{2}_{u}du\right)$.

By using the duality relation with respect to the two independent Brownian
motions for the last two integrals, we have

$\displaystyle CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[R_{[t,T]}^{-1}(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]$
$\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T})\int_{t}^{T}\frac{\eta}{2}R_{[t,s]}^{-1}\sigma_{s}(\partial_{xxx}-\partial_{xx}){\mathcal{B}}(s,X_{s},v_{s})\int_{s}^{T}D^{B^{2}}_{s}\sigma^{2}_{u}du\,ds\right]$
$\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[\int_{t}^{T}D^{B^{1}}_{s}(1-N^{t}_{T})R_{[t,s]}^{-1}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sqrt{1-\eta^{2}}ds\right]$
$\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[\int_{t}^{T}D^{B^{2}}_{s}(1-N^{t}_{T})R_{[t,s]}^{-1}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\eta
ds\right].$

Observing that

$\displaystyle D^{B^{i}}_{s}(1-N^{t}_{T})=$
$\displaystyle-N^{t}_{T}D^{B^{i}}_{s}\left(-\int_{t}^{T}\lambda_{u}du\right)=N^{t}_{T}\int_{s}^{T}D^{B^{i}}_{s}\lambda_{u}du,$
$\displaystyle\eta D^{B^{2}}_{s}\left(\int_{s}^{T}\sigma^{2}_{u}du\right)=$
$\displaystyle\eta\int_{s}^{T}D^{B^{2}}_{s}\sigma^{2}_{u}du=\sum_{i=1}^{2}\int_{s}^{T}D^{B^{i}}_{s}\sigma^{2}_{u}d\langle
B,B^{i}\rangle_{s},$

because the process $\lambda$ is ${\mathcal{F}}_{t}-$adapted and
$D^{B^{1}}\sigma^{2}_{s}=0$ ($\sigma$ depends only on $B^{2}$), we finally get
to

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[\mathrm{e}^{-\int_{t}^{T}r_{u}du}(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]$
$\displaystyle+$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\\!\\!\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]$ $\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right].\hskip 85.35826pt\square$

##  4 The approximation formula

The above representation formula is very compact, but clearly the expectations
cannot be evaluated in closed form and they need being approximated.

Besides, it is evident that the CVA is determined by the default intensity,
but also by the correlations among volatility, price and intensity, still
partially hidden in the factors given by ${\mathcal{B}}$ and its derivatives.
In what follows we construct an approximation that makes those contributions
explicit, so that it might be possible to fully parametrize them, in view of a
sensitivity analysis.

We therefore need to handle the dynamics of $\lambda$ and $\sigma$ less
implicitly and to approximate ${\mathcal{B}}$ and its derivatives, when
evaluated on the anticipating process $v_{t}$. To this aim, we introduce the
martingale $\displaystyle
N^{t}_{s}=\mathbf{E}_{s}\left(\mathrm{e}^{-\int_{t}^{T}\lambda_{u}du}\right)$
(with final value $N^{t}_{T}$), whose dynamics can be explicitly described
resorting to the bond pricing theory, when choosing the intensity $\lambda$ in
the class of affine processes. Indeed,

$N^{t}_{s}=\mathbf{E}_{s}\Big{(}\mathrm{e}^{-\int_{t}^{T}\lambda_{u}du}\Big{)}=\mathrm{e}^{-\int_{t}^{s}\lambda_{u}du}E_{s}\Big{(}\mathrm{e}^{-\int_{s}^{T}\lambda_{u}du}\Big{)}=\mathrm{e}^{-\int_{t}^{s}\lambda_{u}du}\mathrm{e}^{\varphi(T-s)\lambda_{s}+\psi(T-s)},$
(19)

for some deterministic differentiable functions $\varphi$ and $\psi$ of time
to maturity, which implies that

$\displaystyle dN^{t}_{s}=$
$\displaystyle-N^{t}_{s}\Big{\\{}[\lambda_{s}+\varphi^{\prime}(T-s)\lambda_{s}+\psi^{\prime}(T-s)]ds+\varphi(T-s)d\lambda_{s}-\frac{1}{2}\varphi^{2}(T-s)d\langle\lambda,\lambda\rangle_{s}\Big{\\}}$
$\displaystyle=$ $\displaystyle-N^{t}_{s}\varphi(T-s)dM^{\lambda}_{s},$

where $M^{\lambda}$ denotes the martingale part of the process $\lambda$.

As for the anticipating process $v_{t}$, we introduce the martingale

$M_{s}=\mathbf{E}_{s}\left(\int_{0}^{T}\sigma^{2}_{u}du\right),\quad\text{whence}\quad
v^{2}_{s}=\frac{1}{T-s}\left(M_{T}-\int_{0}^{s}\sigma^{2}_{u}du\right),$

while the so called zero-strike variance swap process is given by

$\hat{v}^{2}_{s}=\frac{1}{T-s}\mathbf{E}_{s}\left[\int_{s}^{T}\sigma_{u}^{2}du\right]=\frac{1}{T-s}\left(M_{s}-\int_{0}^{s}\sigma^{2}_{u}du\right)\quad\Rightarrow\quad\mathbf{E}_{s}(v^{2}_{s})=\hat{v}^{2}_{s}.$
(20)

Assuming that $\sigma\in\mathbb{D}^{1,2}$, we may apply the Clark-Ocone-
Haussman formula to deduce the dynamics of $M$

$M_{s}=\mathbf{E}\left[\int_{0}^{T}\\!\\!\\!\sigma^{2}_{u}du\right]+\\!\int_{0}^{s}\\!\\!\mathbf{E}_{s}\left[D^{B^{2}}_{s}\left(\int_{0}^{T}\\!\\!\\!\sigma^{2}_{u}du\right)\right]dB^{2}_{s}=\mathbf{E}\left[\int_{0}^{T}\\!\\!\\!\sigma^{2}_{u}du\right]+\\!\int_{0}^{s}\\!\\!\mathbf{E}_{s}\left[\int_{s}^{T}\\!\\!\\!D^{B^{2}}_{s}\sigma^{2}_{u}du\right]dB^{2}_{s}.$

We remark that this approach does not require the volatility to be a Markovian
process and it can be applied to a rough volatility, if one obtains a
manageable martingale representation. Indeed, in the next section we are going
to provide an explicit CVA approximation in the case of classical Heston and
SABR models, as well as for the Rough Bergomi volatility.

By this new notation, we may rewrite the representation formula (42) as

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[\mathrm{e}^{-\int_{t}^{T}r_{u}du}(1-N^{t}_{T}){\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{T}-\int_{0}^{t}\sigma^{2}_{u}du}{T-t}}\right)\right]$
$\displaystyle+$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\\!\\!\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]$ $\displaystyle+$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right],$

and we suggest a handy approximation formula for the second and third term. By
substituting $v^{2}_{s}$ with its conditional expectation $\hat{v}^{2}_{s}$
(so
$v_{s}=\sqrt{v^{2}_{s}}\approx\sqrt{\mathbf{E}_{s}(v^{2}_{s})}=\hat{v}_{s}$)
in the derivatives of ${\mathcal{B}}$, and freezing the state variables at
their initial values, those terms may be approximated by

$\begin{aligned}
&\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\\!\\!\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]\\\
&+\mathbf{E}_{t}\left[N^{t}_{T}\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right]\\\
\approx&\frac{1}{2}(1-N^{t}_{t})\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\left[\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sum_{i=1}^{2}\mathbf{E}_{s}\left(\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\right)\sigma_{s}d\langle
B,B^{i}\rangle_{s}\right]\\\
+&\partial_{x}{\mathcal{B}}(t,X_{t},v_{t})\mathbf{E}_{t}\left[\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sum_{i=1}^{2}\mathbf{E}_{s}\left(N^{t}_{T}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\right)\sigma_{s}d\langle
B,B^{i}\rangle_{s}\right]\\\
=&\frac{(1-N^{t}_{t})}{2}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\left[\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}d\langle
M,X\rangle_{s}\right]\\\
-&\partial_{x}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\left[\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}d\langle
N^{t},X\rangle_{s}\right],\end{aligned},$

where we used the Clark-Ocone formula to represent the martingales $N^{t}$ and
$M$

$dM_{s}=\sum_{i=1}^{3}\mathbf{E}_{s}\left(\int_{s}^{T}D^{B^{i}}_{s}\sigma^{2}_{u}du\right)dB^{i}_{s},\quad
dN^{t}_{s}=-\sum_{i=1}^{3}\mathbf{E}_{s}\left(N^{t}_{T}\int_{s}^{T}D^{B^{i}}_{s}\lambda_{u}du\right)dB^{i}_{s}$
(21)

and consequently their covariation processes with $X$. Notice that the higher-
order terms in this expansion include the quadratic variation and the
covariance of the martingales $N,X,M$. As these quantities depend on some
small parameters in most models (as we will see in the examples), these terms
will be higher-order as functions of these small parameters. For the first
term, we apply the classical Itô formula to the process $M$ and the function
$F(m)={\mathcal{B}}\left(t,x,\sqrt{\frac{m-\Sigma_{t}}{T-t}}\right)$, that, by
virtue of (13), has derivatives

$\displaystyle\partial_{m}F(m)=$
$\displaystyle\frac{1}{2}(\partial_{xx}-\partial_{x}){\mathcal{B}}\left(t,x,\sqrt{\frac{m-\Sigma_{t}}{T-t}}\right),$
$\displaystyle\partial_{mm}F(m)=$
$\displaystyle\frac{1}{4}(\partial_{xx}-\partial_{x})^{2}{\mathcal{B}}\left(t,x,\sqrt{\frac{m-\Sigma_{t}}{T-t}}\right).$

Recalling the expression of $\hat{v}_{s}$, the PDE (12), and expanding the
first term in (42) by integration by parts, we get

$\displaystyle\mathbf{E}_{t}\left[R_{[t,T]}^{-1}(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]=$
$\displaystyle(1-N^{t}_{t}){\mathcal{B}}(t,X_{t},\hat{v}_{t})$
$\displaystyle+$
$\displaystyle\frac{1}{8}\mathbf{E}_{t}\left[\int_{t}^{T}R_{[t,s]}^{-1}(1-N^{t}_{s})(\partial_{xx}-\partial_{x})^{2}{\mathcal{B}}(t,X_{t},\hat{v}_{s})d\langle
M,M\rangle_{s}\right]$ $\displaystyle-$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[\int_{t}^{T}R_{[t,s]}^{-1}(\partial_{xx}-\partial_{x}){\mathcal{B}}(t,X_{t},\hat{v}_{s})d\langle
N^{t},M\rangle_{s}\right],$

which we may approximate by freezing $N^{t}$ and $\hat{v}$ at the initial
time.

Summarizing, assuming model (10) and (20), the approximated value adjustment
of a vulnerable European call option is given by

$\boxed{\begin{aligned}
CVA(t)=&\mathbf{E}_{t}\left[\mathrm{e}^{-\int_{t}^{T}r_{u}du}(1-N^{t}_{T}){\mathcal{B}}(T,X_{T},v_{T})\right]\\\
\approx&(1\\!-\\!N^{t}_{t}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\\\
&+\frac{1}{8}(1\\!-\\!N^{t}_{t})(\partial_{xx}-\partial_{x})^{2}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\Big{[}\\!\int_{t}^{T}\\!\\!\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}d\langle
M\\!,M\rangle_{s}\Big{]}\\\
&+\frac{1}{2}(1\\!-\\!N^{t}_{t})(\partial_{xxx}-\partial_{xx}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\Big{[}\\!\int_{t}^{T}\\!\\!\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}d\langle
M\\!,X\rangle_{s}\Big{]}\\\
&-\frac{1}{2}(\partial_{xx}-\partial_{x}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\Big{[}\\!\int_{t}^{T}\\!\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}d\langle
N^{t}\\!,M\rangle_{s}\Big{]}\\\
&-\partial_{x}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\Big{[}\\!\int_{t}^{T}\\!\\!\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}d\langle
N^{t}\\!,X\rangle_{s}\Big{]}.\end{aligned}}$ (22)

###### Remark 4.1

In the above formula, we could leave the factor $(1-N^{t}_{s})$ inside the
integrals rather than freezing it also, at the initial time. Eventually, the
two choices give about the same numerical contribution to the final results,
thus we opted for a handier expression, to compute the various terms in a fast
and simple way.

Expression (22) might be interpreted as a first-order approximation formula,
indeed when restricting to some specific models for $\lambda$ and $\sigma$,
the final evaluation of the covariation processes will make the correlation
parameters appear with the first power. Because of its integrability
properties, in what follows, we take $\lambda$ as a CIR process verifying
Feller’s condition to ensure positivity, while for $\sigma$ we make some
different choices, trying to understand also the effect of roughness on the
evaluation. In the last section, we show numerically (comparing with Monte
Carlo simulations) how accurate this approximation can be.

To simplify our discussion, from now on we assume $r\equiv 0$.

##  5 The stochastic volatility-CIR intensity model

Here, we consider a CIR dynamics for the intensity process

$d\lambda_{s}=q(\mu-\lambda_{s})ds+c\sqrt{\lambda_{s}}dW_{s},\qquad\lambda_{t}=\lambda>0,\quad
s>t$ (23)

given that $c^{2}<2q\mu$, with $c,\mu,q>0$ and

$W_{s}=\frac{\rho-\eta\gamma}{\sqrt{1-\eta^{2}}}B^{1}_{s}+\gamma
B^{2}_{s}+\sqrt{\frac{1-(\eta^{2}+\gamma^{2}+\rho^{2})+2\gamma\eta\rho}{1-\eta^{2}}}B^{3}_{s},$
(24)

with parameters $(\gamma,\eta,\rho)$ such that

$\gamma^{2}<1,\,\,\,\eta^{2}<1,\,\,\rho^{2}<1,\quad\gamma^{2}+\rho^{2}+\eta^{2}<1+2\gamma\eta\rho,$
(25)

to obtain the correlations: $\langle B,B^{2}\rangle_{s}=\eta s,\langle
B,W\rangle_{s}=\rho s,\langle B^{2},W\rangle_{s}=\gamma s$.

By Fourier inversion, we know that the martingale $\displaystyle
N^{t}_{s}=\mathrm{e}^{-\int_{t}^{s}\lambda_{u}du}\mathrm{e}^{\varphi(T-s)\lambda_{s}+\psi(T-s)}$
has an explicit expression (see [25]) with

$\varphi(T-s)=-\frac{2(\mathrm{e}^{p(T-s)}-1)}{p-q+(p+q)\mathrm{e}^{p(T-s)}},\qquad\psi(T-s)=-\frac{2q\mu}{c^{2}}\ln\big{[}\frac{2p\mathrm{e}^{(p+q)(T-s)}}{p-q+(p+q)\mathrm{e}^{p(T-s)}}\Big{]},$

with $p^{2}=q^{2}+2c^{2}$ and it has dynamics
$dN^{t}_{s}=-c\varphi(T-s)N^{t}_{s}\sqrt{\lambda_{s}}dW_{s}$.

We are going to couple this intensity process with different choices for the
stochastic volatility.

###  5.1 Heston - CIR

The market model (see [23]) is given by

$\displaystyle dX_{s}=$
$\displaystyle-\frac{1}{2}\sigma^{2}_{s}ds+\sigma_{s}d(\sqrt{1-\eta^{2}}B^{1}_{s}+\eta
B^{2}_{s}),$ (26) $\displaystyle d\sigma^{2}_{s}=$ $\displaystyle
k(\theta-\sigma^{2}_{s})ds+\nu\sigma_{s}dB^{2}_{s},\quad\sigma_{t}=\sigma>0,$

and we assume that Feller’s condition is verified ($2k\theta>\nu^{2}$) also
for the volatility, to ensure the process positivity.

Consequently, for $r<u$, we find that the Malliavin derivative verifies

$D_{s}\sigma^{2}_{u}=\nu\sigma_{r}-\int_{s}^{u}kD_{s}\sigma^{2}_{\xi}d\xi+\nu\int_{s}^{u}\frac{1}{2\sigma_{\xi}}D_{s}\sigma^{2}_{\xi}dB^{2}_{\xi},$

whence taking expectation, and substituting the result in (21), we obtain

$\displaystyle\mathbf{E}_{s}(D_{s}\sigma^{2}_{u})=$
$\displaystyle\nu\sigma_{s}\mathrm{e}^{-k(u-s)}$ (27) $\displaystyle dM_{s}=$
$\displaystyle\frac{\nu}{k}(1-\mathrm{e}^{-k(T-s)})\sigma_{s}dB^{2}_{s}\qquad
M_{0}=\theta(1-\mathrm{e}^{-kT})$

(see [2]), so that we have

$\displaystyle d\langle M,X\rangle_{s}$ $\displaystyle=$
$\displaystyle\eta\nu\frac{1-\mathrm{e}^{-k(T-s)}}{k}\sigma_{s}^{2}ds,$
$\displaystyle d\langle M,M\rangle_{s}$ $\displaystyle=$
$\displaystyle\frac{\nu^{2}}{k^{2}}\left[1-\mathrm{e}^{-k(T-s)}\right]^{2}\sigma_{s}^{2}ds,$
$\displaystyle d\langle N^{t},X\rangle_{s}$ $\displaystyle=$
$\displaystyle-\rho c\varphi(T-s)N^{t}_{s}\sigma_{s}\sqrt{\lambda_{s}}ds,$
$\displaystyle d\langle N^{t},M\rangle_{s}$ $\displaystyle=$
$\displaystyle-\gamma\frac{\nu
c}{k}\varphi(T-s)\left[1-\mathrm{e}^{-k(T-s)}\right]N^{t}_{s}\sigma_{s},\sqrt{\lambda_{s}}ds$

which implies the approximation formula

$\displaystyle CVA(t)\approx$
$\displaystyle(1\\!-\\!N^{t}_{t}){\mathcal{B}}(t,X_{t},\hat{v}_{t})$ (28)
$\displaystyle+$
$\displaystyle\frac{1}{8}(1\\!-\\!N^{t}_{t})(\partial_{xx}-\partial_{x})^{2}{\mathcal{B}}(tX_{t},\hat{v}_{t})\frac{\nu^{2}}{k^{2}}\int_{t}^{T}\left[1-\mathrm{e}^{-k(T-s)}\right]^{2}\mathbf{E}_{t}\left[\sigma_{s}^{2}\right]ds$
$\displaystyle+$
$\displaystyle\frac{\eta}{2}(1\\!-\\!N^{t}_{t})(\partial_{xxx}-\partial_{xx}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\frac{\nu}{k}\int_{t}^{T}\left[1-\mathrm{e}^{-k(T-s}\right]\mathbf{E}_{t}\left[\sigma^{2}_{s}\right]ds$
$\displaystyle+$
$\displaystyle\frac{\gamma}{2}(\partial_{xx}-\partial_{x}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\frac{\nu
c}{k}\int_{t}^{T}\varphi(T-s)\left[1-\mathrm{e}^{-k(T-s)}\right]\mathbf{E}_{t}\left[N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\right]ds$
$\displaystyle+$ $\displaystyle\rho
c\partial_{x}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\int_{t}^{T}\varphi(T-s)\mathbf{E}_{t}\left[N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\right]ds.$

Recalling that
$\mathbf{E}_{t}\left[\sigma^{2}_{s}\right]=\theta+(\sigma_{t}^{2}-\theta)\mathrm{e}^{-k(s-t)}$,
we can compute exactly the first two integrals. For the others, we remark that
empirical experience from the market shows that the dependence between default
and stochastic volatility is rather weak, so that practitioners often consider
them as independent. We assume this viewpoint to approximate
$\mathbf{E}_{t}(N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s})$ by
$\mathbf{E}_{t}(N^{t}_{s}\sqrt{\lambda_{s}})\mathbf{E}_{t}(\sigma_{s})$, even
if $\gamma\neq 0$.

As in [7], to compute $\mathbf{E}_{t}(N^{t}_{s}\sqrt{\lambda_{s}})$ we use the
integration by parts formula, and we write

$d(\sqrt{\lambda_{s}}N^{t}_{s})=\frac{1}{2}\sqrt{\lambda_{s}}N^{t}_{s}\Big{[}\frac{4q\mu-c^{2}}{4\lambda_{s}}-(q+c^{2}\varphi(T\\!-\\!s))\Big{]}ds+cN^{t}_{s}\Big{[}\frac{1}{2}-\varphi(T-\\!s)\lambda_{s}\Big{]}dW_{s},$
(29)

whence, by considering that the martingale part gives null contribution, we
have

$\mathbf{E}_{t}(\sqrt{\lambda_{s}}N^{t}_{s})=\sqrt{\lambda}N^{t}_{t}+\int_{t}^{s}\frac{1}{2}\mathbf{E}_{t}\Big{(}\sqrt{\lambda_{u}}N^{t}_{u}\Big{[}\frac{4q\mu-c^{2}}{4\lambda_{u}}-(q+c^{2}\varphi(T-u))\Big{]}\Big{)}du.$

To approximate this last expectation, we freeze the $\frac{1}{\lambda_{s}}$
factor at the initial value and we solve the resulting ordinary differential
equation, obtaining

$\mathbf{E}_{t}(\sqrt{\lambda_{s}}N^{t}_{s})\approx\sqrt{\lambda}N^{t}_{t}\mathrm{e}^{\int_{t}^{s}g(u)du},\quad\textrm{
where}\quad
g(u)=\frac{4q\mu-c^{2}}{8\lambda}-\frac{1}{2}(q+c^{2}\varphi(T-u)).$ (30)

Finally, the factor $\mathbf{E}_{t}\Big{[}\sigma_{s}\Big{]}$ can be computed
by using a log-normal (moment-matching) approximation as in [7].

####  5.1.1 SABR - CIR

We now consider the SABR model ([22]). For $0\leq\beta\leq 1$ and
$\alpha\in\mathbb{R}$, we have

$\displaystyle dX_{s}$ $\displaystyle=$
$\displaystyle-\frac{\sigma_{s}^{2}}{2}\mathrm{e}^{-2(1-\beta)X_{s}}ds+\sigma_{s}\mathrm{e}^{-(1-\beta)X_{s}}d(\sqrt{1-\eta^{2}}B^{1}_{s}+\eta
B^{2}_{s}),$ (31) $\displaystyle d\sigma_{s}$ $\displaystyle=$
$\displaystyle\alpha\sigma_{s}dB^{2}_{s}.$ (32)

In this case, for any $r<u$

$\sigma_{u}=\sigma_{s}\mathrm{e}^{\alpha(B^{2}_{u}-B^{2}_{s})-\frac{\alpha^{2}}{2}(u-s)}\quad\Rightarrow\quad\sigma^{2}_{s}\mathrm{e}^{2\alpha(B^{2}_{u}-B^{2}_{s})-2\alpha^{2}(u-s)}$

whence

$\displaystyle D_{s}\sigma^{2}_{u}=$ $\displaystyle
2\alpha\sigma^{2}_{u}\quad\Rightarrow\quad\mathbf{E}_{s}(D_{s}\sigma^{2}_{u})=2\alpha\sigma^{2}_{s}\mathrm{e}^{\alpha^{2}(u-s)},\quad\text{and}$
$\displaystyle dM_{s}=$ $\displaystyle
2\alpha\sigma_{s}^{2}\int_{s}^{T}\mathrm{e}^{\alpha^{2}(u-s)}dudB^{2}_{s}=2\sigma_{s}^{2}\frac{\mathrm{e}^{\alpha^{2}(T-s)}-1}{\alpha}dB^{2}_{s},$

and

$\displaystyle d\langle M,X\rangle_{s}$
$\displaystyle=\eta\sigma_{s}^{3}\mathrm{e}^{-(1-\beta)X_{s}}\frac{\mathrm{e}^{\alpha^{2}(T-s)}-1}{\alpha}ds,$
(33) $\displaystyle d\langle M,M\rangle_{s}$
$\displaystyle=4\sigma_{s}^{4}\frac{(\mathrm{e}^{2\alpha^{2}(T-s)}-1)^{2}}{\alpha^{2}}ds,$
$\displaystyle d\langle N^{t},X\rangle_{s}$ $\displaystyle=-\rho
c\varphi(T-s)N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\mathrm{e}^{-(1-\beta)X_{s}}ds,$
$\displaystyle d\langle N^{t},M\rangle_{s}$ $\displaystyle=-2\gamma
c\varphi(T-s)\frac{\mathrm{e}^{\alpha^{2}(T-s)}-1}{\alpha}N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}^{2}ds,$

so that (22) may be approximated by

$\displaystyle CVA(t)\approx$
$\displaystyle(1\\!-\\!N^{t}_{t}){\mathcal{B}}(t,X_{t},\hat{v}_{t})$ (34)
$\displaystyle+$
$\displaystyle\frac{1}{2}(1\\!-\\!N^{t}_{t})(\partial_{xx}-\partial_{x})^{2}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\int_{t}^{T}\\!\\!\frac{(\mathrm{e}^{2\alpha^{2}(T-s)}-1)^{2}}{\alpha^{2}}\mathbf{E}_{t}\big{[}\sigma_{s}^{4}\big{]}ds$
$\displaystyle+$
$\displaystyle\frac{\eta}{2}(1\\!-\\!N^{t}_{t})(\partial_{xxx}\\!\\!-\partial_{xx}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\\!\int_{t}^{T}\\!\\!\frac{\mathrm{e}^{\alpha^{2}(T\\!-s)}\\!-1}{\alpha}E_{t}\Big{[}\sigma_{s}^{3}\mathrm{e}^{-(1\\!-\beta)X_{s}}\Big{]}ds$
$\displaystyle+$ $\displaystyle\gamma
c(\partial_{xx}-\partial_{x}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\Big{[}\\!\int_{t}^{T}\\!\\!\\!\varphi(T-s)\frac{\mathrm{e}^{\alpha^{2}(T-s)}-1}{\alpha}N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}^{2}ds\Big{]}$
$\displaystyle+$ $\displaystyle\rho
c\partial_{x}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\mathbf{E}_{t}\Big{[}\\!\int_{t}^{T}\\!\\!\\!\varphi(T-s)N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\mathrm{e}^{-(1-\beta)X_{s}}ds\Big{]}.$

To get to an actually implementable formula,

  * •

we freeze the term $\mathrm{e}^{-(1-\beta)X_{s}}$ at the initial condition,
obtaining $\mathrm{e}^{-(1-\beta)X_{t}}$;

  * •

for $n=1,2$, as before we treat the terms
$\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}^{n}\Big{]}$ as
regarding independent processes, obtaining
$\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\Big{]}\mathbf{E}_{t}\big{[}\sigma_{s}^{n}\big{]}$;

  * •

we approximate $\mathbf{E}_{t}(N^{t}_{s}\sqrt{\lambda_{s}})$ as in (30);

  * •

for $n=1,2,3,4$, we finally recall that
$\mathbf{E}_{t}(\sigma^{n}_{s})=\sigma_{t}^{n}\mathrm{e}^{n(n-1)\frac{\alpha^{2}}{2}(s-t)}$.

####  5.1.2 rBergomi - CIR

In this subsection, we analyze our formula in the case of the Rough Bergomi
volatility model.

In this case, the market model is still (10), where $\sigma_{s}$ follows

$\sigma_{s}^{2}=\sigma_{0}^{2}\mathrm{e}^{\nu\sqrt{2H}Z_{s}-\frac{\nu^{2}}{2}s^{2H}},\quad\nu,\sigma_{0}>0,\quad
H<\frac{1}{2},$ (35)

with $Z$ a RLfBm of the form $\displaystyle
Z_{s}:=\int_{0}^{s}(s-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}.$

Consequently

$\displaystyle D_{s}\sigma_{u}^{2}=$
$\displaystyle\nu\sqrt{2H}\sigma_{u}^{2}D_{s}\left[\int_{0}^{u}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}-\frac{1}{2}\nu^{2}u^{2H}\right]=\nu\sqrt{2H}\sigma_{u}^{2}D_{s}\left[\int_{s}^{u}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}\right]$
(36) $\displaystyle=$
$\displaystyle\nu\sqrt{2H}\sigma_{u}^{2}(u-s)^{H-\frac{1}{2}},$
$\displaystyle\Rightarrow$
$\displaystyle\qquad\mathbf{E}_{s}(D_{s}\sigma_{u}^{2})=\nu\sqrt{2H}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{s}(\sigma_{u}^{2}).$

It follows that

$dM_{s}=\nu\sqrt{2H}\int_{s}^{T}(u-s)^{H-\frac{1}{2}}E_{s}\big{(}\sigma_{u}^{2}\big{)}du\,dB^{2}_{s},$
(37)

and therefore

$\displaystyle d\langle M,X\rangle_{s}$
$\displaystyle=\eta\nu\sqrt{2H}\sigma_{s}\int_{s}^{T}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}duds,$
(38) $\displaystyle d\langle M,M\rangle_{s}$
$\displaystyle=\nu^{2}2H\left[\int_{s}^{T}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}du\right]^{2}ds,$
$\displaystyle d\langle N^{t},X\rangle_{s}$ $\displaystyle=-\rho
c\varphi(T-s)N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}ds,$ $\displaystyle d\langle
N^{t},M\rangle_{s}$ $\displaystyle=-\gamma
c\nu\sqrt{2H}\varphi(T-s)N^{t}_{s}\sqrt{\lambda_{s}}\int_{s}^{T}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}duds.$

When substituting in the approximation formula (22), using the stochastic
Fubini theorem, we finally obtain

$\displaystyle CVA(t)\approx$
$\displaystyle(1\\!-\\!N^{t}_{t}){\mathcal{B}}(t,X_{t},\hat{v}_{t})$ (39)
$\displaystyle+$
$\displaystyle\frac{\nu^{2}H}{4}(1\\!-\\!N^{t}_{t})(\partial_{xx}\\!\\!-\partial_{x})^{2}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\\!\int_{t}^{T}\\!\\!\\!\mathbf{E}_{t}\left[\Big{[}\int_{s}^{T}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}du\Big{]}^{2}\right]ds$
$\displaystyle+$
$\displaystyle\frac{\eta}{2}\nu\sqrt{2H}(1\\!-\\!N^{t}_{t})(\partial_{xxx}\\!\\!-\partial_{xx}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\\!\int_{t}^{T}\\!\\!\\!\int_{s}^{T}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{t}\Big{[}\sigma_{s}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}\Big{]}duds$
$\displaystyle+$
$\displaystyle\rho\frac{c}{2}(\partial_{xx}\\!\\!-\partial_{x}){\mathcal{B}}(t,X_{t},\hat{v}_{t})\\!\int_{t}^{T}\\!\\!\\!\varphi(T-s)\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\Big{]}ds$
$\displaystyle+$ $\displaystyle\gamma
c\nu\sqrt{2H}\partial_{x}{\mathcal{B}}(t,X_{t},\hat{v}_{t})\int_{t}^{T}\\!\\!\\!\varphi(T-s)\int_{s}^{T}(u-s)^{H-\frac{1}{2}}\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}\Big{]}duds.$

To make the above expression more manageable, further computations and
approximations are required.

  * •

We first remark that

$\displaystyle\mathbf{E}_{t}[\sigma^{2}_{u}]=$
$\displaystyle\mathbf{E}_{t}\big{(}\sigma_{0}^{2}\mathrm{e}^{\nu\sqrt{2H}Z_{u}-\frac{\nu^{2}}{2}u^{2H}}\big{)}=\sigma_{0}^{2}\mathrm{e}^{-\frac{\nu^{2}}{2}u^{2H}}\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{t}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}}E_{t}\big{(}\mathrm{e}^{\nu\sqrt{2H}\int_{t}^{u}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}}\big{)}$
$\displaystyle=$
$\displaystyle\sigma_{0}^{2}\mathrm{e}^{\frac{\nu^{2}}{2}[(u-t)^{2H}-u^{2H}]}\,\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{t}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}},$

since
$\displaystyle\sqrt{2H}\int_{t}^{u}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}\sim
N\big{(}0,(u-t)^{2H}\big{)}$. In the second term we have $t\leq s\leq u$, so
$\mathbf{E}_{t}(\mathbf{E}_{s}(\sigma^{2}_{u}))=\mathbf{E}_{t}(\sigma^{2}_{u})$.

  * •

It follows that in the first integral term we have

$\displaystyle\mathbf{E}_{t}\big{(}\sigma_{s}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}\big{)}=$
$\displaystyle\mathbf{E}_{t}\Big{(}\sigma_{s}\sigma_{0}^{2}\mathrm{e}^{\frac{\nu^{2}}{2}[(u-s)^{2H}-u^{2H}]}\,\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{s}(u-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}}\Big{)}$
$\displaystyle=$
$\displaystyle\mathbf{E}_{t}\Big{(}\sigma_{0}^{3}\mathrm{e}^{\frac{\nu^{2}}{2}[(u-s)^{2H}-u^{2H}-\frac{1}{2}s^{2H}]}\,\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{s}[(u-\xi)^{H-\frac{1}{2}}+\frac{1}{2}(s-\xi)^{H-\frac{1}{2}}]dB^{2}_{\xi}}\Big{)}$
$\displaystyle=$
$\displaystyle\sigma_{0}^{3}\mathrm{e}^{\frac{\nu^{2}}{2}[(u-s)^{2H}-u^{2H}-\frac{1}{2}s^{2H}]}\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{t}[(u-\xi)^{H-\frac{1}{2}}+\frac{1}{2}(s-\xi)^{H-\frac{1}{2}}]dB^{2}_{\xi}}$
$\displaystyle\times$
$\displaystyle\mathbf{E}_{t}\Big{(}\mathrm{e}^{\nu\sqrt{2H}\int_{t}^{s}[(u-\xi)^{H-\frac{1}{2}}+\frac{1}{2}(s-\xi)^{H-\frac{1}{2}}]dB^{2}_{\xi}}\Big{)}$
$\displaystyle=$
$\displaystyle\sigma_{0}^{3}\mathrm{e}^{\frac{\nu^{2}}{2}[(u-s)^{2H}-u^{2H}-\frac{1}{2}s^{2H}]}\mathrm{e}^{\nu^{2}\sqrt{2H}\int_{0}^{t}[(u-\xi)^{H-\frac{1}{2}}+\frac{1}{2}(s-\xi)^{H-\frac{1}{2}}]dB^{2}_{\xi}}\mathrm{e}^{\nu^{2}H\Theta^{2}_{s}},$

where

$\displaystyle\Theta^{2}_{s}=$
$\displaystyle\int_{t}^{s}\big{[}(u-\xi)^{H-\frac{1}{2}}+\frac{1}{2}(s-\xi)^{H-\frac{1}{2}}\big{]}^{2}d\xi=\frac{(s-t)^{2H}}{4H}+\left(\frac{(u-t)^{2H}-(u-s)^{2H}}{2H}\right)$
(40) $\displaystyle+$ $\displaystyle
2\frac{(s-t)^{H+\frac{1}{2}}(u-t)^{H-\frac{1}{2}}}{2H+1}\left(\frac{t-u}{s-u}\right)^{\frac{1}{2}-H}\\!\\!\\!\,_{2}F_{1}\Big{(}\frac{1}{2}-H,H+\frac{1}{2};H+\frac{3}{2};\frac{s-t}{s-u}\Big{)},$

with $\,{}_{2}F_{1}$ the hypergeometric function555Symbolic computation by
means of Wolfram Mathematica 12.2. and numerically tested on MatLab R2019b..

  * •

We observe that the expectation in the second integral term can be written as

$\displaystyle\mathbf{E}_{t}\left[\Big{[}\int_{s}^{T}\\!\\!\\!(u\\!-s)^{H-\frac{1}{2}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}du\Big{]}^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbf{E}_{t}\left[\int_{s}^{T}\\!\\!(u\\!-\\!s)^{H-\frac{1}{2}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}du\\!\\!\int_{s}^{T}\\!\\!\\!(\vartheta\\!-\\!s)^{H-\frac{1}{2}}\mathbf{E}_{s}(\sigma_{\vartheta}^{2})d\vartheta\right]$
$\displaystyle=$
$\displaystyle\int_{s}^{T}\\!\\!\\!\int_{s}^{T}\\!\\!\sigma_{0}^{4}\mathrm{e}^{\frac{\nu^{2}}{2}[(u-s)^{2H}-u^{2H}+(\vartheta\\!-\\!s)^{2H}-\vartheta^{2H}]}(u\\!-\\!s)^{H-\frac{1}{2}}(\vartheta\\!-\\!s)^{H-\frac{1}{2}}$
$\displaystyle\qquad\times\mathbf{E}_{t}\Big{[}\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{s}[(u-\xi)^{H-1/2}+(\vartheta-\xi)^{H-1/2}]dB^{2}_{\xi}}\Big{]}dud\vartheta$
$\displaystyle=$
$\displaystyle\sigma_{0}^{4}\mathrm{e}^{\nu\sqrt{2H}\int_{0}^{t}[(u-\xi)^{H-1/2}+(\vartheta-\xi)^{H-1/2}]dB^{2}_{\xi}}$
$\displaystyle\times\int_{s}^{T}\\!\\!\int_{s}^{T}\\!\\!\mathrm{e}^{\frac{\nu^{2}}{2}[(u-s)^{2H}-u^{2H}+(\vartheta\\!-\\!s)^{2H}-\vartheta^{2H}]}(u\\!-\\!s)^{H-\frac{1}{2}}(\vartheta\\!-\\!s)^{H-\frac{1}{2}}\mathrm{e}^{\nu^{2}H\Phi^{2}_{s}}dud\vartheta,$

where we used again the expression of the volatility second moment and we
set666Symbolic computation by means of Wolfram Mathematica 12.2. and
numerically tested on MatLab R2019b.

$\displaystyle\Phi^{2}_{s}=$
$\displaystyle\int_{t}^{s}[(u\\!-\\!\xi)^{H-\frac{1}{2}}+(\vartheta\\!-\\!\xi)^{H-\frac{1}{2}}]^{2}d\xi=\frac{(u\\!-\\!t)^{2H}-(u\\!-\\!s)^{2H}+(\vartheta\\!-\\!t)^{2H}-(\vartheta\\!-\\!s)^{2H}}{2H}$
(41) $\displaystyle+$ $\displaystyle
2\frac{\Gamma(H+\frac{1}{2})}{\Gamma(H+\frac{3}{2})(u-\vartheta)}\Big{[}(u-s)^{H+\frac{1}{2}}(\vartheta-s)^{H+\frac{1}{2}}{}_{2}F_{1}\Big{(}1,1+2H;H+\frac{3}{2};\frac{u-s}{u-\vartheta}\Big{)}$
$\displaystyle-$
$\displaystyle(u-t)^{H+\frac{1}{2}}(\vartheta-t)^{H+\frac{1}{2}}\\!\,_{2}F_{1}\Big{(}1,1+2H;H+\frac{3}{2};\frac{u-t}{u-\vartheta}\Big{)}\Big{]}.$

  * •

It remains to compute the two expectations

$\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\Big{]},\qquad\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}\Big{]}.$

As before, we consider the volatility and the intensity processes as
independent, hence

$\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\sigma_{s}\Big{]}=\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\Big{]}\mathbf{E}_{t}\big{[}\sigma_{s}\big{]},\quad\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\mathbf{E}_{s}\big{(}\sigma_{u}^{2}\big{)}\Big{]}=\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\Big{]}\mathbf{E}_{t}\big{[}\sigma_{u}^{2}\big{]}.$

We approximate $\mathbf{E}_{t}\Big{[}N^{t}_{s}\sqrt{\lambda_{s}}\Big{]}$ as in
(30), $\mathbf{E}_{t}\big{[}\sigma_{u}^{2}\big{]}$ was computed before, and
similarly

$\displaystyle\mathbf{E}_{t}\big{[}\sigma_{s}\big{]}=$
$\displaystyle\sigma_{0}\mathrm{e}^{\frac{\nu\sqrt{H}}{\sqrt{2}}\int_{0}^{t}(s-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}-\frac{\nu^{2}}{4}s^{2H}}\mathbf{E}_{t}\Big{[}\mathrm{e}^{\frac{\nu\sqrt{H}}{\sqrt{2}}\int_{t}^{s}(s-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}}\Big{]}$
$\displaystyle=$
$\displaystyle\sigma_{0}\mathrm{e}^{\frac{\nu\sqrt{H}}{\sqrt{2}}\int_{t}^{s}(s-\xi)^{H-\frac{1}{2}}dB^{2}_{\xi}+\frac{\nu^{2}}{4}[(s-t)^{2H}-s^{2H}]}$

###### Remark 5.1

By grouping all the terms, we finally get an approximation of the form

$\widehat{CVA}(t)=(1-N_{t}^{t})c^{rf}_{rBergomi}(t,X_{t},\hat{v}_{t})+\rho
g_{1}(t)+\gamma g_{2}(t),$

where $1-N_{t}^{t}$ is the survival probability of the counterparty,
$c^{rf}_{rBergomi}$ is the price of the risk-free call under the rBergomi
model and $g_{1,2}(t)$ are the respective contributions due to underlying-
intensity and volatility-intensity correlations.

###### Remark 5.2

From the previous computations, it appears clear that, by the Markov property,
in the first two cases $CVA(t)$ is a deterministic function of the state
variables, while the same is not verified in the case of the rough volatility
models, indeed the conditional expectations, when computed, leave a stochastic
integral term.

##  6 Error Estimates

In this section, we want to estimate the error bounds theoretically.We take
$r\equiv 0$ for the sake of exposition. Starting from representation We know
from Proposition 3.1 that

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]$
(42) $\displaystyle+$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\\!\\!\int_{t}^{T}\\!\sigma_{s}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]$ $\displaystyle-$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\int_{t}^{T}\\!\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right]$ $\displaystyle=$ $\displaystyle T_{1}+T_{2}+T_{3}.$

Applying again the anticipating Itô formula to the process

$(1-N^{t}_{T})\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(\theta,X_{\theta},v_{\theta})\\!\\!\int_{\theta}^{T}\\!\\!\sigma_{s}\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}$

after taking expectations, keeping in mind that this process is $0$ when
$\theta=T$, we obtain

$\displaystyle 0=$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T})\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(t,X_{t},v_{t})\\!\\!\int_{t}^{T}\\!\\!\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!\sigma_{s}D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right.$ (43) $\displaystyle-$
$\displaystyle\left.(1-N^{t}_{T})\\!\\!\int_{t}^{T}\\!\\!\mathrm{e}^{-\int_{t}^{s}r_{u}du}\sigma_{s}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]+T_{4},$

where $T_{4}$ is a term depending on the product of the Malliavin derivatives
of $\lambda$ and $\sigma$ and higher-order Malliavin derivatives. Similarly,
for the third term we can see that

$\displaystyle 0=$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\int_{t}^{T}\\!\sigma_{s}\partial_{x}{\mathcal{B}}(s,X_{s},v_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right.$ (44) $\displaystyle-$
$\displaystyle\left.N^{t}_{T}\partial_{x}{\mathcal{B}}(t,X_{t},v_{t})\int_{t}^{T}\\!\\!\sigma_{s}\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right]+T_{5},$

where $T_{5}$ is again a term depending on the product of the Malliavin
derivatives of $\lambda$ and $\sigma$ and higher-order Malliavin derivatives.
Substituting in (42), we have

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]$
(45) $\displaystyle+$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(t,X_{t},v_{t})\\!\\!\int_{t}^{T}\\!\\!\sigma_{s}\sum_{i=1}^{2}\\!\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\sigma_{u}^{2}du\,d\langle
B,B^{i}\rangle_{s}\right]$ $\displaystyle-$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\partial_{x}{\mathcal{B}}(t,X_{t},v_{t})\int_{t}^{T}\\!\\!\sigma_{s})\sum_{i=1}^{2}\int_{s}^{T}\\!\\!\\!D_{s}^{B^{i}}\\!\lambda_{u}du\,d\langle
B,B^{i}\rangle_{s}\right]+T_{4}+T_{5}$ $\displaystyle=$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T}){\mathcal{B}}(t,X_{t},v_{t})\right]+\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}(t,X_{t},v_{t})\\!\\!\int_{t}^{T}d\langle
X,M\rangle\right]$ $\displaystyle-$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\partial_{x}{\mathcal{B}}(t,X_{t},v_{t})\int_{t}^{T}\\!\\!d\langle
X,N\rangle_{s}\right]+T_{4}+T_{5}$ $\displaystyle=$ $\displaystyle
T_{1}^{\prime}+T_{2}^{\prime}+T_{3}^{\prime}+T_{4}+T_{5}.$

Denoting by $\displaystyle
A_{r}^{1}:=E_{r}\left[\int_{t}^{T}\\!\\!\\!\\!d\langle
M\\!,X\rangle_{s}\right],A_{r}^{2}:=E_{r}\left[\int_{t}^{T}\\!\\!\\!\\!d\langle
N\\!,X\rangle_{s}\right]$ and recalling the expression of $v$, the above
equality reads as

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle\mathbf{E}_{t}\left[(1-N^{t}_{T}){\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{T}-\int_{0}^{t}\sigma_{u}^{2}du}{T-t}}\right)\right]$
(46) $\displaystyle+$
$\displaystyle\frac{1}{2}\mathbf{E}_{t}\left[(1-N^{t}_{T})\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{T}-\int_{0}^{t}\sigma_{u}^{2}du}{T-t}}\right)\\!\\!A_{T}^{1}\right]$
$\displaystyle-$
$\displaystyle\mathbf{E}_{t}\left[N^{t}_{T}\partial_{x}{\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{T}-\int_{0}^{t}\sigma_{u}^{2}du}{T-t}}\right)A_{T}^{2}\right]+T_{4}+T_{5}.$

The first three terms in the above equality are functions of the martingales
$M,N,A^{1}$ and $A^{2}$ evaluated at time $T$. By a direct application of
Itô’s formula, we finally rewrite (46) as

$\displaystyle\\!\\!\\!CVA(t)=$
$\displaystyle(1\\!-\\!N^{t}_{t})\left[{\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{t}\\!-\\!\int_{0}^{t}\sigma_{u}^{2}du}{T-t}}\right)+\frac{1}{2}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{t}\\!-\\!\int_{0}^{t}\sigma_{u}^{2}du}{T-t}}\right)\\!\\!A_{t}^{1}\right]$
(47) $\displaystyle-$ $\displaystyle
N^{t}_{t}\partial_{x}{\mathcal{B}}\left(t,X_{t},\sqrt{\frac{M_{t}-\int_{0}^{t}\sigma_{u}^{2}du}{T-t}}\right)A_{t}^{2}+T_{1}^{\prime\prime}+T_{2}^{\prime\prime}+T_{3}^{\prime\prime}+T_{4}+T_{5}$
$\displaystyle=$
$\displaystyle(1-N^{t}_{t})\left[{\mathcal{B}}\left(t,X_{t},\hat{v}_{t}\right)+\frac{1}{2}\left(\partial_{xxx}\\!-\\!\partial_{xx}\right){\mathcal{B}}\left(t,X_{t},\hat{v}_{t}\right)\\!\\!A_{t}^{1}\right]-N^{t}_{t}\partial_{x}{\mathcal{B}}\left(t,X_{t},\hat{v}_{t}\right)A_{t}^{2}$
$\displaystyle+$ $\displaystyle
T_{1}^{\prime\prime}+T_{2}^{\prime\prime}+T_{3}^{\prime\prime}+T_{4}+T_{5}.$

All the terms $T_{1}^{\prime\prime}+T_{2}^{\prime\prime}+T_{3}^{\prime\prime}$
can be expressed in terms of the quadratic variation and covariation of the
martingales $N,M,A^{1}$, and $A^{2}$, which, in turn, may be expressed in
terms of the Malliavin derivatives of $\sigma$ and $\lambda$ by using the
Clark-Ocone-Haussman formula.

We may conclude that the error is controlled by the well known bounds on the
Black & Scholes pricing function and its derivatives, while those Malliavin
derivatives, which become explicitly computable once the models are specified.

##  7 Numerical results

In this section, we compare the numerical efficiency of our method with the
standard benchmark Monte Carlo technique to evaluate (9) in the rBergomi
model. As said, we assume $r=0$ and $R=0$. We simulated $N=10^{6}$ sample
paths of the joint model rBergomi-CIR with a discretization of the time
interval $[0,T]$ obtained by $100$ equi-spaced points. The exact covariance
matrix of the three gaussian processes $(B,Z,W)$ has been used to get the
sample paths (see [5]). Hence, the rBergomi is simulated by standard
multivariate normal generation, while for the intensity component we
implemented the Euler discretization scheme with full truncation. We chose two
sets of parameters for the intensity CIR process as in [13], set A and set B
(see Table 1), which agree with those of calibrated default intensities. Then
we fixed the rBergomi model parameters to $\nu=0.1$, $H=0.1$,
$\sigma_{0}=0.08$ and the underlying-volatility correlation to $\eta=-0.2$,
and we considered a set of equi-spaced underlying-intensity and volatility-
intensity correlations, $\rho$ and $\gamma$ respectively. In fact we would
like to assess the quality of the our approximation which highlights the
dependence of the credit adjustment on the default intensity correlations,
known as right/wrong-way-risk. All the computational procedures were
implemented by using MatLab R2021b on a Intel(R) Core(TM) i7-4700MQ CPU @
2.40GHz. The results of the approximations are reported in Tables (2)-(7) and
shown in the Figures (1)-(6) for an ATM call option with $S_{0}=100$, $K=100$
and maturities $T=0.25,0.5,1$ years, respectively. The time to compute the set
of $45$ adjustments with the Monte Carlo algorithm was approximately one hour,
while about $30$ secs are needed for the computation of our approximation,
which requires the numerical evaluation of multiple integrals. This step was
done by using the MatLab build-in function. We may notice that the
approximation errors increase with maturity, ranging from $1.18e-04$ to
$1.7e-03$ for set A and from $1.19e-06$ to $4.6e-03$ for set B. This was to be
expected since in our approximation we freeze some terms at their initial
value and therefore their contribution in time dimension is lost. We notice
that the error due to the volatility-intensity correlation remains essentially
constant across all simulations, while most of the variation is due to the
underlying-intensity correlation, which proves to be the most relevant
parameter. We notice the the errors pattern is different between the two
intensity sets of parameters: for set A the larger errors come from small
(underlying-intensity) correlations, almost the opposite for set B.

Finally, we considered also the sensitivity of our approximation to the Hurst
parameter $H$. Figure (7) reports the errors between Monte Carlo estimate and
the first order approximation as a function of the correlation $\rho$ and the
Hurst parameter $H$. Even in this case, the dependence of the error from $H$
is almost constant and depends mainly on the correlation parameter, however
remaning of order $10^{-4}$.

| $\lambda_{0}$ | $q$ | $\mu$ | $c$
---|---|---|---|---
Set A | 0.035 | 0.35 | 0.035 | 0.1
Set B | 0.01 | 0.8 | 0.02 | 0.2
Table 1:  Parameter sets for the CIR default intensity.

###  7.1 Results - Set A

$\rho\backslash\gamma$ | -0.3 | -0.15 | 0 | 0.15 | 0.3
---|---|---|---|---|---
-0.8 | 1.186e-04 | 1.195e-04 | 1.199e-04 | 1.206e-04 | 1.214e-04
-0.6 | 1.501e-04 | 1.504e-04 | 1.508e-04 | 1.513e-04 | 1.520e-04
-0.4 | 1.748e-04 | 1.750e-04 | 1.753e-04 | 1.758e-04 | 1.764e-04
-0.2 | 1.933e-04 | 1.934e-04 | 1.937e-04 | 1.940e-04 | 1.945e-04
0.0 | 2.056e-04 | 2.056e-04 | 2.057e-04 | 2.060e-04 | 2.064e-04
0.2 | 2.116e-04 | 2.116e-04 | 2.116e-04 | 2.118e-04 | 2.121e-04
0.4 | 2.115e-04 | 2.113e-04 | 2.112e-04 | 2.114e-04 | 2.116e-04
0.6 | 2.051e-04 | 2.048e-04 | 2.047e-04 | 2.048e-04 | 2.050e-04
0.8 | 1.926e-04 | 1.922e-04 | 1.920e-04 | 1.921e-04 | 1.920e-04
Table 2: Approximation errors, set A for CIR parameters. Here $T=1/4$, the risk-free call price is $1.5877$. The lengths of the Monte Carlo confidence intervals for each price are in the range $(6.7651\mathrm{e}-05,9.7582\mathrm{e}-05)$. Figure 1: CVA’s $T=1/4$ \- set A: Monte Carlo ($*$-) vs approximation (o-.). $\rho\backslash\gamma$ | -0.3 | -0.15 | 0 | 0.15 | 0.3
---|---|---|---|---|---
-0.8 | 1.737e-04 | 1.783e-04 | 1.811e-04 | 1.845e-04 | 1.887e-04
-0.6 | 3.288e-04 | 3.307e-04 | 3.330e-04 | 3.357e-04 | 3.391e-04
-0.4 | 4.479e-04 | 4.492e-04 | 4.509e-04 | 4.531e-04 | 4.559e-04
-0.2 | 5.332e-04 | 5.339e-04 | 5.350e-04 | 5.366e-04 | 5.388e-04
0.0 | 5.846e-04 | 5.847e-04 | 5.853e-04 | 5.864e-04 | 5.881e-04
0.2 | 6.022e-04 | 6.017e-04 | 6.018e-04 | 6.024e-04 | 6.037e-04
0.4 | 5.861e-04 | 5.851e-04 | 5.847e-04 | 5.849e-04 | 5.858e-04
0.6 | 5.365e-04 | 5.349e-04 | 5.340e-04 | 5.339e-04 | 5.347e-04
0.8 | 4.536e-04 | 4.513e-04 | 4.501e-04 | 4.500e-04 | 4.493e-04
Table 3: Approximation errors, set A for CIR parameters. Here $T=1/2$, the risk-free call price is $2.2450$. The lengths of the Monte Carlo confidence intervals for each price are in the range $(1.8252\mathrm{e}-04,3.0014\mathrm{e}-04)$. Figure 2: CVA’s $T=1/2$ \- set A: Monte Carlo ($*$-) vs approximation (o-.). $\rho\backslash\gamma$ | -0.3 | -0.15 | 0 | 0.15 | 0.3
---|---|---|---|---|---
-0.8 | 2.1015e-04 | 1.8692e-04 | 1.7121e-04 | 1.5363e-04 | 1.3338e-04
-0.6 | 5.2435e-04 | 5.3558e-04 | 5.4804e-04 | 5.6225e-04 | 5.7867e-04
-0.4 | 1.0766e-03 | 1.0844e-03 | 1.0937e-03 | 1.1047e-03 | 1.1179e-03
-0.2 | 1.4554e-03 | 1.4600e-03 | 1.4662e-03 | 1.4743e-03 | 1.4845e-03
0.0 | 1.6613e-03 | 1.6627e-03 | 1.6660e-03 | 1.6712e-03 | 1.6787e-03
0.2 | 1.6947e-03 | 1.6930e-03 | 1.6934e-03 | 1.6960e-03 | 1.7011e-03
0.4 | 1.5560e-03 | 1.5513e-03 | 1.5489e-03 | 1.5491e-03 | 1.5523e-03
0.6 | 1.2457e-03 | 1.2379e-03 | 1.2331e-03 | 1.2313e-03 | 1.2332e-03
0.8 | 7.6526e-04 | 7.5409e-04 | 7.4722e-04 | 7.4475e-04 | 7.4117e-04
Table 4: Approximation errors, set A for CIR parameters. Here $T=1$, the risk-
free call price is $3.1742$. The lengths of the Monte Carlo confidence
intervals for each price are in the range
$(4.9156\mathrm{e}-04,9.4118\mathrm{e}-04)$. Figure 3: CVA’s $T=1$ \- set A:
Monte Carlo ($*$-) vs approximation (o-.).

###  7.2 Results - Set B

$\rho\backslash\gamma$ | -0.3 | -0.15 | 0 | 0.15 | 0.3
---|---|---|---|---|---
-0.8 | 1.6115e-04 | 1.5887e-04 | 1.5728e-04 | 1.5560e-04 | 1.5378e-04
-0.6 | 6.7428e-05 | 6.6354e-05 | 6.5250e-05 | 6.4024e-05 | 6.2642e-05
-0.4 | 1.1889e-06 | 1.7841e-06 | 2.4470e-06 | 3.2382e-06 | 4.2046e-06
-0.2 | 4.5520e-05 | 4.5687e-05 | 4.5942e-05 | 4.6338e-05 | 4.6921e-05
0.0 | 6.5712e-05 | 6.5452e-05 | 6.5321e-05 | 6.5344e-05 | 6.5567e-05
0.2 | 6.1873e-05 | 6.1201e-05 | 6.0678e-05 | 6.0338e-05 | 6.0237e-05
0.4 | 3.4001e-05 | 3.2921e-05 | 3.2032e-05 | 3.1366e-05 | 3.0991e-05
0.6 | 1.7912e-05 | 1.9391e-05 | 2.0617e-05 | 2.1556e-05 | 2.2129e-05
0.8 | 9.3819e-05 | 9.5716e-05 | 9.7231e-05 | 9.8323e-05 | 9.9291e-05
Table 5: Approximation errors, set B for CIR parameters. Here $T=1/4$, the risk-free call price is $1.5877$. The lengths of the Monte Carlo confidence intervals for each price are in the range $(1.4796e-05,4.6716e-05)$. Figure 4: CVA’s $T=1/4$ \- set B: Monte Carlo ($*$-) vs approximation (o-.). $\rho\backslash\gamma$ | -0.3 | -0.15 | 0 | 0.15 | 0.3
---|---|---|---|---|---
-0.8 | 8.2741e-04 | 8.1698e-04 | 8.0882e-04 | 8.0048e-04 | 7.9149e-04
-0.6 | 3.8350e-04 | 3.7768e-04 | 3.7173e-04 | 3.6546e-04 | 3.5854e-04
-0.4 | 6.5879e-05 | 6.2609e-05 | 5.8946e-05 | 5.4764e-05 | 4.9882e-05
-0.2 | 1.2894e-04 | 1.2985e-04 | 1.3128e-04 | 1.3333e-04 | 1.3618e-04
0.0 | 2.0179e-04 | 2.0055e-04 | 1.9992e-04 | 1.9991e-04 | 2.0072e-04
0.2 | 1.5314e-04 | 1.4983e-04 | 1.4707e-04 | 1.4509e-04 | 1.4416e-04
0.4 | 1.7027e-05 | 2.2576e-05 | 2.7401e-05 | 3.1240e-05 | 3.3776e-05
0.6 | 3.0968e-04 | 3.1744e-04 | 3.2417e-04 | 3.2963e-04 | 3.3338e-04
0.8 | 7.2598e-04 | 7.3603e-04 | 7.4441e-04 | 7.5091e-04 | 7.5625e-04
Table 6: Approximation errors, set B for CIR parameters. Here $T=1/2$ and the price of the risk-free call is $2.2450$. The lengths of the Monte Carlo confidence intervals for each price are in the range $(4.0954\mathrm{e}-05,1.6693\mathrm{e}-04)$. Figure 5: CVA’s $T=1/2$ \- set B: Monte Carlo ($*$-) vs approximation (o-.). $\rho\backslash\gamma$ | -0.3 | -0.15 | 0 | 0.15 | 0.3
---|---|---|---|---|---
-0.8 | 3.0554e-03 | 3.0129e-03 | 2.9782e-03 | 2.9431e-03 | 2.9057e-03
-0.6 | 1.3068e-03 | 1.2805e-03 | 1.2548e-03 | 1.2278e-03 | 1.1989e-03
-0.4 | 1.1803e-04 | 1.0141e-04 | 8.4800e-05 | 6.7575e-05 | 4.7144e-05
-0.2 | 5.2478e-04 | 5.3076e-04 | 5.3764e-04 | 5.4572e-04 | 5.5712e-04
0.0 | 6.2528e-04 | 6.2068e-04 | 6.1724e-04 | 6.1659e-04 | 6.1960e-04
0.2 | 1.8308e-04 | 1.6810e-04 | 1.5603e-04 | 1.4691e-04 | 1.4118e-04
0.4 | 8.0508e-04 | 8.2955e-04 | 8.5101e-04 | 8.6852e-04 | 8.8145e-04
0.6 | 2.3448e-03 | 2.3792e-03 | 2.4097e-03 | 2.4354e-03 | 2.4540e-03
0.8 | 4.4474e-03 | 4.4929e-03 | 4.5321e-03 | 4.5636e-03 | 4.5870e-03
Table 7: Approximation errors, set B for CIR parameters. Here $T=1$ and the
price of the risk-free call is $2.2450$. The lengths of the Monte Carlo
confidence intervals for each price are in the range
$(1.2015e-04,6.0014e-04)$. Figure 6: CVA’s $T=1$ \- set B: Monte Carlo ($*$-)
vs approximation (o-.).

(a)

(b)

Figure 7: Approximation error as a function of the correlation $\rho$ and the
Hurst parameter $H$ for sets A and B, $T=0.25$.

## References

  * [1] E. Alòs: A decomposition formula for option prices in the Heston model and applications to option pricing approximation, Finance and Stochastics, 16, 3, 403–422, (2012).
  * [2] E. Alòs, C. Ewald: Malliavin differentiability of the Heston volatility and applications to option pricing, Advances in Applied Probability, 40, 1, 144-162, (2008).
  * [3] E. Alòs, J. Leòn: An Intuitive Introduction to Fractional and Rough Volatilities, Mathematics, 9(9), 994, (2021).
  * [4] E. Alòs, J. Leòn, J. Vives: On the short-time behavior of the implied volatility for jump-diffusion models with stochastic volatility, Finance and Stochastics, Springer, 11, 4, 571–589, (2007).
  * [5] E. Alòs, E.G. Lorite: Malliavin Calculus in Finance. Theory and Practice, CRC Press, Taylor & Francis, (2021).
  * [6] E. Alòs, K. Shiraya: Estimating the Hurst parameter from short term volatility swaps: a Malliavin calculus approach, Finance and Stochastics, 23, 423–447, (2019).
  * [7] E. Alòs, F. Antonelli, A. Ramponi, S. Scarlatti: CVA and vulnerable options in stochastic volatility models, International Journal of Theoretical and Applied Finance, 24, 2, 1–34, (2021).
  * [8] F. Antonelli, A. Ramponi, S. Scarlatti: CVA and vulnerable options pricing by correlation expansions, Annals of Operation Research, 299, 401–427, (2021).
  * [9] C. Bayer, P. Friz, J. Gatheral: Pricing under rough volatility, Quantitative Finance, 16, 6, 887–904, (2016).
  * [10] T. R. Bielecki, D. Brigo, S. Crepey, Counterparty Risk and Funding: A Tale of Two Puzzles. Chapman and Hall/CRC, (2014).
  * [11] T. R. Bielecki, M. Rutkowski: Credit Risk: Modeling, Valuation and Hedging, Springer Finance Series, (2002).
  * [12] T. R. Bielecki, M. Jeanlanc, M. Rutkowski: Valuation and Hedging of Credit Derivatives, Lecture notes CIMPA- UNESCO Morocco School, (2009).
  * [13] D. Brigo, F. Vrins:  Disentangling wrong-way risk: pricing credit valuation adjustment via change of measures, In Innovations in Insurance, Risk- and Asset Management, WSPC Proceedings ,(2018).
  * [14] D. Brigo, T. Hvolby, F. Vrins: Wrong-Way Risk adjusted exposure: Analytical Approximations for Options in Default Intensity Models, European Journal of Operational Research, 269, 1154–1164, (2018).
  * [15] F. Comte, E. Renault: Long memory in continuous-time stochastic volatility models, Mathematical Finance, 8, 291-323, (1998).
  * [16] D. Duffie, K. Singleton: Modeling Term Structures of Defaultable Bond, Review of Financial Studies, 12, 3, (1999).
  * [17] F. A. Fard: Analytical pricing of vulnerable options under a generalized jump-diffusion model, Insurance Mathematics and Economics, 60, 19–28, (2015).
  * [18] P. V. Gapeev: Some extensions of Norros’ lemma in models with several defaults. Inspired by Finance, The Musiela Festschrift. Kabanov Yu. M., Rutkowski M., Zariphopoulou Th. eds. Springer, 273–281, (2014).
  * [19] P. V. Gapeev, M. Jeanblanc, L. Li, M.Rutkowski: Constructing random measures with given survival processes and applications to valuation of credit derivatives, Contemporary Quantitative Finance, Essays in Honour of Eckhard Platen. Chiarella, C., Novikov, A. eds. Springer, 255–280, (2010).
  * [20] J. Gatheral, T. Jaisson, M. Rosenbaum: Volatility is rough, Quantitative Finance, 18, 6, 933–949, (2018).
  * [21] A. Gnoatto, C. Reisinger, A. Picarelli, Athena: Deep xVA Solver ? A Neural Network Based Counterparty Credit Risk Management Framework http://dx.doi.org/10.2139/ssrn.3594076, (2021).
  * [22] P. S: Hagan, D. Kumar, A.S. Lesniewski, D. E. Woodward: Managing Smile Risk, Wilmott Magazine, 84–108, (2002).
  * [23] S. L. Heston, A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options: The Review of Financial Studies, 6, 2, 327–343, (1993).
  * [24] P. Klein: Pricing Black-Scholes options with correlated credit risk, Journal of Banking & Finance, 20, 1211–1229, (1996)
  * [25] D. Lamberton, B. Lapeyre: Introduction to Stochastic Calculus applied to Finance, Chapman and Hall/CRC Financial Mathematics Series, 2nd ed., (2007).
  * [26] D. Lando: On Cox Processes and Credit Risky Securities, Review of Derivatives Research, 2, 99–120, (1998).
  * [27] A. Nikeghbali: An essay on the general theory of stochastic processes, Probability Surveys, 3, 345-412, (2006).
  * [28] D. Nualart: The Malliavin Calculus and Related Topics, Probability and its Applications series, 2nd ed. Springer, (2006).
  * [29] S. Zhu, M. Pykhtin: A Guide to Modeling Counterparty Credit Risk, GARP Risk Review, July/August (2007).

##  Appendix A Malliavin Calculus

Here we introduce the main concepts of Malliavin Calculus used in this work.
We refer to [28] for a deeper approach to this field.

Let ${\mathcal{S}}$ be the set of random variables of the form

$F=f(W(h_{1}),\ldots,W(h_{n})),$ (48)

with $h_{1},\ldots,h_{n}\in L^{2}([0,T])$, $W(h_{i})$ denotes the Wiener
integral of the function $h_{i}$, for $i=1,..,n$, and $f\in
C_{b}^{\infty}(\mathbb{R}^{n})$ (i.e., $f$ and all its partial derivatives are
bounded). The Malliavin derivative of $F$, $DF$, is defined as the stochastic
process in $L^{2}(\Omega\times[0,T])$ given by

$D_{s}F=\sum_{j=1}^{n}{\frac{\partial f}{\partial
x_{j}}}(W(h_{1}),\ldots,W(h_{n}))(s)h_{j}(s).$

We can also define the iterated Malliavin derivative operator as

$D_{s_{1},\ldots,s_{m}}^{m}F=D_{s_{1}}\ldots D_{s_{m}}F,\qquad
s_{1},\ldots,s_{m}\in[0,T].$

These operators are closable in $L^{p}(\Omega)$ for any $p\geq 1$ and we
denote by ${\mathbb{D}}^{n,p}$ the closure of ${\mathcal{S}}$ with respect to
the norm

$||F||_{n,p}=\left(E\left|F\right|^{p}+\sum_{i=1}^{n}E||D^{i}F||_{L^{2}([0,T]^{i})}^{p}\right)^{\frac{1}{p}}.$

We will also make use of the notation
$\mathbb{L}^{n,p}=\mathbb{D}_{W}^{n,p}(L^{2}([0,T])).$

The adjoint of the derivative operator $D$ is the divergence operator
$\delta$, that coincides with the Skorohod integral. That is, the domain of
$\delta$, denoted by Dom $\delta$, is the set of processes $u\in
L^{2}(\Omega\times[0,T])$ such that there exists $\delta(u)\in L^{2}(\Omega)$
such that

$E(\delta(u)F)=E\left(\int_{0}^{T}(D_{s}F)u_{s}ds\right),\qquad\hbox{\rm for
every}\;\;F\in{\mathcal{S}}.$ (49)

We use the notation $\delta(u)=\int_{0}^{T}u_{s}dW_{s}$. It is known that
$\delta$ is an extension of the classical Itô integral in the sense that
$\delta$, applied to adapted and square integrable processes, coincides with
the classical Itô integral. We also know that the space ${\mathbb{L}}^{1,2}$
is included in the domain of $\delta$.

The following result gives us an explicit martingale representation for random
variables in $\mathbb{D}_{W}^{1,2}$.

######  Theorem A.1 (The Clark-Ocone formula)

Consider a random variable $F\in\mathbb{D}_{W}^{1,2}$.Then

$F=E(F)+\int_{0}^{T}E_{r}(D_{r}F)dW_{r}.$

Finally, we recall the anticipating Itô’s formula

###### Proposition A.1

Assume the model (10) where $\sigma\in\mathbb{L}^{1,2}$ and consider the
processes $Y,Z$ defined $Y_{t}=\int_{t}^{T}\sigma^{2}_{s}ds$. Let
$F:[0,T]\times\mathbb{R}^{3}\rightarrow\mathbb{R}$ be a function in
$C^{1,2}([0,T]\times\mathbb{R}^{2})$ such that there exists a positive
constant $C$ such that, for all $t\in\left[0,T\right],$ $F$ and its partial
derivatives evaluated in $\left(t,X_{t},Y_{t},\right)$ are bounded by $C.$
Then it follows that

$\displaystyle F(t,X_{t},Y_{t})$ $\displaystyle=$ $\displaystyle
F(0,X_{0},Y_{0})+\int_{0}^{t}{\partial_{s}F}(s,X_{s},Y_{s})ds$ (50)
$\displaystyle-\frac{1}{2}\int_{0}^{t}{\partial_{x}F}(s,X_{s},Y_{s})\sigma^{2}_{s}ds$
$\displaystyle+\int_{0}^{t}{\partial_{x}F}(s,X_{s},Y_{s})\sigma_{s}dW_{s}$
$\displaystyle-\int_{0}^{t}{\partial_{y}F}(s,X_{s},Y_{s})\sigma^{2}_{s}+\int_{0}^{t}{\partial_{xy}^{2}F}(s,X_{s},Y_{s})\Theta_{s}ds$
$\displaystyle-\int_{0}^{t}{\partial_{z}F}(s,X_{s},Y_{s})a^{2}_{s}$
$\displaystyle+\frac{1}{2}\int_{0}^{t}{\partial_{xx}^{2}F}(s,X_{s},Y_{s},Z_{s})\sigma^{2}_{s}ds$

where $\Theta_{s}:=(\int_{s}^{T}D_{s}\sigma_{r}^{2}dr)\sigma_{s}$.