A Gamma Ornstein–Uhlenbeck model driven by a Hawkes process

We propose an extension of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}Γ-OU Barndorff-Nielsen and Shephard model taking into account jump clustering phenomena. We assume that the intensity process of the Hawkes driver coincides, up to a constant, with the variance process. By applying the theory of continuous-state branching processes with immigration, we prove existence and uniqueness of strong solutions of the SDE governing the asset price dynamics. We propose a measure change of self-exciting Esscher type in order to describe the relation between the risk-neutral and the historical dynamics, showing that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}Γ-OU Hawkes framework is stable under this probability change. By exploiting the affine features of the model we provide an explicit form for the Laplace transform of the asset log-return, for its quadratic variation and for the ergodic distribution of the variance process. We show that the proposed model exhibits a larger flexibility in comparison with the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}Γ-OU model, in spite of the same number of parameters required. We calibrate the model on market vanilla option prices via characteristic function inversion techniques, we study the price sensitivities and propose an exact simulation scheme. The main financial achievement is that implied volatility of options written on VIX is upward shaped due to the self-exciting property of Hawkes processes, in contrast with the usual downward slope exhibited by the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}Γ-OU Barndorff-Nielsen and Shephard model.


Introduction
In recent years, the implied volatility indices, such as VIX in the US and V2X in Europe, proved themselves to be key financial instruments for investment strategies, hedging (see [9,55]) and also indicators of the "stress" on the market [35]. This growing importance has shed light on the peculiar form of the implied volatility dynamics: the occurrence of large variations on very short periods, with a tendency to form clusters of spikes. Moreover, new derivatives products appear, due to market regulation and standardization, like volatility options, with a significant and increasing demand. Implied volatility of these product exhibits an upward slope, this is a clear evidence of market volatility risk aversion pushing people to buy options to cover this risk [53]. The existence of such stylized facts accounts for the emergence of a wide range of financial models taking into account these features. Heston [34] proposes a model where the variance is stochastic and follows a CIR [17] diffusion. By using Poisson processes, Bates [7] adds jumps in asset dynamics, while Sepp [58] includes jumps in both the asset returns and the variance. Duffie et al. [22] generalize the jump-diffusion framework by including a stochastic intensity for the jump processes both in the asset returns and the volatility dynamics, however no self-exciting example is proposed in their paper. Kallsen et al. [43] consider the case where jumps are added in the stock via a time-changed Lévy process. More recently, a large literature develops rough volatility models ( [8,23] and [24]), while Abi Jaber et al. [2] propose a model based on the generalization of affine Volterra processes.
The aim of this paper is to propose a new model for asset pricing able to capture the main stylized features but preserving both mathematical and numerical tractability. Our idea is to build up our model as an extension of the Barndorff-Nielsen and Shephard (BNS) model in order to include jumps clusters in both the volatility and the stock return dynamics. In their celebrated papers Barndorff-Nielsen and Shephard [5,6] propose a stochastic volatility model of Ornstein-Uhlenbeck type driven by a subordinator; by considering as a concrete specification for the subordinator a compound Poisson process with exponentially distributed jumps size, they obtain a model where both the variance and the log-returns are driven by the same jump process. This type of stochastic volatility model, named -OU model and originally introduced in order to fit empirical data, was later shown to be suitable for option pricing by Nicolato and Venardos [54]. As pointed out by a large amount of literature (e.g. [59]), single factor models, and, more specifically, the BNS model, are not flexible enough and they are outperformed by multi-factor stochastic volatility models in practice (a multifactor extension of the BNS model is proposed in [52]). On the other hand, multi-factor models have a main drawback in the huge number of parameters required, giving birth to fitted parameters instability (e.g. [7]).
We propose to extend the -OU model in order to include jump clustering features. The clustering effect of jumps of the Hawkes [33] process is well suited to take into account the periods of turmoil in the implied volatility, typically observed in implied volatility indices. These clustering features have been studied throughout various financial asset models: see Fiura [28] for FX rates, Hainaut [32] and Jiao et al. [39] for interest rates, Abergel and Jedidi [1], Horst and Xu [36], Zheng et al. [61] for microstructure and limit order books, Errais et al. [25] for credit risk, Jiao et al. [40] for energy prices and Granelli and Veraart [30] for risk premium and contagion. In this paper, we shall insist on the stylized facts related to the clustering effects of the implied volatility indices, such as VIX index, and the related options.
In the model we are going to propose, the jump component of the dynamics is taken into account through a marked Hawkes process, with exponentially distributed jump sizes. Thus, the model has three state variables, the stock, the variance process and the intensity of the jumps. In order to prevent parameters instability, we assume that the intensity process coincides with the variance process up to an additive constant, representing the minimal intensity of jumps arrival. Our model, that will be referred to as -OU Hawkes volatility model, shares the same number of parameters of -OU since the constant intensity is replaced by the shifting constant linking the intensity and the variance processes. In spite of this, the -OU Hawkes model is intrinsically multi-factor, since both variance and intensity are stochastic processes, so it exhibits more flexibility than the -OU model. Parsimony is preserved by the fact that the two processes coincide up to an additive constant. By exploiting a measure change of Esscher type in a self-exciting setting, as introduced in Jiao et al. [39,40], we describe the relation between the risk-neutral and the historical dynamics, showing that the -OU Hawkes framework is stable under this probability change unlike the α-stable case in Jiao et al. [39,41].
We show that it is possible to fit the -OU Hawkes model to obtain the skew, that is the slope and the left wing of the implied volatility. We then perform a calibration of our model using plain vanilla market data and provide a sensitivity analysis in order to illustrate the flexibility of the model. This sensitivity analysis shows, in particular, that two parameters, namely the speed of mean reversion and the minimal jump intensity, have a small impact on the implied volatility. By using extreme values, we observe that the main impact of both the speed of mean reversion and the minimal jump intensity is detected around the money and consists, roughly speaking, in increasing the smile and preserving the slope for the left side when the minimal jump intensity decreases or the speed of mean reversion increases. This result could appear counterintuitive at first sight since, without jumps, the -OU Hawkes model, as the BNS, reduces to a Black Scholes diffusion and then the smile is expected to disappear. The explanation is contained in the self-exciting structure of the -OU Hawkes model since, when the minimal jump intensity decreases, the endogeneity of the jumps increases and that explains why extreme events occur more often.
We also point out that moment explosion does not take place in -OU Hawkes model. This interesting feature is crucial for several reasons: first, when Monte Carlo simulations are needed in order to perform derivatives evaluation with no close-form solutions available, an infinite variance can constitute a serious drawback. A similar problem arises when considering derivatives with a nonlinear payoff. Third, in dynamic portfolio optimization moment explosion problems can give rise to an infinite value function in many usual frameworks, which represents a main drawback in financial applications.
By focusing on the options written on VIX, the same self-exciting effect gives birth to an upward sloping implied volatility. This upward slope behavior is coherent with market data (e.g. [53]), but is exhibited by very few models in the literature. Nicolato et al. [53] study the implied volatility of variance options for different stochastic volatility models with jumps and show that only inverse-Gamma Ornstein-Uhlenbeck is able to reproduce an upward slope. Indeed, in the case of exponential law for jumps size, the implied volatility for options written on VIX is down-sloping. We then conclude that, in the -OU Hawkes model, the upward slope is a consequence of self-exciting structure of jumps. We also provide a sensitivity analysis for VIX implied volatility.
We conclude the analysis of options by performing a calibration of the model. The calibration is performed first on Eurostoxx 50 options for two maturities before the COVID crisis. This first set of parameters is exploited for sensitivities analysis. Then, we perform a second round of calibrations on both options on S&P500 and VIX index for a very short maturity (1 month). Calibrations are performed separating S&P500 and VIX options giving birth to two new sets of calibrated parameters. These parameters are similar excepting by the leverage parameter ρ that is close to 0 for VIX option. We link this effect to the definition of VIX index given in the Chicago Board Options Exchange [16] white paper which exploits a static replication valid for continuous paths underlying [21]. Surprisingly the -OU Hawkes volatility model calibrated on VIX options tries to reproduce the continuous paths forcing a negligible leverage.
Finally, we address the problem of the simulation of the proposed model, which is a very important aspect for its practical applicability. Exact simulation schemes for stochastic volatility models have been first introduced in literature by Broadie and Kaya [12] which show how to simulate exactly the transition density of the Heston (and its extensions with jumps) model. Cai et al. [13] and Li and Wu [49] develop a similar approach for, respectively, SABR [31] and Ornstein-Uhlenbeck stochastic volatility [57] models. The main drawback of those methods is that they rely on the numerical inversion of conditional Laplace transforms which make the simulation schemes slow and hard to implement (see, e.g., [46]). An important advantage of the proposed -OU Hawkes volatility model is that it can be simulated exactly without resorting to any numerical method. We propose an exact simulation scheme for the model and evaluate the performances by a comparison with the classic Euler scheme.
The rest of the paper is organized as follows. Section 2 details a statistical analysis of VIX values pointing out the existence of jump clusters and giving a first fit to data. In Sect. 3 we introduce the -OU Hawkes volatility model by providing theoretical results about existence and characterization in a general framework. After presenting the affine properties of the -OU Hawkes model, we discuss the moment explosion issue. Moreover, we compute the Laplace transform of the quadratic variation of the process involved, and provide a closed form expression for the variance swaps rates. After introducing an Esscher type transform for the present model, we provide a class of equivalent martingale measures. Section 4 is devoted to numerical applications, i.e. European pricing via characteristic function inversion techniques, calibration, perfect simulation and VIX options analysis. Finally, Sect. 5 summarizes the main properties of -OU Hawkes volatility model, and presents a systematic comparison with both the BNS and the Heston models.

Clusters in VIX: stylized facts
In this section, we illustrate and discuss the clustering effects of the VIX. A general analysis of this index is performed by Avellaneda and Papanicolaou [4]. This point provides us with an empirical justification of our approach involving clustering effects both on the jumps of the asset price and on the volatility. Moreover, we also see the particular importance of downward jumps. The VIX index represents the square-root of the implied variance extracted from short dated options on the S&P 500 index. Hence, we consider the square of VIX index since it coincides with the variance swap rate, that is a linear functional with respect to time, converging to the quadratic variation of the logarithm of the equity process (see, e.g. [43] for more details).
Our  Fig. 1. Since the sample starts in a period of mild volatility we can reasonably assume that the Hawkes intensity, at this time, is close to its minimal level, which we denote by λ. The jumps are detected as the larger positive fluctuations using the algorithm detailed in Callegaro et al. [14].
The first analysis concerns the distribution of gaps between two jumps in order to confirm that the jumps are clustered. Given the total number of jumps in a given period and assuming a constant intensity, the conditional distribution of gap lengths is uniform. It is now easy to test this hypothesis using a Kolmogorov-Smirnov (K-S) test. The renormalised K-S statistics gives a value of 4.45 that is really large compared to all critical values usually considered: for instance, for a rejection of the null hypothesis at level of 0.01, the critical value is 1.62, for 0.001 it is 1.95. We then reject the Poisson arrival rate of jumps. In contrast, we can test the hypothesis of an intensity proportional to the index VIX 2 itself and the related renormalised K-S statistics gives a value of 1.04 that is small compared to the critical values usually considered, for instance 1.22 for a significance threshold 0.1. As a consequence, we cannot reject the pure self-excited framework. The top right subplot of Fig. 1 illustrates the result of the goodness of fit. It can be remarked that the values obtained assuming a constant intensity are really far from the ideal diagonal whereas the ones assuming an intensity proportional to the VIX 2 itself are much closer.
We now turn on the joint analysis of S&P and VIX 2 . We identify, in an independent way, the jump times and the size of the relative increments of S&P and the absolute increments of VIX 2 . Our study shows that almost all large negative jumps of S&P coincide with the positive jumps in VIX 2 , whereas only one half of positive jumps of S&P coincides with the negative ones in VIX 2 . By studying the arrival times of the negative jumps of VIX 2 , we see that they coincide with the very large values of the VIX 2 itself and follow large positive jumps. We can reproduce this effect with only positive jumps in VIX and an exponential mean-reversion speed. Moreover, the method adopted to identify jumps is unable to detect relatively small jumps (see [14]). In particular, jumps smaller than three standard deviations of the other increments are classified as usual Brownian noise. The fact that we identify almost 2.5 more positive jumps in VIX 2 than negative jumps in S&P could be explained by the presence of a larger Brownian contribution, covering a part of relatively small jumps. Indeed, we stop our iteration at the fifth loop since the ratio of the variance at the fourth and fifth loop in S&P is 0.975, i.e. really near to 1, showing that the split between the jumps and the Brownian component is done. However, the same ratio for VIX 2 increments is only 0.83 showing that the Brownian part is certainly overestimated. In order to reach a variance ratio for VIX 2 similar to the one of S&P more than eleven loops are required. The total number of large fluctuations reaches a frequency of one every ten days that pushes us to cut the Brownian component in the variance process in agreement with the BNS model.
We now focus on the law of jumps sizes in VIX 2 and S&P. According to the Kou [47,48] setup we obtain a quite acceptable fit for an exponential law truncated at the threshold used to split increments between jumps and Brownian oscillations. The Kolmogorov-Smirnov test gives values of 1.06 and 1.17 respectively for the positive jumps of VIX 2 and negative ones of S&P, that could not authorize to reject the exponential law hypothesis. Results on the joint S&P and VIX 2 analysis are reported in the bottom subplots of Fig. 1.

The 0-OU Hawkes volatility model
As the BNS model, the -OU Hawkes volatility model is a stochastic volatility model with the same jumps driving both the volatility and the log-return processes, including a leverage effect. We construct the variance process by subtracting a strictly positive constant λ to the Hawkes intensity process (which is bounded from below by λ itself), in this way the variance process exhibits the suited mean-reverting form and can visit all the positive values. Since the model is constructed in order to reproduce the VIX 2 behavior, which is naturally defined with respect to the risk neutral probability, we start by considering the risk neutral framework, the historical probability will be then introduced in Sect. 3.5.

Definitions
Let ( , F , Q) be a probability space supporting a standard Brownian motion W := {W t } t≥0 and a marked Hawkes process independent on W . The Hawkes process can be characterized by its counting measure μ(dt, dz) defined on (R + ) 2 and by the sequence of jump times {τ i } i∈N and marks {Z i } i∈N , representing the jump times and jump sizes respectively. The compensator of μ can be written as λ t θ(dz)dt, where θ(dz) is a distribution on the measurable space R + , and the compensated measure is denoted byμ(dz, dt) = μ(dz, dt) − λ t θ(dz)dt.
In this paper, we assume that θ is of exponential type θ(dz) := ναe −αz dz. The intensity process reads: where λ < λ 0 and β are positive constants. From now on, we require that the following integrability assumption holds: The following condition holds: Under our assumption on the jump size distribution, (2) can be written asβ := β− ν α > 0. This condition could be considered as the usual non-explosion condition for Hawkes processes in a marked Hawkes framework (see [10]).
The stock price process {S t } t∈R + is defined via its log-return process: {X t } t∈R + : X t := ln(S t /S 0 ). In order to describe the leverage effect mentioned before, the jump process driving the log-return dynamics is multiplied by −ρ, with ρ > 0. Without loss of generality, we will assume that the interest and dividend rates are null. We assume that X satisfies the following SDE: which under our assumption on the jump size distribution becomes: where γ = νρ ρ+α and the variance process is given by Note that, thanks to the negative sign in front of the jump term and the domain of Z i , jumps in the log-return process are negative and the compensator is positive. The natural filtration of the point process N will be denoted by F N := F N t t≥0 , while the natural filtration of (λ, S) will be denoted by F := {F t } t≥0 .

Existence and basic properties
Next, we shall prove existence and regularity of the solution of the SDE, by using the tools of continuous state branching processes with immigration (CBI). In particular, we will show the existence of a strong solution under mild hypotheses including the -OU Hawkes volatility model. We extend in particular the evolution by assuming a possible infinite activity of the jump component. The -OU Hawkes volatility model will be obtained under the hypothesis of finite activity and an exponential law for the jump size. The advantage of this generalisation is to provide a unified framework as in Barndorff-Nielsen and Shephard [5,6], by including both finite and infinite activity. Moreover, the representation exploiting Branching processes features provides valuable insights and ideas for interpretation of the model behaviour. We first write the evolution of the compensated version of the SDE satisfied by the triplet (λ, σ, X ). We then characterize the evolution as a multi-dimensional exponential affine process where λ is an autonomous CBI. A byproduct of this characterization is that our model could be rewritten in the so-called Dawson and Li representation of CBI (see [19,20], [51]). This last characterization provides access to the ergodic distribution of the process λ and σ . Under hypothesis 1, the process λ is exponential ergodic. The law of the invariant distribution and its parameters will be also derived.

Lemma 1 (Compensated representation)
The triplet (λ, σ, X ) satisfies the following SDE where the jump component is compensated.
The proof is a straightforward computation. The crucial point is that the compensator is proportional to λ. Hereafter, in this section, we assume that the compensated measure satisfies the integrability condition This SDE admits a unique solution, which coincides almost surely with the solution of (5).
The main advantage of this representation is to highlight that the speed of mean reversion and level of the intensity between two jumps change due to the self-exciting property. This representation is needed in order to obtain the change of probability result in Sect. 3.5. Moreover, one of the consequences of the previous result is the affine structure of the couple (λ, X ). The next result characterizes its Fourier-Laplace transform.

Proposition 2 (Fourier-Laplace transform)
The couple (λ, X ) is an exponential affine process. That is the Laplace transform of (λ, X ) satisfies: denotes the real part of a complex number. The functions ψ and φ satisfy: with starting conditions ψ u,w (0) = w and φ u,w (0) = 0. In particular, considering the Fourier transform, the pseudo Riccati operator reads R(iu, iw) = − 1 Proof A necessary condition for the Laplace transform to exist is the integrability of the jumps, that is In Proposition 3, we shall show that this is actually a sufficient condition. In the following, we consider the operator R and F only on this domain. The main statement is a direct application of Duffie et al. [22,Proposition 1]. By adopting the notation of Kallsen et al. [43], we can write the differential characteristics (b, c, K ) of the process (λ, X ). Since we consider the compensated version, we do not need to specify any particular truncation function. (5) and recall that the relations between the differential characteristics (b, c, K ) and the affine characteristics (F, R) are the following: where (u, w) T denotes the transposed of the column vector (u, w). By remembering that the moment generating function evaluated at y of an exponential random variable with parameter α is y/(y − α), we obtain the Riccati system of ODE (6) for ψ and φ.
Remark 1 Proposition 2 can be further extended to obtain the Laplace transform of t 0 X s ds, t 0 λ s ds by using the results in Brignone and Sgarra [11].
For sake of readability, we will write again the generalized Riccati operator in the following form on D where p(u) and q(u) are the two following quadratic polynomials: The next result is a direct consequence of the affine structure of the model.

Corollary 2 (Ergodic distribution) Under Assumption 1, the intensity process λ is exponential ergodic and the moment generating function of the invariant distribution is given by
Simplifying, we get: from which the thesis follows.

Explicit Laplace transform and moments explosion
Next, we obtain the explicit form of the Laplace transform of the -OU Hawkes model, but, in order to evaluate its domain, we first deal with the explosion of moments. The next proposition shows that the -OU Hawkes model has moments of every positive order and the negative moments of orders larger than u > −α/ρ.
we then deduce that proving the result.
We now focus on the explicit solution of the Laplace transform of the couple (X , λ). Looking at the ODE satisfied by ψ, we remark that it is non-linear but of first order and separable. Then, we can formally solve it as indicated in the following corollary.

Corollary 3 (Explicit form for the inverse function of the Laplace transform)
Proof The first equality of (8) is a direct consequence of the ODE satisfied by ψ u,w (t), see Proposition 2. We have then to compute explicitly dy .
The first integrand has the form f (y)/ f (y) then we have That is the term L(ψ u,w (t), u, w) in (8). For the second term, we need to distinguish three cases, that is if the polynomial β y 2 + p(u)y + uq(u) has two real roots, i.e. (u) > 0, only one root, i.e. (u) = 0 or complex roots. We first consider the case (u) > 0, i.e. u ∈ I , then β dy .
By splitting and integrating the integral and recalling that w Suppose now that the two solutions coincide that is (u) = 0 and w − (u) = w + (u) = − p(u)/(2β). Then, we have 1 2β Finally, when (u) < 0, we have 1 2β From which the thesis follows.

Remark 2
The function H coincides, up to some constants, to the corresponding function of the Heston model. The function L, that does not appear in Heston case, is the logarithm of a quadratic function of ψ.

Variance swap and VIX
Next, we aim to provide closed form expressions for the variance swap rates. The VIX index could then be obtained easily. VIX is the forecasted average volatility of S&P during the next month, the equivalent index based on Eurostoxx 50 is called V2X. In the rest of the paper, the term "VIX" is used broadly to refer to VIX or V2X or all other volatility indices based on a particular underlying, where no risk of confusion exists. The next proposition, gives the explicit form of the variance swap rates.

Proposition 4 (Variance swap rates) Under -OU Hawkes volatility model, the variance swap rate at time t with a time to maturity t + T reads
where Proof Recalling that the -OU Hawkes volatility model belongs to exponential affine class, we can obtain the Laplace transform of the quadratic variation thanks to Kallsen et al. [43,Lemma 4.2]. It yields that, for (u) ∈ R − , where U denotes the confluent hypergeometric function of the second kind. We can obtain the expected value of the quadratic variation by differentiation of the Laplace transform and setting u = 0. As in Kallsen et al. [43,Lemma 4.2], we obtain, under Assumption 1, 0). In order to obtain the explicit form of 1 (T ), we differentiate (10) and, taking u = 0, we obtain the differential equation satisfied by 1 . It reads By solving the previous linear equation and recalling the initial condition 1 (0) = 0, we obtain the explicit form for 1 (T ) given in the statement. Similarly, differentiating the relation satisfied by 0 and taking u = 0, we obtain 0 (T ) = T 0 βλ 1 (s)ds + λ T and a direct integration gives the explicit form for 1 (T ) in the statement.
We can obtain the following corollary since the VIX index could be expressed as the square root of the variance swap rate with maturity one month (see e.g. [4,41,53,55] and references therein).

Corollary 4
Let T equal to one month, i.e. T := 1/12, then VIX index reads V I X t = √ 1 (T )λ t + 0 (T ). Moreover, the forward rate of VIX index reads The proof of the last statement is based on the Laplace transform of square root function using the Laplace transform given in Proposition 2. Equation (11) is almost closed given the result of Corollary 3.

Change of probability
In this section, we investigate the change of probability in the framework defined above. One of the main advantages of the integral representation detailed in Proposition 1 is that it gives rise to a natural extension of the Esscher transform in the jump clustering framework. The main result is detailed in the following proposition showing that the class of -OU Hawkes models is stable under a self-exciting Esscher type change of probability (see [39,40] for similar results in the α-stable case, which is not stable under the same kind of probability change). We point out that the Esscher transform that we are going to introduce is not the Esscher transform for the semimartingale S t , as defined by Kallsen and Shiryaev [44], but it is the Esscher transform of the driving Hawkes process, strictly analogous to the transform defined by Nicolato and Venardos [54, eq. 3.14]. For a critical presentation of the Esscher transform for BNS models we refer to Hubalek and Sgarra [37]. In this section, in order to avoid ambiguity, we add a superscript P or Q on the various quantities of (5) depending on the reference probability.
Proposition 5 (Self-exciting Esscher transform) Let (λ, σ, X ) be as in Proposition 1 under the risk neutral probability Q. Fix (η, ξ ) ∈ R × (−α Q , ∞) and define Then, the Doléans-Dade exponential E(U ) is a martingale and the probability measure P defined by dP/dQ| F t := E(U ) t is equivalent to Q. Moreover, under P, the couple (λ, X ) satisfies the evolution of exponential affine class (see Eq. 5) with parameters and the dynamics with respect to Q takes the following form: Proof First, we remark that R + e −ξ z θ P (dz) < ∞, since (ξ ) ∈ (−α Q , ∞). It is easy to show that the triplet (λ, X , Y ), where Y := E(U ), is Markovian and exponential affine by applying the same argument of Proposition 1. Thanks to the integrability property we can then apply Kallsen and Muhle-Karbe [42,Corollary 3.9] and this will imply that Y is a martingale, that P exists and it is equivalent to Q. We have easily that dY t = Y t− dU t . For the second statement, let f ∈ C 2 b (R + × R), we apply Itō formula to H t := f (λ t , X t ) Y t . Let's denote by f λ (resp. f X ) the first derivative of f with respect to λ (resp. X ). A standard but tedious computation gives By identifying the terms, we obtain the evolution of (λ, X ) under P.

Numerical applications
This section deals with numerical studies of the -OU Hawkes volatility model. In the first part, we show how to price European call options using characteristic function inversion techniques. A calibration on a recent set of market data is proposed. Moreover, we estimate the behaviour of the slope of the implied volatility as time to maturity goes to zero and show that it exhibits a power decay with parameter 0.2. The second part deals with options written on VIX index. The main result is that VIX implied volatility increases with strike.Nicolato et al. [53] have shown that -OU model exhibits a VIX implied volatility decreasing with the strike. We deduce that the increasing VIX implied volatility is a consequence of the selfexciting structure. Then, we present some results on the calibration on S&P 500 and VIX options on the same date and maturity and compare the resulting parameters. In the third part we show how to simulate efficiently the -OU Hawkes model. In particular, we will show that the model can be simulated exactly, allowing for unbiased estimation of derivatives prices. Finally, we study numerical efficiency of the proposed method through a comparison with the Euler simulation scheme.
Computations are done using Matlab ® (Version R2019b) in Microsoft Windows 10 ® running on a machine equipped with Intel(R) Core(TM) i7-9750HQ CPU @2.60GHz and 16 GB of RAM.

European option pricing and calibration
In Sect. 3.3 we obtained the Laplace transform of (X , λ) in terms of the solution of an ODEs system. The characteristic function of log-returns under the proposed model is given by ϒ(u) := E e iu X T and can be obtained by Proposition 2. This result opens the doors to option pricing via standard characteristic function inversion algorithms such as, for example, FFT [15] and COS [26] methods. The latter has the advantage to present exponential convergence to the true solution, while preserving linear computational complexity. The main consequence is that option prices can be estimated through a smaller number of evaluations of the characteristic function, which is particularly important when it has to be computed through time consuming numerical techniques such as solution of ODEs systems (as in the present case). See Brignone and Sgarra [11] for more details. Hence, given an initial price S 0 , a strike K and maturity T , European option prices are computed as where Following Fang and Oosterlee [26] we set where L can be chosen arbitrary large and c i denotes the i-th cumulant of ln S T K . Since the model is affine, cumulants can be computed analytically, for this purpose, we adapt  Table 1 (Set A) the procedure outlined in Feunou and Okou [27]. According to Fang and Oosterlee [26], if L = 10, (14) gives a truncation error around 10 −12 , which is negligible for our purposes. Given European options prices, the corresponding implied volatility can be computed by inverting the Black-Scholes formula. Equipped with an efficient procedure for computing the implied volatility surface, we calibrate our model on the Eurostoxx 50 market data of the 19 November 2019. We minimize the differences (in absolute value) between market and model implied volatilities. The resulting parameters are displayed in Table 1 (Set A). Calibrated risk neutral density and implied volatility are plotted for different maturities in Fig. 2. We remark that the density functions have a left tail heavier than the right one especially for the shorter maturity.
We now focus on the slope of ATM implied volatility, i.e. the derivative of implied volatility with respect to the at the money strike, and its dependency on the maturity. Figure 3 confirms that the implied volatility exhibits a steeply negative ATM slope. This phenomenon is magnified for short maturities over an extremely large range, down to 5 minutes that could be considered as a limit of microstructure. According to empirical studies, we test if the slope of implied volatility skew exhibits a power decay with time to maturity. The ATM implied volatility slope in the -OU Hawkes volatility model has a power decay behaviour with parameter −0.20 and an R 2 coefficient of determination of 0.81. We stress that in the Heston model it presents a different behavior: it converges to a finite value as time to maturity goes to zero. Then, the behaviour of the -OU Hawkes volatility model is extremely different from that of the Heston model (for short maturities) and replicates better the empirical facts.
Next, we focus on the sensitivity of the implied volatility with respect to model parameters. In Fig. 4 we show the model implied volatility for different sets of parameters. We decide to consider σ 0 as an initial condition for the process. Then sensitivities with respect to λ are obtained shifting both λ and λ 0 by the same value. We will focus on the impact of a change of a parameters on the level and the slope of implied volatility for negative moneyness, the position of the minimum of the implied volatility and the smile (i.e. convexity around this  Table 1 (Set A) minimum). We first remark that an important change of λ or β is needed to distinguish their sensitivities. We will see that these parameters impact more the implied volatility for options written on VIX. For extreme changes on λ, we observe that the main impact of an increase in λ is that the position of the minimum of implied volatility moves on right and the implied volatility goes up, the slope for negative moneyness is unchanged. For extreme changes on β, we observe that the main impact of an increase in β is to accentuate the smile for positive moneyness. There is also a small negative effect on the level of implied volatility. The position of the minimum and the slope are roughly unchanged. A small change on α and/or ρ has an important effect on the slope with opposite directions, i.e. a positive change in α (resp. ρ) decreases (resp. increases) the steepness of the slope. Convexity is roughly unchanged in both the cases. ρ has a positive impact on the position of the minimum whereas α has no effect.

Options on VIX
We focus now on options written on VIX index that coincides with the square root of the variance swap rate. Of course, this analysis also applies to the equivalent index based on the Eurostoxx 50, the V2X index. According to Proposition 4, the variance swap rate is affine in λ and given by (9). Call and Put options on VIX or V2X are very popular with trading volumes increasing year after year. The prices of these options is generally expressed in terms of Black-Scholes implied volatility even if the structure of the model is far from log-normality and the index itself is not computed in a linear way. These options will be referred to as options on volatility, whatever the underlying index (VIX or V2X). In order to use Black-Scholes formula, we need to define the options and their associated forward  Table 1 (Set A), initial volatility σ 2 0 is fixed at 0.0079 contract, thanks to (9), we have where the variance time to maturityT := 1 12 is one month in coherence with the definition of VIX and V2X , the expiry of the option is denoted by T . Thus, the underlying volatility swap goes from T to T +T . Figure 5 shows the implied volatility of options written on V2X. The yellow curve is obtained using the calibrated parameters in Table 1 (Set A). We observe that the implied volatility is increasing and then coherent with the usual behavior of VIX implied volatility. Comparing with the exponential law case studied in Nicolato et al. [53], we highlight that the self-exciting property of the -OU Hawkes volatility model changes the slope of VIX implied volatility from negative to positive. Nicolato et al. [53] obtain a positive slope but they need to consider an inverse-gamma law with ν = 1.2 and then a model with infinite variance. The α-Heston model by Jiao et al. [41], which accounts both for power decay and self-exciting features, exhibits an implied VIX volatility which is convex around the money but without any moment larger than the first one. From a financial point of view, the increasing shape of the VIX smile can be explained by the important role played by the options on VIX for high strikes, i.e. protection against turmoils in the equity markets. The strong demand for this protection is responsible for the high levels of the implied volatility for elevated strikes. Figure 5 also details the sensitivities of the implied volatility of options written on VIX with respect to α, β, λ and ρ. First, we note that a change in α impacts the level of implied volatility, acting as a translation of the smile. Without surprise the influence of ρ is the more glaring. Second, we can observe that the VIX implied volatility also exhibits a sensitivity with respect to β and λ for standard perturbations. Given the small sensitivity of implied volatility of the underlying with respect to β and λ, it is natural to choose β and λ to fit the smile of the VIX volatility.  Table 1 (Set A), λ is computed such that σ 2 0 is fixed at 0.0079

Calibration on S&P500 and VIX options
Next, we calibrate the proposed model on S&P500 and VIX options. We find model parameters minimizing the difference in absolute value between the market and model implied volatilities (1 month maturity). Resulting parameters are reported in Table 1, with parameter set B resulting from the calibration on S&P500 options and set C resulting from the calibration on options on VIX. Results show that the -OU Hawkes volatility model calibrates well on S&P500 and VIX options separately. We highlight that the options sets used for the calibration refer to a crisis period, i.e. November 6 2020, for both the COVID-19 pandemic and the close presidential contest in US. However, a joint calibration does not work well. It is easy to remark that the main issue is related to the leverage parameter ρ. The other parameters are similar, in particular comparing with the values estimated out of the crisis (i.e. set A).
In particular we remark that the calibrated leverage ρ is negligible for VIX options. This effect is probably due to the VIX definition. As indicated in the Chicago Board Options Exchange [16] white paper, VIX index is computed as the square root of the static replication of Variance Swap rates. This static replication is introduced in Demeterfi et al. [21] and is based on the crucial hypothesis of continuous paths of the underlying. In Demeterfi et al. [21, pp. 29-35] the impact of jumps on the underlying is detailed showing that a bias exists as long as underlying jumps. In our model, jumps in the underlying disappear as long as the leverage parameter ρ goes to zero. It is then quite natural that a joint S&P500 and VIX calibration is unattainable in practice. Moreover, the VIX is based on an average of short dated options of different expiries and not a single expiry. More surprisingly, the two distinct calibrations on the same day give similar parameters with the substantial exception of the leverage ρ. We could conclude that the -OU Hawkes volatility model fits the behaviour of implied volatilities on both S&P500 and VIX albeit a more detailed theoretical analysis of the jumps corrections is needed to confirm this goodness.

Exact simulation and its performance
We start by pointing out that the model simulation crucially depends on the Hawkes process. Several exact simulation schemes have been proposed in literature for the simulation of this kind of process. Gonzato et al. [29] provides an extensive literature review and finds the exact simulation scheme proposed by Dassios and Zhao [18] is the most efficient, being the fastest among the exact methods. Hence, we decide to simulate the Hawkes process using the exact method of Dassios and Zhao [18]: given an initial date t 0 = 0, a final date T , an initial value for λ 0 and the parameters of the model, the algorithm allows to obtain a sample of the quadruplet where N T is the total number of jumps in the period [0, T ], τ k is the k−th jump time, Z τ k is the k−th jump size which has an exponential distribution with parameter α. Given the quadruplet, we can compute: where Hence, is possible to simulate X T given λ 0 exactly and efficiently (no numerical methods or approximations are required at any step), implying that unbiased estimators for path independent derivatives can be obtained. Moreover, is also possible to extend the methodology in order to generate paths for log-returns process observed at a discrete time grid, we summarize the procedure in Algorithm 1. Compute number of jumps in the interval: A = #{k : t j−1 < τ k ≤ t j }

5:
if A = 0 then 6: compute root mean square errors, RMSE = bias 2 + standard error 2 , where the bias is given by the difference between the true price of the ATM European call option computed using (12) and the simulated one. To estimate the bias precisely and remove the variance, we employ 4 × 10 8 simulation trials for both the Euler and the exact scheme. For the Euler scheme, we consider different numbers of time discretization steps N . Then, using the estimated biases, we compute RMSEs for different number of simulation trials M and record the CPU time necessary to complete the simulation. Results are displayed in Table 2 for the three different parameter sets reported in Table 1. For the exact scheme we don't report any bias estimate since it is unbiased. Numerical results show that the exact scheme outperforms the Euler scheme both in terms of accuracy and CPU time. It presents constantly smaller RMSEs and is faster. This is particularly evident for high number of simulation trials, indeed, in the case with 256 × 10 4 simulations, it is more than 20 times faster for parameters set A and more than 10 times faster for parameters set C than the Euler scheme. Finally, Fig. 7 presents a graphical speed accuracy comparison on a log-log scale. From this plot is possible to see that the proposed method exhibits the optimal convergence rate of 0.5 typical of the exact methods (see e.g. [12]), contrarily to the Euler scheme, whose convergence rate is smaller (around 0.32 for all the parameter sets).

Conclusion
In this concluding section we want to provide a systematic comparison of the properties of the proposed -OU Hawkes volatility model and the two most similar stochastic volatility models available in the literature, i.e. the BNS and Heston models. This comparison can be resumed in the following points.
Ergodic distribution of variance process the three models share the same ergodic distribution, that is of Gamma type. However, a main difference between Heston and -OU Hawkes models is that the volatility process can not reach 0. This property can be easily deduced by the evolution given by (1) with the assumption that λ 0 > λ. Parsimony Heston model is characterized by four parameters, namely long run mean, speed of mean reversion, correlation and volatility of volatility. -OU BNS and Hawkes volatility models are characterised by five parameters: background intensity (playing also the role of a reverting level in the Hawkes model framework), mean reverting speed, leverage coefficient and two parameters for the jump size specifying, respectively, the exponential law intensity and a scaling factor. The scaling factor plays a minor role. Thus, we can could consider that the three models are all parsimonious. Fourier-Laplace transform the three models are exponential affine and the associated pseudo-Riccati differential equation can be solved explicitly. In particular, the inverse of the joint Fourier-Laplace transform of (X , λ) is an explicit function containing only functions of ln and arctan type as in Heston model. It exists a huge amount of literature devoted to numerical methods for option pricing based on Fourier transform inversion methods (e.g. [15,26]). All these methods can be adapted to the present framework of the -OU Hawkes model, see Sect. 4. Moreover, it is not surprising that the term-structure and the evolution in time of the implied volatility for intermediate and longer maturities is similar to the one described by the Heston model. However, for very short maturities, the implied volatility behaviour is definitively more similar to that exhibited by the BNS -OU. Exact simulation scheme we show that the proposed model can be simulated exactly without resorting to any numerical method. This is a significant advantage with respect to the Heston model whose transition can be still simulated exactly but only resorting to time consuming numerical techniques (i.e. repeated numerical inversion of Laplace transforms and root finding algorithms, see [12]). The simulation scheme we propose is unbiased and can be easily exploited not only for European options but also for more complex options such as barrier and other path dependent options. Vanilla option pricing can be performed in a very fast way also by simulating the Hawkes process and then applying a conditional Black-Scholes formula in the same spirit of Hull and White [38], Romano and Touzi [56], Willard [60] and Broadie and Kaya [12]. The same approach can be used also to obtain unbiased estimates for the Greeks. Explosion of moments -OU Hawkes volatility model has not the undesirable property that moments of order higher than 1 can explode in finite time [3]. This feature is important because moment explosion can create problems in some occasions, for example when performing derivatives pricing through Monte Carlo simulation (finite second moment is necessary), or in dynamic portfolio optimization problems with power utility (the value function could be infinite when the expected value of a power of the underlying is infinite). As long as the moment generating function of the jumps size distribution exists, the moments of any order exist for all finite maturities. This is the main difference with the Heston and many other stochastic volatility models (see for instance the inverse-Gamma OU in [53]). From this point of view, -OU Hawkes volatility model is more similar to -OU BNS. Leverage effect it is guaranteed by the negative sign of the parameter ρ multiplying the jump driver as in the usual -OU BNS. In Sect. 4 we show that the slope of ATM implied volatility has a power decay as a function of the time to maturity. Volatility and jumps clusters this is the main difference between our model and the two other models. As a matter of fact, a Hawkes driver implies clusters by construction. It means that the equity prices can experience turmoil periods with spikes and high volatility levels followed by very long periods with a persistence of flat volatility. Options written on VIX -OU Hawkes volatility exhibits an increasing implied volatility for options written on VIX, in coherence with empirical evidence, whereas -OU BNS and Heston model are down-sloping. This result is particularly important since, comparing BNS and Hawkes volatility models, the only difference is in the self-exciting property of the jumps in Hawkes framework. An increasing implied volatility could be obtained in a Lévy framework as for instance inverse-Gamma OU in Nicolato et al. [53], but without preserving the existence of all moments.
In conclusion, the -OU Hawkes volatility model introduced in this paper is parsimonious but can reproduce a broad range of empirical evidences and, in particular, the behavior of VIX options. For numerical purposes, the explicit Laplace transform enables to exploit characteristic function inversion techniques for vanilla pricing and calibration as well. An efficient exact simulation scheme can be exploited for more complex payoffs.