Modelling the clustering of jumps in equity returns: a bifactor self-exciting jump diffusion approach∗ Donatien Hainaut† Univ. Catholique de Louvain

Franck Moraux‡ Univ. Rennes 1 and CREM

October 7, 2017

Abstract It is now well documented that price jumps in stock markets tend to cluster. This paper introduces a new self-exciting jump diffusion process, which can capture many dynamic features of asset returns, among which the jump clustering effect. Here, a shock immediately increases the instantaneous probability of observing a new jump and it has a lasting effect by increasing the average frequency of future jumps in price. The short-term influence can explain why shocks cluster and arrive in groups. The long-term influence can explain why the clustering lasts over time. In our specification, both the occurrence, the size and the sign of past jumps influence the instantaneous probability of observing a jump and its associated long-run mean value. The intensity of jump arrival is a bifactor self-excited process whereas the jump size is a random variable. We show that our specification remains analytically tractable by deriving many analytical results, useful to exploit the proposed framework. Among other things, we derive closed and semi-closed form expressions for the mean and the variance of the intensity as well as for the moment generating of log returns. We find a class of changes of measure that preserves the dynamics of the process under the risk neutral measure and price financial options. We fit the model to the S&P 500 stock index and use our bifactor model to investigate whether only positive, only negative, or both types of price jumps matter in the propagation of shocks. We filter state variables by sequential Monte Carlo and calibrate the model by a peak over threshold approach. We finally evaluate European options and analyze the sensitivity of implied volatilities to parameters and factors. JEL Codes: G10, G13. ∗ Acknowledgements: This paper has been presented at the 2016 EUROFIDAI conference in Paris and various seminars. Many thanks to Elisa Luciano, Alexandre Gohin, Guillaume Bagnarosa, Christos Alexakis , Saqib Aziz ,Michael Dowling, Roman Matkovskyy, Dima Tawil, Dieter Vanwalleghem, Pavel Shevchenko for comments, discussions or suggestions. † Postal address: ISBA, UCL,Voie du Roman Pays 20, B-1348 Louvain-la-Neuve (Belgique). Email to: [email protected]. ‡ The CREM is the Centre de Recherche en Economie et Management (UMR 6211). Postal address: IGR-IAE Rennes, 11 rue Jean Mac´e, 35000 Rennes (France). Email to: [email protected].

1

1

Introduction

It is now well documented that price jumps in stock markets tend to cluster. It is, for instance, very common that stock price shocks arrive in groups during economic downturns. The phenomenon has first been highlighted by Yu (2004) and by Maheu and McCurty (2004) who consider DJIA returns and a list of individual stock returns, respectively. More recently, Ait-Sahalia et al (2015) examine daily returns of international equity market indices and conclude that “Jump clustering in time is a strong effect in the data”. All these investigations evidence that past shocks can affect the probability of observing a new shock in the next future. Moreover, another contribution of this burgeoning literature is that the jumps in price should not be only associated to the arrival of rare and large unexpected news. Actually, frequent and small jumps associated to modest unexpected news may also impact significantly the price dynamics of financial assets. E.g., Maheu and McCurdy (2004) explain that “Market crashes can be realized in a series of jumps over a short period ”. Ait-Sahalia et al (2015) remark that from “mid-September to midNovember 2008, the US stock market jumped by more than 5% on 16 separate days” (see their footnote 2, page 587). In broader terms, the jump clustering dramatically questions the way we understand equity markets and their dynamics. This phenomenon calls for changes in common practices, tools and habits of managers, academics and regulators1 ... With no surprise, a recent strand of the literature makes significant efforts to understand the cause of the clustering of jumps. The why we observe some clustered jumps in equity prices has been especially questioned and investigated with high-frequency data. There is a broad consensus in this literature that news releases impact a lot at both the individual and market levels (see e.g. Lee and Mykland (2008) and Evans (2011))2 . More generally, it is widely believed that the information flow matters (see Rangel (2011), Lee (2012) and Fulop et al. (2015)). The information flow can explain the jump clustering effect because the information process itself tend to be clustered (Maheu and McCurdy (2004)). At this stage, the modelling of the jump clustering remains at an early stage. Some quantitative methods are needed both to model this phenomenon and to investigate consequences of it on asset pricing, option pricing and risk management. As illustrated by previous investigations, self-exciting processes offer a natural solution to model jump clustering in the dynamics of stock prices3 . In such an approach, the 1

E.g., portfolio managers are highly interested in estimating how often unexpected jumps can occur. Lee and Mykland (2008) conclude that the individual stock jumps are linked to company-specific news events (such that scheduled earnings announcements but also unscheduled news), while the “S&P 500 Index jumps are associated with general market news announcement”. More recently, Lee (2012) goes a step further by analyzing the predictability of jumps in individual stock returns, using both macroeconomic and firm-specific news releases. 3 The clustering of jumps is also very frequent in high-frequency data (HFD). As a result, self-exciting processes are commonly used to model HFD. E.g., Engle and Russell (1998) point out the fundamental role of Hawkes self-exciting processes in modelling the duration between two transactions. Readers interested 2

2

occurrence of a shock at a given point in time depends on the history of price jumps. This sort of dynamic was initially introduced by Hawkes (1971a, b) and Hawkes and Oakes (1974) to model earthquake aftershocks. In the most common and simplest specification, the intensity of jumps, akin to the instantaneous probability of a shock, increases as soon as a jump in price occurs. The influence of this jump on the intensity then decays exponentially over time. Hence, these processes may have a shorter or longer term memory of past shocks depending on the value of the ’decay’ parameter. The simplest Hawkes process is nowadays extensively used to explain the daily dynamics of financial assets. However, A¨ıt-Sahalia et al. (2014, 2015) develop a multivariate setting with self and mutual excitations to examine whether some jumps result from past shocks and/or a contagion effect. And Carr and Wu (2016) and Fulop et al. (2015) develop mono-asset settings, where only negative jumps matter in evaluating asset risk. In all the cited papers, the impact of a jump over time is limited because, after a shock, the intensity decreases and reverts to the same fixed and constant level at an exponential rate. The available specifications fail to explain extended periods of jump clustering4 . This article contributes to the literature by proposing a new self-excited model in which jumps in a stock price have also a lasting effect of increasing the frequency of future shocks. Our specification relies on the traditional view that the intensity of jumps instantaneously increases when a shock occurs and then decreases and reverts to a long-run level. But, the novelty we introduce is that this long-run level is no more a constant, as it also depends on past jumps. Hence, the jump arrival intensity reverts to a mean level that is itself increased by jumps in price. We call this model the bifactor self-excited jump diffusion (SEJD) process, because the jumps are led by two dependent state variables: the jump intensity, and its long-run mean level5 . Shocks that occur have therefore an immediate consequence on the instantaneous probability of a shock and a lasting effect. Another distinctive feature of our model is that the different amplitudes and different signs of past shocks can influence differently both the jump intensity and the long-run mean level. As a result, our specification allows us to investigate whether large, small, positive and negative shocks have a similar influence as postulated in most existing research. Most often, it is assumed that jumps whatever their size and their sign, immediately increase the probability of a new shock by a constant amount. A rare exception is Carr and Wu (2016) who consider that only negative shocks modify the probability of future jumps. We by HFD applications may also consult Giot (2005), Bowsher (2007), Large (2007), Chavez-Demoulin and McGill (2012), Bacry et al. (2013), and Da Fonseca and Zaatour (2014) for more recent contributions. For a very interesting overview of these works, one may consult Bauwens and Hautsch (2009). From now on, we keep this issue out of the scope of the present research. It must be noticed that Hawkes’ type processes are now applied to different financial markets. Hainaut (2016) proposes a bivariate Hawkes process for modelling interest rates. Kilic (2017) adapts the multivariate model of A¨ıt-Sahalia et al. (2015) for the exchange markets, or more precisely to study contagion effects between the U.S. dollar and Chinese yuan. 4 Often, the constant long-run intensity is empirically found to be low. Note that a larger value for the constant long-run intensity is not a solution, because it would apply for ever. 5 This bifactor approach should not be confused with the multivariate approach of A¨ıt-Sahalia et al. (2015). There, the multifactor feature refers to the number of financial assets considered simultaneously.

3

use our model to test this hypothesis on the S&P 500 index and investigate whether only negative shocks, only positive shocks or both negative and positive jumps matter. Our specification can consider different distributions for the jump size. Our analysis therefore determines the most appropriate distribution for jumps among the positive exponential distribution, the negative exponential distribution, the double exponential distribution and the double Bernoulli distribution. We undertake a full analysis of our bifactor selfexciting jump diffusion process by developing a number of analytical and quasi-analytical results. The rest of paper proceeds as follows. Section 2 presents the continuous time model. We provide the key properties of this new process and demonstrate, among other things, that our model captures a larger spectrum of dynamics than those of the more common self-exciting processes. Moreover, our model is tractable from a stochastic calculus perspective. We obtain analytical expressions for the first two centered moments of state variables, and semi-closed-form expressions for their moment-generating functions. We show that the expectation and variance of the jump arrival intensity also depend upon a sum of exponential functions. Then we turn to the pricing of some financial and derivative contracts. We identify a family of changes of measure that preserve the structure of our specification, which we then use to define a risk-neutral measure. Here, we provide analytical pricing formulae for realized variances. Section 3 discusses calibration issues, and includes an empirical analysis. To estimate the parameters, we develop an asymmetric peak-over-threshold procedure. Next, a sequential Monte Carlo procedure is used to filter the state variables accurately, and to assess their likelihood. We explore the empirical performance of our model on daily returns of various financial time series collected over a period of 10 years (S&P 500 index). Then, section 4 demonstrates the applicability of the framework to the pricing of derivatives. Lastly, we conclude the paper in section 5.

2 2.1

A new framework for modelling the clustering of jumps The price dynamics

Consider a probability space (Ω, F, P) equipped with a right-continuous filtration {Ft }t≥0 , on which is defined the price process of a stock index, denoted by S = (St )t . In the proposed model, the stock index price is governed in continuous time by a diffusion and a jump process, as follows:   Nt X

dSt = µdt + σdWt + d  eJj − 1  − λt E eJ − 1 dt (1) St− j=1

= µdt + σdWt + eJ − 1 dNt − λt E eJ − 1 dt , where µ is a constant drift term, dWt denotes the increment of the Brownian motion W = (Wt )t , and σ is a constant diffusion coefficient (σ ∈ R+ ). The third term of the 4

above equation highlights the influence of shocks to the price that can occur in the interval (t− , t + dt). The last term is the jump compensator: its presence ensures that the average stock return is equal to µdt. The second equality arises naturally because the probability of observing more than one jump in an infinitesimal period is negligible. The solution to the previous equation is     Z t Nt 2 X
σ St = S0 exp  µ − t − E eJ − 1 λs ds + σWt + Jj  . 2 0 j=1

The process N = (Nt )t counts the number of jumps observed before time t, the sequence of i.i.d. random variables (Jj )j on the size of the return jumps. Each random size (Jj ) is an independent copy of a random variable J. We denote by µJ the expected jump size in absolute terms (i.e. E (|J|) = µJ ), which is, of course, a positive real number. In our numerical applications, we test several distributions for J: positive and negative exponential, double exponential, and double Bernoulli distributions. Double exponential jumps (DEJ)6 are defined by three parameters (p, ρ+ , ρ− ). The features of these parameters are provided in appendix A. Strictly positive or negative jumps are a particular case of the DEN. The intensity of the jump process N is akin to the instantaneous probability of observing a new shock. We assume that the instantaneous probability of jumps i.e. the jump arrival intensity, denoted by λ = (λt )t , depends on past jumps, as follows: ! Nt X (2) dλt = α (θt − λt ) dt + ηd |Ji | i=1

= α (θt − λt ) dt + η |J| dNt . With such a specification, the intensity reverts with speed α ∈ R+ to a long-run mean level θt , which varies over time. When a shock on price occurs, the intensity instantaneously increases by an amount η |J|. This mechanism allows for contagion between jumps, and explains why jumps may arrive in groups in economic downturns. The random variation of intensity is proportional to the absolute value of the jump in price. This is a desirable property, because jumps can be positive or negative, while the jump arrival intensity must remain positive. In A¨ıt-Sahalia et al. (2015), the shock intensity caused by a jump is constant, and depends only on the counting process Nt . Instead, our model links this shock to the amplitude of jumps. If θt were constant and, say, equal to θ, then the solution to equation (2) would simply be Z t −αt λt = θ + (λ0 − θ) e + ηe−α(t−s) dLs (3) 0 6

Assuming a constant jump arrival intensity for the point process N makes our model specification similar to that introduced by Kou (2002). In fact, the latter work introduced the double exponential distribution to model jumps in returns that can be positive or negative. For convenience, the appendix provides some general results for this, now, standard distribution. It is worth emphasizing that the Levy process considered by Kou (2002) is, by design, unable to accommodate the clustering of jumps.

5

P t where L = (Lt )t is a marked point process defined by Lt = N i=1 |Ji | i.e. the sum of the absolute values of jumps in the asset price up to time t. The integrand in the last term of equation (3) is called the kernel function. It captures the decreasing influence of the time s−shock over time. The kernel is here a decreasing exponential function, reflecting that the influence of past jumps on the current intensity decays exponentially. Thus, the influence of jumps is limited in time, and λt is representative of the short-term influence of jumps. Note that while this latter expression may seem familiar to those who know the Hawkes process, this process is not standard at all. Recall that the process L used in this study depends on the realized jumps in absolute terms. Consequently, in the remainder of this paper, we provide results for the general process, as well as for this simpler, but new single-factor process. Figure 1 shows the simulated paths of λt and JdNt in the single-factor model, illustrating that past jumps are forgotten at an exponential rate. Insert Figure 1 Consider now that a shock has a lasting effect of increasing the average frequency of future jumps in price. We assume that the behaviour of the time-varying long-run mean θt is itself mean-reverting, self-exciting, and described by ! Nt X dθt = β (γ − θt ) dt + δd (4) |Ji | . i=1

= β (γ − θt ) dt + δ |J| dNt , where β ∈ R+ is the speed of mean reversion and γ ∈ R+ is the level towards which the long-run mean reverts. A shock in price then increases θt by δ |J|. However, because θt is the mean reversion level of λt , the impact of this shock on the probability of new jumps is more persistent than in the single-factor case (this point is proven mathematically in the next section). A direct integration of equation (4) leads to Z t −βt θt = γ + (θ0 − γ) e + δe−β(t−s) dLs . 0

This solution can be substituted into equation (2) and, in the next section, we derive the analytical expression for the jump arrival intensity when the long-run mean level is stochastic. The system of stochastic differential equations (1), (2), and (4) fully characterizes what we call the bifactor self-exciting jump diffusion (SEJD) model. Here, the two factors are λt and θt . This setting deserves some general comments. First, it is clear that a shock impacts directly and indirectly the equity price process. Second, equations (2) and (4) show that the realized jumps in the underlying price process can simultaneously modify the jump arrival intensity λ = (λt )t and its long-run mean level θ = (θt )t . Nevertheless, η ∈ R+ and δ ∈ R+ are two parameters tuning the influence of past realized 6

jumps on these two processes7 . Third, equations (7) and (9) are both satisfied because the probability to observe more than one jump in an infinitesimal period is almost surely zero.

2.2

Analyzing the jump feature of the bifactor SEJD model

In this section, we provide various results and analytical expressions to aid in the understanding of the main features of the general model. Before entering the core of this section, note that the jump arrival intensity process λ = (λt )t is not a Markov process, but that the joint process (λ, θ, N ) is. Therefore, in order to apply the standard results of stochastic calculus, we must often consider the processes λ = (λt )t , θ = (θt )t , and N = (Nt )t jointly. We first study λ = (λt )t and θ = (θt )t , which are the state R variables. Then, we consider the joint process (λ, θ, N ), which is equivalent to (λ, θ, λs ds). Lastly, we discuss how this applies to any regular function of our state variables. Proposition 1 provides an analytical expression for the value of the jump arrival intensity at time t in the bifactor SEJD model. It shows that the kernel function of the jump arrival intensity is a weighted sum of two exponential functions having different exponential rates. Proposition 1: The jump arrival intensity of the bifactor self-excited process admits the following expression:  Z t α  −βt λt = γ + (λ0 − γ) e−αt + (θ0 − γ) φbi (u) dLu , (5) e − e−αt + α−β 0 where the kernel function is: αδ −β(t−u) φ (u) = e − α−β bi



 αδ − η e−α(t−u) α−β

if α 6= β and λt = γ + (λ0 − γ) e−αt + (θ0 − γ) αte−αt +

Z

t

φ (u) dLu 0

φ (u) = δαe−α(t−u) (t − u) . otherwise. 7

Consequently, we can investigate their respective value and whether they are equal. If these parameters are different from zero, then any realization of a jump immediately increases the instantaneous probability of observing an other jump and the long-run level toward which this intensity process tends to meanrevert. Testing whether η and δ are different from zero is equivalent to testing whether self-excitation matters in these processes. If these parameters are equal, then any realization of a jump impact λ and θ simultaneously and similarly.

7

The kernel function of the bifactor jump process highlighted i Proposition 1 can be rewritten as follows:  α  −β(t−u) φbi (u) = ηe−α(t−u) + δ e − e−α(t−u) α−β # " 1 − e−(α−β)(t−u) −α(t−u) −β(t−u) , for u < t, = ηe + δe α α−β where the functions ηe−α(t−u) and δe−β(t−u) are meaningful. These expressions call for a couple of remarks. First, the bifactor kernel function, φbi (.), nests the function ηe−α(t−u) as a special case (set δ to zero), meaning that the bifactor jump process can capture a larger spectrum of stochastic behaviours than the single-factor jump process can. Second, the
α second term of the r.h.s. (δ α−β e−β(t−u) − e−α(t−u) ) is positive, irrespective of the values of α and β. Thus, assuming that all parameter estimates can be similar, the bifactor kernel function will be at least as large as the single-factor function, meaning that the influence of a past jump on the intensity may last longer in our bifactor setting. Our model may then display a longer influence of shock than that of a basic Hawkes process. Obviously, all parameter estimates cannot be similar because both models aim at capturing the same reality so that parameter estimates will compensate one another. Proposition 2: If (β − (ηµJ − α))2 + 4 (δµJ α + β(ηµJ − α)) ≥ 0, the expectations of θt and λt are given by !
    1 γ1 t − 1 e 0 E (θt |F0 ) γβ −1 γ1
=V (6) V 1 γ2 t − 1 E (λt |F0 ) 0 0 γ2 e     γt θ0 e1 0 −1 V , +V 0 eγ2 t λ0 where γ1 , γ2 are constant, and given by q 1 1 γ1 = ((ηµJ − α) − β) + (β − (ηµJ − α))2 + 4 (δµJ α + β(ηµJ − α)) 2 2q 1 1 γ2 = ((ηµJ − α) − β) − (β − (ηµJ − α))2 + 4 (δµJ α + β(ηµJ − α)), 2 2

(7)

and V is a matrix defined by  V =

−δµJ −δµJ −β − γ1 −β − γ2



the determinant and inverseof which are respectively equal to Υ = δµJ (γ2 − γ1 ) and  −β − γ2 δµJ V −1 = Υ1 . β + γ1 −δµJ An important consequence of this proposition is that the model is stable (in the sense that the limits of λt and θt exist when t → +∞) if and only if γ1 and γ2 are negative. 8

Hereafter, we constrain the structural parameters so as to imply negative values for γ1 and γ2 , and to ensure a stable jump model. We see later that α is also constrained to be larger than ηµJ . Throughout the rest of the paper, we assume that α 6= β. Under these conditions, the above asymptotic expressions become !     − γ11 0 E (θt |F0 ) γβ −1 lim = V V (8) t→∞ E (λt |F0 ) 0 − γ12 0     1) γβ −β γ11γ2 − (γγ21+γ γ2  , =  γβ − (β + γ2 ) (β + γ1 ) δµJ γ11γ2 where γ1 and γ2 are defined in proposition 2.2. By substituting in (β + γ2 ) (β + γ1 ) = −δµJ α and γ1 γ2 = β (α − ηµJ ) − δµJ α (computed in the appendix), we have !   β(α−ηµJ ) γ β(α−ηµ E (θt |F0 ) J )−αδµJ lim = . t→∞ E (λt |F0 ) γ β(α−ηµαβ J )−αδµJ Note that because we constrain γ1 and γ2 to be negative, γ1 γ2 = β (α − ηµJ ) − δµJ α is positive, as is (α − ηµJ ). Therefore, these asymptotic values are well defined. In addition, they share a common term, which means we can write lim (E (θt |F0 )) = lim E (λt |F0 ) −

t→∞

t→∞

βηµJ , β (α − ηµJ ) − αδµJ

so that the asymptotic value for θ is lower than that for λ, becauseµJ is strictly positive. Finally, equation (6) may be rewritten as follows:     γβ 1 (β + γ2 ) (β + γ1 ) + θ0 − (β + γ1 ) δµJ λ0 eγ1 t − E (λt |F0 ) = Υ γ1     1 γβ (β + γ2 ) (β + γ1 ) + θ0 − (β + γ2 ) δµJ λ0 eγ2 t + Υ γ2   1 γβ γβ (β + γ2 ) (β + γ1 ) − Υ γ2 γ1 γ1 t γ2 t : = υ1 e − υ2 e + υ3 , (9) highlighting that the expected value of λt is the difference between two decreasing exponential functions. This expression emphasizes that for certain sets of parameters, the curve of E (λt |F0 ) is increasing, decreasing, and then hump-shaped. As we will see in subsection 2.3, the expected intensities in a single-factor model do not present the same richness of term structure, and are either strictly increasing or decreasing, without the hump. The next proposition provides a general decomposition (in two terms) and an explicit expression of the variance of the jump arrival intensity. The decomposition highlights the 9

key role played by the kernel function in the computation of the variance. Therefore, we expect the variance to be significantly influenced by the shape of the kernel function. Proposition 3: The variance of the jump arrival intensity λt may be written as follows: Z t φbi (u)2 E (λu |F0 ) du, V (λt |F0 ) = E(J 2 ) 0

where Z

t

 υ3  φbi (u)2 E (λu |F0 ) du = ε3 1 − e−(α+β)t + (10) α+β 0    υ   υ1  γ1 t υ2  γ2 t 3 ε1 e − e−2βt − e − e−2βt + 1 − e−2βt + γ1 + 2β γ2 + 2β 2β  

υ3
υ2 υ1 γ1 t −2αt γ2 t −2αt −2αt ε2 e −e − e −e + 1−e + γ1 + 2α γ2 + 2α 2α      υ2 υ1 γ1 t −(α+β)t γ2 t −(α+β)t − , ε3 e −e e −e γ1 + α + β γ2 + α + β  2 2     αδ αδ αδ αδ with ε1 = α−β η − α−β . When the jump , ε2 = η − α−β , and ε3 = 2 α−β size has a double Bernoulli distribution, a simple exponential distribution or a double exponential distribution, its second moment E(J 2 ) admits a closed form (see appendix A). If the process is stable and γ1 and γ2 are negative, the limit of V (λt |F0 ) as t → ∞ is constant as all exponential terms in equation (10) decay to zero: Corollary 1: The asymptotic variance of λt is independent of the initial value of the process, and is equal to γβ (β + γ2 ) (β + γ1 ) lim V (λt |F0 ) = −E(J 2 ) × δµJ γ2 γ1 # "     2  2 αδ αδ 1 αδ 1 αδ η− . 2 + + η− α2 − β 2 α−β α−β 2β α−β 2α

t→∞

The probability of no jump occurring within an interval of time [0, T ] may be computed by considering the following expression:   Z T   P [NT = 0] = E exp − λu du |F0 . 0

In fact, this expression is just a special   case of the joint moment-generating function (mgf) Rt of the triplet Λt ≡ 0 λu du, λt , θt . The next proposition provides a semi-closed-form expression for this mgf. 10


Proposition 4: Define ψ (z1 , z2 ) := E ez1 J+z2 |J| . The joint mgf of Λs , λs , and θs is given by   E eω0 Λs +ω1 λs +ω2 θs | Ft = exp (A(t, s) + B(t, s)λt + C(t, s)θt + ω0 Λt ) , where A(t, s), B(t, s) , and C(t, s) satisfy the following system of ODEs:  ∂A   ∂t = −βγC ∂B = αB − [ψ (0, Bη + Cδ) + ω0 − 1] , ∂t   ∂C = −αB + βC ∂t

(11)

with the terminal conditions A(s, s) = 0, B(s, s) = ω1 , and C(s, s) = ω2 . The closed-form expression of ψ (z1 , z2 ) is provided in the appendix. Closed-form expressions for A(t, s), B(t, s), and C (t, s) are not available, but we can compute them numerically using Euler’s method. The result of Proposition 4 is especially useful for a number of reasons, such as the probability of no jump occurring in certain period (see above), the probability density function of λt or of θt , and so on. Moreover, Proposition 4 is useful, as shown in the following section, to identifying the dynamics of jumps under affine equivalent measures. Before closing this passage on the state variables, it is worth mentioning some results for any stochastic process Y = (Yt )t , defined by Yt = f (t, λt , θt ), where f is a regular function of time, the intensity, and the mean reversion level. Firstly, the infinitesimal generator associated with the underlying jump process is Af

= ft + α(θt − λt ) fλ + β (γ − θt ) fθ Z +∞ +λt [f (t, λt + η |z|, θt + δ |z|) − f (t, λt , θt )] dν(z). −∞

Second, the dynamics of Y can be described by dYt = [ft + α(θt − λt ) fλ + β (γ − θt ) fθ ] dt + [f (t, λt + η |J|, θt + δ |J|) − f (t, λt , θt )] dNt . Coming back to the equity price, the next proposition derives the mgf of the log-return of St , denoted by Xt := ln SS0t . Proposition 5:


If ψ (z1 , z2 ) := E ez1 J+z2 |J| , the mgf of ω1 Xs , for s ≥ t, is given

by E e

ω1 Xs




| Ft =



St S0

ω1 exp (A(t, s) + B(t, s)λt + C(t, s)θt ) ,

where A(t, s), B(t, s), and C(t, s) are solutions to the system of ODEs    ∂A σ2 2 σ2    ∂t = −ω1 µ − 2 − ω1 2 − βγC ∂B = αB + ω1 (ψ (1, 0) − 1) − [ψ (ω1 , B η + C δ) − 1] , ∂t    ∂C = −αB + βC ∂t 11

(12)

with the terminal conditions A(s, s) = 0, B(s, s) = 0, and C(s, s) = 0. The next proposition provides an explicit expression for the realized variance and its expectation, as measured by the quadratic variation of instantaneous returns. Rt s Proposition 6: The realized variance, measured as the quadratic variation of 0 dS Ss over the period [0, T ] is given by Z T
2 2 2 eJu − 1 dNu . σR (T ) = σ T + (13) 0

Its expectation is equal to E

2 σR (T ) | F0




2

=σ T +E




2  e −1 J

Z

T

E (λu | F0 ) du,

(14)

0


2  RT where 0 E (λu | F0 ) du is given by equation (9). Appendix A provides E eJ − 1 for double Bernoulli, simple exponential and double exponential jumps. The previous proposition reveals that the curve of the expected realized variances is a mixture of exponential functions. The expected realized volatility does not admit any semi-closed-form expression. However Brockhaus and Long (2000) suggest evaluating it √ using the following approximation, which is the second-order expansion of x:
q  q 2 (T ) | F
V σ 0 R 2 (T ) | F 2 (T ) | F E σR E σR = 0 0 −
3 , 2 (T ) | F 2 8 E σR 0 where the variance of the realized variance is provided by the next proposition. 2 (T ) is equal to Proposition 7: The variance of σR Z 

4  2 V σR (T ) | F0 = E eJ − 1

T

E (λu | F0 ) du,

0


4  RT where 0 E (λu | F0 ) du is provided by equation (9). Appendix A provides E eJ − 1 for double Bernoulli, simple exponential and double exponential jumps. This result is a direct consequence of the expression (13) for the realized variance.

2.3

The single-factor model

We now examine the difference between our bifactor model, in terms of its short- and long-term influence, and the model where θt is a constant θ. Let’s examine the different expressions. The expected jump arrival intensity is given by   E (λt | F0 ) = λ0 e−(α−ηµJ )t + λ∞ 1 − e−(α−ηµJ )t , (15) 12

α where λ∞ = α−ηµ θ. This expression is well defined and meaningful if α − ηµJ > 0. As J expected, the expected value of the jump arrival intensity tends to a constant as t becomes large, and the long-term value of λt is equal to θ only if the expected size of the jumps is zero. The expression (15) is far simpler than that associated with the bifactor jump process. In this univariate setting, the variance of the jump arrival intensity satisfies Z t 2 φsingle (u)2 E (λu |F0 ) du, V (λt |F0 ) = E(J ) 0

where E (λu |F0 ) is given by equation (15) and φsingle (u) = ηe−α(t−u) . This expression is similar to that in Proposition 3. In addition, after simplification, the variance has the following form:  2   η2λ
η (1 − λ∞ )  −(α−ηµJ )t ∞ 2 −2αt −2αt V (λt |F0 ) = E(J ) e + , −e 1−e α + ηµJ 2α which converges asymptotically to the following constant: η 2 λ∞ . t→∞ 2α When the mean reversion level of λt is constant and equal to θ, the mgf of the log-return Xt is given by
E eω1 Xs | Ft = eω1 Xt exp (A(t, s) + B(t, s)λt ) , s ≥ t, lim V (λt |F0 ) = E(J 2 )

where A(t, s) and B(t, s) are solutions to the following system of ODEs:   ( 2 ∂ σ2 A = −ω µ − − ω12 σ2 − αθB 1 ∂t 2 , ∂ ∂t B = αB + ω1 (ψ (1, 0) − 1) − [ψ (ω1 , B η) − 1]

(16)

with terminal conditions A(s, s) = 0, B(s, s) = 0. Finally, the expression of the realized variance is similar to equation (14), except that the integral of the expected intensity is given by: Z T  λ0 − λ∞  E (λu | F0 ) du = 1 − e−(α−ηµJ )T + λ∞ T. α − ηµJ 0

3

Estimating the bifactor SEJD

In this section, we fit our SEJD models to the S&P 500 index. Our main goal is to test in-sample whether considering a short- and a long-term influence of jumps can help to match the dynamics of the data. The other purpose of this analysis is to select the most appropriate distributions for jump sizes from among the positive, negative, double exponential, and double Bernoulli distributions. Before presenting the empirical results, we describe the data and introduce the econometric method we develop to estimate the structural parameters and to filter the state variables. Note that we only provide a sketch of our econometric approach. The method is described fully in the appendix. 13

3.1

Data description

We collect a sample of S&P 500 daily data from Bloomberg. The sample period ranges from September 2005 to October 2015, and the time series contains 2544 returns. Table 1 provides the relevant statistics. The yearly volatility reaches 20.64% and the very high sample kurtosis indicates that the distribution of returns has fat-tails. The Jarque–Bera and Lillie tests reject normality, whereas the Durbin–Watson statistic reveals serial dependence.

Insert Table 1

Figure 2 plots the prices and returns of the S&P 500 index over the sample period. Grouped jumps in returns are clearly visible from September 2008 to end 2009 (the U.S. credit crunch period), and from September 2011 to February 2012 (the second period of the double-dip recession). Shocks during these periods do not display any clear trend: negative movements alternate regularly with large positive moves. This observation corroborates the link between the frequency of jumps and their absolute values, as assumed in our modelling approach.

Insert Figure 2

3.2

Structural parameter estimation

Estimating the parameters of a jump-diffusion process using a time series is challenging, and requires advanced econometric techniques. Our model involves, in particular, two latent processes (λ = (λt )t and θ = (θt )t ) that we need to back out using a filtering technique. An Euler discretization makes it possible to rewrite our model into a state-space form.8 Nevertheless, the state-space formulation we obtain is highly non-Gaussian and nonlinear, so that we cannot use either the popular Kalman filter, or any variants in the present context. There is a considerable body of literature on how to perform filtering of non-Gaussian and nonlinear state-space models. In this study, we employ a peak-overthreshold approach, similar to that in Embrechts et al. (2011). This approach is robust and sufficient to identify the main advantages of the bifactor model over the single-factor 8

This discretization induces a bias. However, the size of the bias needs to be kept in perspective. Highfrequency models, such as that of Bacry et al. (2013), explain the variation of stock prices exclusively using a bivariate jump process. Over a trading day, the number of jumps is then significantly high, and the price process behaves like a Brownian motion at a macro-level, as proven by Muzy et al. (2013). In our low-frequency approach, the jump process represents only large and unusual shocks, while small jumps arriving at a high frequency are replicated by the diffusion process. Thus, it is not unrealistic to assume that we observe at most one jump per day, at least during periods of normal economic activity.

14

model. The discrete record of T observations of log returns, equally spaced, with a lag ∆ of one day of trading, is {x1, x1 , x2 , ..., xT }. A jump is believed to occur when this return is above or below certain thresholds. These thresholds, denoted by g(α1 ) and g(α2 ) depend on the lag between observations and on the confidence levels, α1 α2 . To define these, we use a log-likelihood maximization to fit a pure Gaussian process: xi ∼ µ∆ + σW∆ . If Φ(.) denotes the probability density function of a standard normal, g(α1 ), g(α2 ) are set √ −1 to the α1 and α2 percentiles of the Brownian term: g(αi ) = σ ∆Φ (αi ). When a jump is detected, the dynamics of the stock index are approximated by: n (xi − µ∆) ∼ Ji (xi − µ∆) > g(α1 ) or (xi − µ∆) < g(α2 ) . The levels of confidence, α1 and α2 , are optimized such that the skewness and the kurtosis of xi for periods without jumps are close to those of a normal distribution. For the S&P index, we find that α1 and α2 are equal to 94% and 91%, respectively. The skewness and kurtosis of the returns for days without detected jumps are equal to 0.047 and 3.28, respectively. The volatilities of the sample from which we draw positive, negative, and both types of jumps are 18%, 16% and 12%, respectively. Once jumps are detected, the sample paths of θt and λt for a given set of parameters are approximated by: ∆θi = β(γ − θi−1 )∆ + δ Ji Ijump at ti ∆λi = α(θi − λi−1 )∆ + η Ji Ijump at ti . Then, the jumps and intensities are estimated by maximizing two log likelihoods: one for the distribution of the jumps, and one for the dynamics of λt . For example, for the DEJD with double exponential jumps, ( P (ρ− , ρ+ , p) = arg max ni=1 log ν (xi | ρ− , ρ+ , p) Ijump at ti P (α, η, β, δ, γ, λ0 ) = arg max ni=1 log P (∆i N | λi ∆) Ijump at ti . After calibrating the parameters, a sequential Monte Carlo procedure (also called particles filter) is applied to estimate the log likelihood and to filter hidden state variables with enhanced accuracy. 3.2.1

Filtering and log-likelihood calculations

As in Doucet et al. (2000), our estimation procedure relies on simulations of a timediscretized version of the SEDJ model, characterized by equations (1), (2), and (4). Denote by ∆ the length of the time interval. The ex ante continuous return (over the period ∆) at time t = j∆, defined by Xj = ln S(j+1)∆ − ln Sj∆ , then satisfies the following equation in discrete time:   √
σ2 0 J Xj = µ − − λj E e − 1 ∆ + σ ∆εj + ∆Lj , (17) 2 | {z } µj

15

PN(j+1)∆ 0 where εj denotes a standard normal random variable and ∆Lj = k=Nj∆ Jk . Here, 9 N(j+1)∆ − Nj∆ is approximated by a Poisson random variable, with parameter λj−1 ∆. The Euler approximation of equations (2) and (4) then provide the discretized dynamics of the latent variables λ = (λt )t and θ = (θt )t , given by

where ∆Lj =

PNj

k=Nj−1

λj = λj−1 +

α(θj−1 − λj−1 )∆ + η ∆Lj

(18)

θj = θj−1 +

β (γ − θj−1 ) ∆ + δ ∆Lj ,

(19)

|Jk |. Note that, at this stage, the model parameters are assumed 0

to be known. Denote vj = (λj , θj , Nj , ∆Lj , ∆Lj ) as the particle that puts together all necessary information about the jump process at time t = j∆. The above system of equations admits a useful state-space representation, where equation (17) provides a measurement equation, or system (the space) that defines the relationship between the (possibly observed) return and the underlying state variables. The 0 vector vj = (λj , θj , Nj , ∆Lj , ∆Lj ) can help in finding the transition system (the state) that describes the dynamics of the state variables. These dynamics depend on equations (18) and (19). Denote the sample of observed continuous returns by {x1 , x2 , ..., xn } and let Gj be the set of data associated with the subset {x1 , x2 , ..., xj }. Then, G = (Gj )j forms a filtration. Conditional contained in vj , the return density p(xj |vj ) is Gaussian  on information √  0 p(xj |vj ) = N µj ∆ − ∆Lj , σ ∆ . On the one hand, it is possible to simulate the transition density p(vj+1 | vj ) using equations (18) and (19). The density of v0 is p(v0 ), and the posterior distribution of vj is denoted by p(vj | y1:j ). Because P (A|B) = P P(A∪B) (B) , this posterior distribution can be rewritten as follows: p(vj | x1:j ) =

p(x1:j , vj ) p(x1:j )

(20)

and, according to the Bayes’ rule, the denominator satisfies the equality p(x1:j ) = p(x1:j−1 , xj ) = p(xj | x1:j−1 )p(x1:j−1 ) On the other hand, the numerator of equation (20) is developed as follows: p(x1:j , vj ) = p(xj | vj )p(vj | x1:j−1 )p(x1:j−1 ). Finally, we obtain the following expression for the posterior distribution: p(vj | x1:j ) = R

p(xj | vj ) p(vj | x1:j−1 ), p(xj |vj )p(vj | x1:j−1 )dvj

9

(21)

Approximating N(j+1)∆ − Nj∆ with a Poisson random variable of parameter λj−1 ∆ introduces a bias in the simulation. An alternative to reduce this bias consists of using the simulation algorithm of Dassios and Hongbiao (2013).

16

where

Z p(vj | x1:j−1 ) =

p(vj | vj−1 )p(vj−1 | x1:j−1 )dvj−1 .

(22)

To summarize, the calculation of p(λj , θj | x1:j ) is done in two steps, similarly to those used in the Kalman Filter. The first is a prediction step, in which we estimate p(vj | x1:j−1 ) by the relation (22). In the correction step, we next calculate the probabilities p(vj |x1:j ) using equation (21). In practice, the integral in the prediction step is replaced by a Monte Carlo 0 (i) (i) (i) (i) (i) simulation of N particles, vj = (λj , θj , ∆Lj , ∆Lj ), for i = 1, . . . , N . In addition, the structure of the particle filter algorithm is as follows: (i)

1. Initial step: Draw N values of v0 , for i = 1, . . . , N , from an initial distribution p(v0 ). 2. For j = 1 : T 0 (i) (i) (i) (i) Prediction step: Draw samples ∆Lj and ∆Lj , and update λj , θj using the relations (i)

λj

(i)

θj

(i)

(i)

(i)

(i)

= λj−1 + α(θj−1 − λj−1 )∆ + η∆Lj   (i) (i) (i) = θj−1 + β γ − θj−1 dt + δ∆Lj . (i)

(i)

Correction step: The particle vj has a probability of occurrence equal to wj =  √  0 (i) (i) P p(xj | vj ) , where p(x | v ) ∼ N µ ∆ − ∆L , σ ∆ . j j j j i=1:N p(xj | vj ) Resampling step: Resample, with replacement, N particles according to the im(i) (i) portance weights wj . The new importance weights are set to wj = N1 . Note that the log-likelihood function yield by this algorithm is obtained by simulations, and is then not smooth enough to estimate the model using log-likelihood maximization directly. Then, a gradient method fails to find a maximum. A way to avoid this drawback consists of using the Particle Markov Chain Monte Carlo (PMCMC) method. However, to limit the computation time, this approach requires accurately determining the prior and transition distributions of parameters, which we do not know.

3.3

Empirical results

The sequential Monte Carlo procedure runs with 5000 particles to evaluate the loglikelihood Akaike Information Criterion (AIC) for the single-factor and bifactor models, with simple or double exponential and double Bernoulli jumps. The results are reported in Table 2. Comparing the AIC statistics suggests that the bifactor model outperforms the single-factor model in fitting the data if jumps have positive exponential distributions. The same conclusion arises by observing larger log-likelihoods for the bifactor model. This suggests that the shocks on the U.S. stock market have a lasting effect, a kind of longrun effect. However, are these likelihoods significantly larger? To answer this question, it is tempting to ‘test’ whether this result is true. Note, however, that these are just 17

average log-likelihoods, which means the test is only approximate. In the final row of Table 2, we routinely apply the likelihood ratio test to assess whether the bifactor model outperforms the single-factor model. The likelihood ratio statistic is commonly defined as LRT = −2 (ln L0 − ln La ), where ln L0 is the log-likelihood of the single-factor specification, and ln La is the log-likelihood of the bifactor specification. This test statistic behaves asymptotically like a random variable with a chi-squared distribution, with degrees of freedom d, where d is the difference between the number of parameters involved in the two specifications. Therefore, we have LRT → χ2d , where d = 2. This analysis leads to the same conclusion. The best fit is obtained with DEJ jumps. This suggests that both positive and negative jumps matter in the U.S. market, and not just negative shocks, as assumed by Carr and Wu (2016). Insert Figure 3 Figure 3 shows the log-likelihood of the model with double exponential jumps, for different sizes of the particles set. We observe that the log-likelihood is stable for the filter with more than 1000 particles. This confirms the reliability of log-likelihood estimates. The particle filter works with few scenarios, but at each time step, the filter retains only the most significant particles through the resampling step. This approach is particularly efficient and thrifty in terms of computations in comparison to a classic Monte Carlo simulation. Insert Table 2 The asymptotic frequency is around 22 jumps per year for models with double exponential jumps. Table 3 reports the parameters for the single-factor and bifactor specifications. Positive shocks occur with a probability of 37%, and their average amplitude is 1 =3.28%. Negative jumps are more frequent, with a slightly lower average amplitude ρ+ 1 ( ρ− =-2.95%). By construction, the parameter η tunes the feedback effect of jumps on λt . The η of the bifactor models are all above the equivalent values in the single-factor model. The speeds of mean reversion of λt are higher in the bifactor models then in the single-factor models to ensure a mean reversion to a long-run stable target. The levels of mean reversion γ and θ are quite similar. Table 4 presents the asymptotic expectation and standard deviations of λt and θt for the single-factor and bifactor models. Insert Table 3

Insert Table 4 18

Figure 4 provides the QQ plots to assess the ability of the DEJ bifactor model to fit the data. The upper, left graph is a standard normality plot used to assess whether residuals are Gaussian. The other graph examines whether the filtered jumps behave according to a double exponential distribution. Insert Figure 4

Insert Figure 5 Figure 5 compares the kernel functions of the bifactor and single-factor models when jumps are double constant and double exponential. The graphs in the second and third rows show the spreads between the double- and single-factor kernels and their ratios, respectively. The kernel function at time zero informs on the instantaneous impact of a price jump on the jump arrival intensity. Furthermore, the graphs reveal that this value is much higher for the bifactor model (η = 381) than it is for the single-factor model (η = 337). As time passes, the bifactor kernel function decreases faster than the singlefactor kernel function does. After 0.08 years, the bifactor kernel function decreases at a lower rate than the single-factor function does, so that the influence of a realized jump in the stock index or the individual index lasts longer in the bifactor model. This is confirmed by ratio of double- to single-factor kernels, which emphasizes the long-term effect. Insert Figure 6 Figure 6 shows the term structure of the expected realized volatility yield by the bifactor model. We observe that the closer λ0 is to the asymptotic level limt→∞ E(λt ), the flatter the curve becomes. As we expect, deviations of λ0 mainly modify the shortterm slope of the curve. Deviations of θ0 from limt→∞ E(λt ) have a bigger impact on the long-run term structure. For values of θ0 above (resp. below) limt→∞ E(θt ), the curve is convex and decreasing (resp. concave and increasing). Insert Figure 7 Using our econometric methodology, we can filter out values of λt and θt over time. We present these time series in Figure 7. Because both λt and θt react to the same jumps, their evolution over time seems similar. They both experience peaks around 2008, during the U.S. credit crunch, and between September 2011 and February 2012, the second period of the double-dip recession. However, because their relative sensitivities to the realized jumps are of different orders, as shown by the estimates of η vs δ in Table 3, we can rationalize several differences. In particular, θt is less affected by jumps than is λt . 19

4

Pricing with the SEDJ model

Empirical evidence suggests that the bifactor SEDJ is suitable for modelling the return time series of the S&P 500 index. This section shows that it is also useful for valuation purposes. Because prices depend on two hidden processes, representative of short- and long-term influence, a market composed of cash and a market composed of stock are incomplete, by construction. Hence, several equivalent measures are eligible as definitions of a risk-neutral measure. In this study, we focus on a family of changes of measure that preserves the dynamics of the process. These are induced by exponential martingales of the form     λt Mt (ξ, ϕ) := exp (κ1 (ξ), κ2 (ξ)) + ξ Lt − κ3 (ξ)t (23) θt   Z Z t 1 t 2 × exp − ϕ ds − ϕ(s)dWs , 2 0 s 0 where (ϕs )s is Fs adapted and ξ is constant. Then, κ1 (ξ), κ2 (ξ), and κ3 (ξ) are functions of ξ. Zhang et al. (2009) use a similar change of measure to simulate rare events of a one-dimensional Hawkes process, without a Brownian component and with constant jumps only. In our framework, jumps are random and the affine change of measure modifies both the frequencies and the distribution. The next proposition details the conditions under which Mt is a local martingale: Proposition 8: If for any parameter ξ, there exist suitable solutions κ1 (.), κ2 (.), and κ3 (.), for the system of equations   =0 κ1 α − κ2 β (24) κ2 βγ − κ3 =0,   κ1 α − (ψ (0, κ1 η + κ2 δ + ξ) − 1) = 0 then Mt (ξ) is a local martingale. Assuming the existence of suitable solutions to the system (24), it is possible to show that Mt is a martingale. An equivalent measure Qξ,ϕ is then defined by: dQξ,ϕ Mt (ξ, ϕ) = . (25) dP Ft M0 (ξ, ϕ) The next proposition develops the dynamics of λt and θt under Q.   Proposition 9: Let NtQ be the counting process associated with the jump arrival t

ξ,ϕ . Then, J Q , the magnitude of the jumps, intensity λQ t = ψ(0, κ1 η + κ2 δ + ξ)λt , under Q is defined by its mgf:   ψ (z , z + (κ η + κ δ + ξ)) Q Q 1 2 1 2 , ψ Q (z1 , z2 ) := EQ ez1 J +z2 |J | = ψ(0, κ1 η + κ2 δ + ξ)

20

PNtQ Q Q Q and (Lt )t , the process such that LQ t = j=1 |Jj |. Then, the dynamics of λt and θt under Qξ are as follows:   Q dλQ = α θ − λ dt + η Q dLQ (26) t t t t
dθtQ = β γ Q − θt dt + δ Q dLQ (27) t , where the parameters verify γ Q = γψ(0, κ1 η + κ2 δ + ξ) η Q = ηψ(0, κ1 η + κ2 δ + ξ) δ Q = δψ(0, κ1 η + κ2 δ + ξ).

The next proposition shows that the jump distribution is preserved under the chosen equivalent measure.   Proposition 10: Under Qξ,ϕ , jumps JiQ are double-exponential random variables, i with density ν Q (z) = pQ ρ+Q e−ρ

+Q z

1{z≥0} − (1 − pQ )ρ−Q e−ρ

−Q z

1{z<0} ,

where the parameters are adjusted as follows: ρ+Q = ρ+ − (κ1 η + κ2 δ + ξ) , ρ−Q = ρ− + (κ1 η + κ2 δ + ξ) , pρ+ ρ−Q pQ = . (pρ+ ρ−Q + (1 − p)ρ− ρ+Q ) Because we have identified a class of equivalent measures, we can define at least one risk-neutral measure under which the expected return on St is equal to the risk-free rate. The dynamics of (Xt )t , under Q, are driven by the following SDE:  
σ2 dXt = µ− − λt E eJ − 1 − ϕt σ dt + σdWtQ + J Q dNtQ 2 !
E eJ − 1 σ2 Q − λt − ϕt σ dt + σdWtQ + J Q dNtQ . = µ− 2 ψ(0, κ1 η + κ2 δ + ξ) Then, because dSt = d(eXt ), we infer that dSt = (µ − ϕt σ) St dt + σ St dWtQ # "
 Q  J −1 E e Q Q dt + St eJt − 1 dNt − λt ψ(0, κ1 η + κ2 δ + ξ) 21

and 

EQ

dSt | Ft St



"



= µ − ϕt σ + λQ EQ e t

JQ

#
E eJ − 1 −1 − dt . ψ(0, κ1 η + κ2 δ + ξ) 

Consequently, the drift is equal to the risk-free rate if and only if     E(eJ −1) Q Q Q J µ + λt E e − 1 − ψ(0,κ1 η+κ2 δ+ξ) − r ϕt = . σ The next corollary summarizes our result. Corollary 2: The Ft adapted process (ϕt )t    Q Q µ + λt EQ eJ − 1 − ϕt = σ

E(eJ −1) ψ(0,κ1 η+κ2 δ+ξ)

 −r

makes the expected instantaneous return of the equity index equal to the risk-free interest rate, and the risk-neutral dynamics of St are described by the following SDE: dSt = rSt dt + σ St dWtQ  Q   i  h Q dt . +St eJt − 1 dNtQ − EQ eJ − 1 | Ft λQ t

(28)

Under the risk-neutral measure, the stock price is ruled by the geometric jump diffusion process presented in equation (28). According to Itˆo’s lemma for semi-martingales, the log-return Xt = log SS0t under Q is then equal to Z Xt = 0

t

Z t Nt  Q   X 1 2 Q J Q e − 1 | Ft λs ds + σ dWtQ + JiQ . r − σ ds − E 2 0 i=1

The remainder of this section focuses on European options written on the stock log-return. The option payoff is denoted by V (T, XT ), where T is the expiry date of the contract. Payoffs of call and put options with a strike price K are given by V (T, XT ) = [S0 eXT −K]+ and V (T, XT ) = [K − S0 eXT ]+ , respectively. The option prices are equal to the product of a discount factor and the integral of the payoff, weighted by the density of XT |Ft . If this density is denoted by fXT |Ft (x), the option value is given by: Q

E



−r(T −t)

e

V (T, XT ) | Ft



= e

−r(T −t)

Z

+∞

−∞

V (T, x) fXT |Ft (x)dx.

(29)

The density function does not admit any closed-form expression, but it can be approximated numerically by inverting the mgf of the log return provided in proposition 5. The next discrete Fourier transform (DFT) is very useful in this case: 22

Proposition 12: Let M be the number of steps used in the DFT, and let ∆x = be the discretization step. We denote ∆z = M2π∆x , and zj

2xmax M −1

= (j − 1)∆z ,

for j = 1 . . . M . The values of fXt (.) at points xk = − M 2 ∆x + (k − 1)∆x are approximated by the sum   M   X 2π 2 Q fXT |Ft (xk ) ≈ Re  Ij ϕ i zj , Xt , λQ (−1)j−1 e−i M (j−1)(k−1)  , (30) t , θt M ∆x j=1

  Q where Ij = 21 1{j=1} + 1{j6=1} , and ϕ i zj , Xt , λQ , θ is the mgf of XT , conditional on Ft . t t Once the density is obtained numerically, the option is evaluated by discretizing the integral in equation (29): M   X V (T, xk ) fXT |Ft (xk .) ∆x ≈ e−r(T −t) EQ e−r(T −t) V (T, XT ) | Ft k=1

Note that this method is different those that proposed by Carr and Madan (1999), who calculate the option value using the DFT directly, and for a series of strike prices. To illustrate this section, we price a European call on the S&P 500. The maturities range from two weeks to three months, and the strike prices run from 90% to 110% of the spot S&P value. We use the parameter estimates provided in Table 3. The intensity and its mean reversion level are set to their asymptotic averages, λ0 = limt→∞ E(λt ) and θ0 = limt→∞ E(θt ), as shown in Table 4. The first graph of Figure 8 shows the surface of implied volatilities. For short-term maturities, we observe a pronounced smile of volatilities, although it is slightly asymmetric. For longer maturities, the smile is flatter and around 25%. The sensitivity of the volatility surface to the parameters is shown in the next three subplots of Figure 8. Increasing γ shifts the volatility surface up by a few percentage points. Reducing the speed of the mean reversion or increasing the feedback parameters η and δ both raise the implied volatilities. The smile is clearly more sensitive to modifications of the dynamics of λt than it is to those related to θt . The two last graphs of Figure 8 illustrate the influence of the initial values λ0 and θ0 on the short- and medium-term smiles. Shifting λ0 up increases the instantaneous probability of observing a series of new jumps, and then increases the implied volatilities for all maturities. On the other hand, increasing θ0 mainly affects the medium-term implied volatilities, and has nearly no impact on the one-month smile. Note that the asymmetry of the smile is controlled by the parameters defining the double exponential jumps. Insert Figure 8

23

5

Conclusion

This article proposes a new model with a short- and a long-term influence of past shocks. The occurrence of a jump in price immediately increases the probability of observing new jumps, and their long-run frequency. We derived closed- and semi-closed-form expressions for the mean and variance of the intensity, for the mgf of the log returns, and for the expected realized variance under various specifications. We showed that the influence is defined by a rich kernel function, made of two exponential functions. We also proved that this influence function will be at least as large as that of a single-factor model. Thus, the influence of past jumps on the intensity may last longer in our setting than it does in a classic self-exciting process. We propose a peak-over-threshold method to the estimate parameters. We also develop a sequential Monte Carlo algorithm to filter hidden state variables and to calculate the log-likelihoods. The model is fitted to the time series of S&P 500 daily returns, and the results corroborate the existence of a long-term influence. In particular, the influence function displays a higher resilience to shocks in the long term than it does in a singlefactor model. We also test whether the stock market remembers only negative shocks, or both positive and negative jumps. Our findings contradict the hypothesis of Carr and Wu (2016), showing that self-excitation is induced by past positive and negative jumps. From an economic point of view, this means that when the stock market is shackled by a movement of large amplitude, investors should anticipate an increase in the risk of immediate jump occurrences, but should also adjust their long-term expectation about the frequency of jumps. Finally, to investigate the pricing capability of our framework, we introduced a class of changes of measure that preserve the dynamics of the process under the risk-neutral measure. Next, we priced different options, and showed that the smile of implied volatilities generated by the model is sensitive to the influence parameters. All in all, we think that our specification represents an interesting compromise between the simplest Hawkes process and some very complicated Hawkes’ type processes. Actually, the standard Hawkes process is not always supported by financial data (see Embrechts et al. (2011), Bacry et al. (2013)). And more intricate processes, which can in principle better match empirical data, have a jump intensity that is, in general, no more than a Markov process. This is, of course, a serious problem if one wants to adopt and deploy a setting with self-exciting feature for financial applications, because in that case we can then no longer use most tools of financial calculus.

24

6 6.1

Appendix Exponential and double-exponential jumps

The double exponential distribution was popularized in the field of finance by the work of Kou (2002). Here, we recall several well-known results. Note that our notation is slightly different from that used by Kou (2002). A double-exponential distributed random variable J may take positive or negative values. Its probability density function (defined on R) is given by + − ν (z) = pρ+ e−ρ z 1{z≥0} − (1 − p) ρ− e−ρ z 1{z<0} , while the associated cumulative distribution function is h  i − + P [J ≤ z] = (1 − p) e−ρ z 1{z≤0} + (1 − p) + p 1 − e−ρ z 1{z>0} . Consequently, this distribution depends on three parameters: ρ+ ∈ R+ , ρ− ∈ R− , and p ∈ (0, 1), where p (resp. (1 − p)) denotes the probability of observing an upward (resp. downward) exponential jump, and ρ1+ (resp. ρ1− ) gives the size of an average positive (resp. negative) jump. When only unidirectional jumps are considered, all developments remain valid with p = 1 or p = 0, for positive and negative exponential jumps. The expected value of the size of jump (J) is the weighted sum of these average sizes; that is, E(J) = p ρ1+ + (1 − p) ρ1− . The expected value of the absolute value of the size of jump (|J|) is E (|J|) = p ρ1+ + (1 − p) |ρ1− | ≡ µJ . Important expressions for our research are the second moments of these two random variables:

2 2 E J 2 = E |J|2 = p + (1 − p) . 2 (ρ+ ) (ρ− )2 Finally, we have ρ− ρ+ + (1 − p) , (31) ρ+ − (z1 + z2 ) ρ− − (z1 − z2 )
+ − if (z1 + z2 ) < ρ+ and (z1 − z2 ) > ρ− , so that E eJ = ψ (1, 0) = p ρ+ρ −1 + (1 − p) ρ−ρ −1 . On the other hand, we have the following useful expressions for the calculation of the expected realized variance and volatility.  +  
2  ρ ρ+ + 1 J E e −1 =p − (32) ρ+ − 2 ρ+ − 1  −  ρ ρ− + 1 + (1 − p) − , ρ− − 2 ρ− − 1   ψ (z1 , z2 ) := E ez1 J+z2 |J| = p

and  E

 ρ+ ρ+ ρ+ 3ρ+ + 1 −4 + +6 + − + ρ+ − 4 ρ −3 ρ −2 ρ −1  −  − ρ ρ ρ− 3ρ− + 1 + (1 − p) −4 − +6 − − − . ρ− − 4 ρ −3 ρ −2 ρ −1


4  e −1 = p J



25

6.2

Double Bernoulli jumps

When J has a double Bernoulli distribution, its probability density function (defined on R) is given by 1 1 ν (z) = p + δ{z=1/ρ+ } + (1 − p) − δ{z=1/ρ− } ρ ρ , where the distribution depends on three parameters: ρ+ ∈ R+ , ρ− ∈ R− , and p ∈ (0, 1), where p (resp. (1 − p)) denotes the probability of observing an upward (resp. downward) constant jump of size ρ1+ (resp. ρ1− ). The expected value of the size of jump (J) is the weighted sum of these average sizes; that is, E(J) = p ρ1+ + (1 − p) ρ1− . The expected value of the absolute value of the size of jump (|J|) is E (|J|) = p ρ1+ + (1 − p) |ρ1− | ≡ µJ , and

E J 2 = E |J|2 = p

1 (ρ+ )2

+ (1 − p)

1

. (ρ− )2

The mgf of the jumps and of their absolute value are given by   (z +z ) 1 (z −z ) 1 ψ (z1 , z2 ) := E ez1 J+z2 |J| = pe 1 2 ρ+ + (1 − p) e 1 2 ρ− ,

(33)

whereas we have the following useful expressions for the calculation of the expected realized variance and volatility:    2  2  1 1
2  J + − ρ ρ E e −1 =p e −1 + (1 − p) e − 1 ,    4  4  1 1
4  J + − ρ ρ E e −1 =p e −1 + (1 − p) e − 1 .

6.3

Proofs

Proof of proposition 1: Integrating equation (2) leads to the following expression for λt : Z t Z t −α(t−s) −αt θs ds + e λ0 + η e−α(t−s) dLs . λt = αe 0

0

Now, insert the expression (4) for θs into this relation. Then, we have   Z t Z s Z t −α(t−s) −βs −β(s−u) −αt λt = αe γ+e (θ0 − γ) + δ e dLu ds+e λ0 +η e−α(t−s) dLs , 0

0

0

or Z λt = γ

t

−α(t−s)

αe

Z ds + (θ0 − γ)

0

+e−αt λ0 + η

t

−α(t−s) −βs

αe

e

t

ds +

e−α(t−s) dLs .

0

26

t

−α(t−s)

Z

αδe 0

0

Z

Z

0

s

e−β(s−u) dLu ds

Then, assuming α 6= β, changing the order of integration, and integrating, gives  
α λt = γ 1 − e−αt + (θ0 − γ) e−βt − e−αt (α − β) Z t   αδ +e−αt λ0 + eβ(u−t) − eα(u−t) + ηe−α(t−u) dLu , 0 α−β as expected. Now, if α = β, then   Z t Z s Z t −αt −α(s−u) −αs −α(t−s) e−α(t−s) dLs , e dLu ds + e λ0 + η γ+e (θ0 − γ) + δ αe λt = 0 0 0 Z t Z t = γ + (λ0 − γ) e−αt + η e−α(t−s) dLs + αt (θ0 − γ) e−αt + δα e−α(t−u) (t − u) dLu . 0

0

 Proof of proposition 2: Because E (dLs |F0 ) = E (|J| |F0 ) × E (λs− |F0 ) , we have from equation (4) that E (θt |F0 ) = γ + e−βt (θ0 − γ) + δµJ

Z

t

e−β(t−s) E (λs− |F0 ) ds.

0

If we derive this last expression with respect to time, we find that E (θt |F0 ) is the solution of an ODE: Z t ∂ E (θt |F0 ) = −βe−βt (θ0 − γ) + δµJ E (λt |F0 ) − βδµJ e−β(t−s) E (λs− |F0 ) ds ∂t 0 = δµJ E (λt |F0 ) − β E (θt |F0 ) + γβ . (34) On the other hand, the expectation of λt is also given by the relation Z t Z t −α(t−s) −αt E (λt |F0 ) = α e E (θs |F0 ) ds + e λ0 + ηµJ e−α(t−s) E (λs− |F0 ) ds. 0

0

Furthermore, if we derive this last expression with respect to time, we obtain the following ODE for E (λt |F0 ): ∂ E (λt |F0 ) = αE (θt |F0 ) − α ∂t

Z

t

αe−α(t−s) E (θs |F0 ) ds − αe−αt λ0 + 0 Z t ηµJ E (λt |F0 ) − αηµJ e−α(t−s) E (λs− |F0 ) ds 0

= αE (θt |F0 ) + (ηµJ − α) E (λt |F0 ) .

27

In matrix form, E (θt |F0 ) and E (λt |F0 ) are solutions to a system of ODEs:  ∂       −β δµJ E (θt |F0 ) γβ ∂t E (θt |F0 ) = + . ∂ α (ηµJ − α) E (λt |F0 ) 0 ∂t E (λt |F0 ) | {z }

(35)

M

Finding a solution requires determining the eigenvalues γ and the eigenvectors (v1 , v2 ) of matrix M , multiplied by the expectations. These eigenvalues cancel the following determinant:   −β − γ δµJ det =0 α (ηµJ − α) − γ and are roots of a second-order polynomial: 0 = (−β − γ) ((ηµJ − α) − γ) − αδµJ = −αδµJ − β(ηµJ − α) + γ [β − (ηµJ − α)] + γ 2 0 = −αδµJ − β(ηµJ − α) + (β − (ηµJ − α)) γ + γ 2 . If the discriminant, ∆, given as ∆ = (β − (ηµJ − α))2 + 4 (δµJ α + β(ηµJ − α)) > 0, is positive, then the roots γ1 and γ2 are given by: 1 γ1 = − (β − (ηµJ − α)) + 2 1 γ2 = − (β − (ηµJ − α)) − 2

1√ ∆, 2 1√ ∆. 2

The product of the roots is the constant of the second-order polynomial γ1 γ2 = −αδµJ − β(ηµJ − α)

(36)

and, by passing, we have (β + γ2 ) (β + γ1 ) h i 1 = [2β + ((ηµJ − α) − β)]2 − (β − (ηµJ − α))2 + 4 (δµJ α + β(ηµJ − α)) 4  1 2 4β + 4β ((ηµJ − α) − β) + ((ηµJ − α) − β)2 − (β − (ηµJ − α))2  = | {z } 4 =0

−δµJ α − β(ηµJ − α) = −δµJ α.

(37)

The eigenvectors of the matrix M are orthogonal, and such that    −β δµJ v1 = 0. v2 α (ηµJ − α) 28

Then, we infer that 

v1i v2i



 =

−δµJ −β − γi

 : i = 1, 2.

Finally, if D = diag(γ1 , γ2 ) and V is the matrix of eigenvectors, the matrix M in equation (35) admits the following decomposition:   −β δµJ = V DV −1 . α (ηµJ − α) If we define



u1 u2

 =V

−1



E (θt |F0 ) E (λt |F0 )

 ,

the system (35) can be rewritten as two independent ODEs, as follows:  ∂       γ1 0 u1 γβ −1 ∂t u1 = . + V ∂ 0 γ2 u2 0 ∂t u2 If we introduce the following notation,     1 γβ −1 = V 0 2 the solutions the ODEs are given by
1 γ1 t e − 1 + d1 eγ1 t γ1
2 γ2 t u2 (t) = e − 1 + d2 eγ2 t , γ2   θ0 0 −1 . Then, where d = (d1 , d2 ) is such that d = V λ0 !
    γt     1 γ1 t − 1 e 0 θ0 e1 0 γβ u1 −1 −1 γ1
V + , V = 1 γ2 t − 1 0 0 eγ2 t λ0 u2 0 γ2 e u1 (t) =

which concludes the proof.  Proof of proposition 3: We define R t the martingale Mt := λt − E(λt |F0 ). Then, [λ, λ]t = [M, M ]t and [M, M ]t = Mt2 −2 0 Ms− dMs . Furthermore, V (λt |F0 ) = E ([M, M ]t |F0 ) = E ([λ, λ]t |F0 ). The expected quadratic variation of λ is given by Z t  αδ  β(u−t) [λ, λ]t = e − eα(u−t) + ηe−α(t−u) dLu , 0 α−β  Z t  αδ  β(u−t) α(u−t) −α(t−u) e −e + ηe dLu 0 α−β t 2 Z t   αδ eβ(u−t) − eα(u−t) + ηe−α(t−u) |Ju |2 dNu . (38) = α−β 0 29


From this last equation, and because E |J|2 = p (ρ+2 )2 + (1 − p) (ρ−2 )2 , we infer that: 2

Z t

V (λt |F0 ) = E(|J| ) 0

2   αδ β(u−t) αδ α(u−t) E (λu |F0 ) du. e e + η− α−β α−β

If we substitute the expression (9) for E (λu |F0 ) into this last equation, this concludes the proof. . Proof of Proposition
4: To simplify the notation, we temporarily denote f = E eω0 Λs +ω1 λs +ω2 θs | Ft as t ≤ s. This is the solution to the following Itˆo equation: 0 = ft + fΛ λt + α(θt − λt ) fλ + β (γ − θt ) fθ Z +∞ +λt f (t, λt + η |z|, θt + δ |z|) − f (.) dν(z) .

(39)

−∞

In the remainder of this proof, f is assumed to be an exponential affine function of λt , θt , and Λt : f = exp (A(t, s) + B(t, s)λt + C(t, s)θt + D(t, s)Λt ) , where A(t, s), B(t, s), C(t, s), and D(t, s) are functions of time, with the terminal conditions A(s, s) = 0, B(s, s) = ω1 , C(s, s) = ω2 , and D(s, s) = ω0 . Under this assumption, the partial derivatives of f are given by:   ∂ ∂ ∂ ∂ ft = A(t, s) + B(t, s)λt + C(t, s)θt + D(t, s)Λt f, ∂t ∂t ∂t ∂t fΛ = D(t, s)f fλ = B(t, s)f : fθ = C(t, s)f. In addition, the integrand in equation (39) becomes: Z +∞ f (t, λt + η |z|, θt + δ |z|) − f (.) dν(z) = f [ψ (0, B(t, s)η + C(t, s)δ) − 1] . 0

Substituting these expressions into equation (39) leads to the following relation:     ∂ ∂ ∂ 0 = A(t, s) + βγC(t, s) + D(t, s)Λt + θt C(t, s) + αB(t, s) − βC(t, s) ∂t ∂t ∂t   ∂ B(t, s) + D(t, s) − αB(t, s) + [ψ (0, B(t, s)η + C(t, s)δ) − 1] . +λt ∂t This relation implies that D(t, s) = ω0 and that A(t, s), B(t, s), and C(t, s) satisfy the system of ODEs:  ∂   ∂t A = −βγC ∂ . (40) ∂t B = αB − [ψ (0, Bη + Cδ) + ω0 − 1]  ∂ ∂t C = −αB + βC 30


Proof of Proposition 5: As usual, denote f = E eω1 Xs | Ft , with t ≤ s. This is a solution to the following Itˆ o equation for a semi-martingale:  
σ2 σ2 J 0 = ft + fX µ − − λt E e − 1 + fXX + α(θt − λt ) fλ + β (γ − θt ) fθ(41) 2 2 Z +∞ f (t, Xt + z, λt + η |z|, θt + δ |z|) − f (.) dν(z) . +λt −∞

If f is assumed to be an exponential affine function of λt , θt , and Xt , then: f = exp (A(t, s) + B(t, s)λt + C(t, s)θt + D(t, s)Xt ) , where A(t, s), B(t, s), C(t, s), and D(t, s) are time dependent functions. The result is proven in a similar way to proposition 4.  Proof of Proposition 8: Let us denote by Yt the exponent of Mt :   λt + ξ Lt − κ3 (ξ)t Yt = (κ1 (ξ), κ2 (ξ)) θt Z Z t 1 t 2 − ϕ(s) ds − ϕ(s)dWs 2 0 0

(42)

Then, the infinitesimal dynamics are given by dYt = κ1 α (θt − λt ) dt + κ2 β (γ − θt ) dt 1 + (κ1 η + κ2 δ + ξ) |J|dNt − κ3 dt − ϕ(t)2 dt − ϕ(t)dWt . 2 In the remainder of this proof, the random measure of J is denoted as Ξ(.), and is such R∞ that J = 0 Ξ(dz). Applying Ito’s lemma for semi-martingales to Mt leads to the next relation: 1 dMt = Mt dYt + Mt d [Yt , Yt ]ct Z ∞ 2  +Mt e(κ1 η+κ2 δ+ξ)|z| − 1 − (κ1 η + κ2 δ + ξ) |z| Ξ(dz)dNt , −∞

which is developed as follows: 1 1 dMt = Mt ((κ1 α − κ2 β) θt + κ2 βγ − κ3 ) dt − ϕ(t)2 Mt dt + Mt ϕ(t)2 dt 2 2   Z ∞  (κ1 η+κ2 δ+ξ)|z| −ϕ(t)Mt dWt − Mt λt κ1 α − e − 1 ν(dz) dt −∞ Z ∞  +Mt e(κ1 η+κ2 δ+ξ)|z| − 1 (Ξ(dz)dNt − λt ν(dz)dt) . −∞

31

Since the integrals with respect to Ξ(dz)dNt − λt ν(dz)dt are local martingales, Mt is also a local martingale if and only if the following relations hold:   =0 κ1 α − κ2 β κ2 βγ − κ3 =0 
R∞  (κ η+κ δ+ξ)|z| 1 2 κ1 α − −∞ e − 1 ν(dz) = 0 .  Proof of Proposition 9: If Yt is the exponent of Mt , as defined by equation (42), and if we denote ψ b = ψ(0, κ1 η + κ2 δ + ξ), the expectation under Q of the mgf of λQ T and θTQ is equal to     Q Q b b = E eYT −Yt +ω1 ψ λT +ω2 ψ θT |Ft EQ eω1 λT +ω2 θT |Ft   b b = e−Yt E eYT +ω1 ψ λT +ω2 ψ θT |Ft .   b b If f (.) denotes E eYT +ω1 ψ λT +ω2 ψ θT |Ft , according to Itˆo’s lemma, this solves the following differential equation: 0 = ft + (κ1 α (θt − λt ) + κ2 β (γ − θt ) − κ3 ) fY + 1 1 α (θt − λt ) fλ + β (γ − θt ) fθ − ϕ2t fY dt + ϕ2 fY Y + 2 2 Z +∞ λt f (t, λt + η|z| , θt + δ|z| , Yt + (κ1 η + κ2 δ + ξ) |z|) − f dν(z) .

(43)

−∞

As usual, we assume that f (.) is an exponential affine function of state variables: f = exp (A(t, T ) + ψ(0, κ1 η + κ2 δ + ξ) (B(t, T )λt + C(t, T )θt ) + D(t, T )Yt ) , in which case,  ft =

 ∂ ∂ ∂ b b ∂ A + ψ λt B + ψ C θt + D Y t f ∂t ∂t ∂t ∂t

fY = D f

fY Y = D 2 f

fλ = Bψ b f

fθ = Cψ b f.

Substituting these expressions into equation (43) leads to the following relation (after grouping terms): ∂ 1 1 A + κ2 βγD − κ3 D + βγψ b C − ϕ2t Df + ϕ2t D2 f ∂t 2 2   Z +∞ h i ∂ b b b Bψ η|z|+Cψ b δ|z|+D(κ1 η+κ2 δ+ξ)|z| +λt ψ B − κ1 αD − αψ B + e − 1 dν(z) ∂t 0     ∂ b ∂ b b +θt ψ C + κ1 αD − κ2 βD + αψ B − βψ C + Yt D . ∂t ∂t 0=

32

Thus, we infer that D(t, s) = 1 as

∂ ∂t D(t, s)

= 0. The, we have that

∂ A + κ2 βγ − κ3 + βγψ b C ∂t Z +∞ h i b b b ∂ b 0=ψ B − κ1 α − αψ B + eBψ η|z|+Cψ δ|z|+(κ1 η+κ2 δ+ξ)|z| − 1 dν(z) ∂t 0 ∂ 0 = ψ b C + κ1 α − κ2 β + αψ b B − βψ b C : . ∂t 0=

Using conditions (24), this system is simplified as follows: ∂ A + βγψ b C ∂t " #
ψ 0, Bψ b η + Cψ b δ + (κ1 η + κ2 δ + ξ) ∂ −1 0 = B − αB + ∂t ψb

0=

0=

∂ C + αB − βC, ∂t

and we can conclude the proof by comparing these results to those of proposition 2.  Proof of Proposition 10: By construction, the mgf for jumps under the risk-neutral measure is the ratio ψ (z, (κ1 η + κ2 δ + ξ)) ψ Q (z, 0) = . ψ (0, κ1 η + κ2 δ + ξ) If we denote ζ = (κ1 η + κ2 δ + ξ), the numerator and denominator in this equation are given by ψ (z, ζ) = ψ (0, ζ) =

pρ+ (ρ− + ζ − z) + (1 − p)ρ− (ρ+ − ζ − z) (ρ+ − ζ − z) (ρ− + ζ − z) + − pρ (ρ + ζ) + (1 − p)ρ− (ρ+ − ζ) . (ρ+ − ζ) (ρ− + ζ)

Then, since Q

ψ (z, 0) =

pρ+ ρ−Q (pρ+ ρ−Q +(1−p)ρ− ρ+Q )

(1−p)ρ− ρ+Q ρ−Q (pρ+ ρ−Q +(1−p)ρ− ρ+Q ) z) (ρ−Q − z)


ρ−Q − z ρ+Q + (ρ+Q −

ρ+Q − z


,

the proof is complete. 

References [1] A¨ıt-Sahalia Y., Cacho-Diaz J., Laeven R.J.A. 2015. Modeling financial contagion using mutually exciting jump processes. Journal of Financial Economics, 117, 585606. 33

[2] A¨ıt-Sahalia Y., Laeven R. J. A., Pelizzon L.. 2014. Mutual excitation in eurozone sovereign CDS. Journal of Econometrics, 183(2):151–167. [3] Bacry E., Delattre S., Hoffmann M., Muzy J. F. 2013. Modelling microstructure noise with mutually exciting point processes. Quantitative Finance, 13(1), 65-67. [4] Bauwens L. and N. Hautsch, 2009: Modelling Financial High Frequency Data Using Point Processes, in: T. Mikosch, J.-P. Kreiß, R. A. Davis, T. G. Andersen (Eds.), Handbook of Financial Time Series, Springer Berlin, 953–979. [5] Bowsher C. 2007. Modelling security market events in continuous time: Intensity based, multivariate point process models, Journal of Econometrics 141(2), 876–912 [6] Brockhaus O., Long D. 2000. Volatility swaps made simple. Risk, Jan., 92-96. [7] Carr P., Madan D. 1999. Option valuation using the fast Fourier transform. Journal of Computational Finance 2 (4), 61–73. [8] Carr P., Wu L. 2016. Leverage Effect, Volatility Feedback, and Self-Exciting Market Disruptions. Forthcoming in Journal of Financial and Quantitative Analysis. [9] Chavez-Demoulin V., McGill J. 2012. High-frequency financial data modelling using Hawkes processes. Journal of Banking and Finance. 36, 3415–3426. [10] Da Fonseca J., Zaatour R. 2014. Hawkes Process: Fast Calibration, Application to Trade Clustering, and Diffusive Limit. Journal of Futures Markets. 34(6), 548-579 [11] Dassios A., Hongbiao Z. 2013. Exact simulation of Hawkes process with exponentially decaying intensity. Electronic Communications in Probability, 18(62), 1-13. [12] Doucet A., De Freitas J.F.G. and Gordon N. . Sequential Monte Carlo Methods in Practice. Cambridge University Press, Cambridge 2000. [13] Embrechts P., Liniger T., Lu L., 2011, Multivariate Hawkes processes: an application to financial data. Journal of Applied Probability, 48 (A), 367–378. [14] Engle R. F., Russell J. R., 1998. Autoregressive Conditional Duration: A New Model for Irregularly Spaced Transaction Data, Econometrica, 66, 1127-1162 [15] Fulop A., Li J., Yu J. 2015. Self-Exciting Jumps, Learning, and Asset Pricing Implications. Review of Financial Studies 28, 876-912. [16] Giot P., 2005. Market risk models for intraday data. European Journal of Finance. 11 (4), 309–324. [17] Hainaut D. 2016. A bivariate Hawkes process for interest rate modelling. Economic modelling (57), 180-196.

34

[18] Hawkes A., 1971(a). Point sprectra of some mutually exciting point processes. Journal of the Royal Statistical Society Series B, 33, 438-443. [19] Hawkes A., 1971(b). Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58, 83–90. [20] Hawkes A., Oakes D., 1974. A cluster representation of a self-exciting process. Journal of Applied Probability, 11, 493-503. [21] Kilic E. 2017. Contagion effects of U.S. Dollar and Chinese Yuan in forward and spot foreign exchange markets. Economic modelling (62), 51-67. [22] Kou S.G., 2002. A Jump-Diffusion Model for Option Pricing.Management Science, 48(8), 1086-1101. [23] Large J., 2007. Measuring the resiliency of an electronic limit order book, Journal of Financial Markets 10(1), 1-25 [24] Maheu J., McCurdy T., 2004. News arrival, jump dynamics and volatility components for individual stock returns. Journal of Finance 59, 755-793 [25] Muzy J-F, Delattre S., Hoffmann M., Bacry E. 2013. Some limit theorems for Hawkes processes and application to financial statistics. Stochastic Processes and their Applications, 123(7), 2475-2499. [26] Yu J. 2004. Empirical characteristic function estimation and its applications. Econometric reviews 23 (2) 93-123. [27] Zhang X.W. , Glynn P.W., Giesecke K., Blanchet J. 2009. Rare event simulation for a generalized Hawkes process. in Proceedings of the 2009 Winter Simulation Conference.

35

7

Figures Figure 1: Simulated sample paths for λt and JdN Simulated sample path of t

18

16

Simulated sample path of JdNt

0.1

t

, t=

JdNt

0.05

14

12

0

10

8

6

-0.05 0

1

2

3

4

5

0

1

years

2

3

4

5

years

This figure presents simulated sample paths for λt and dLt , over five years, in the single-factor model.

Figure 2: Prices and returns of the S&P 500 S&P 500

2500

2000

1500

1000

500 Jan-06

Jan-08

Jan-10

Jan-12

Jan-14

S&P 500 daily return

0.1

0.05

0

-0.05

-0.1 Jan-06

Jan-08

Jan-10

Jan-12

Jan-14

This figure presents S&P 500 daily quotes and log-returns from the 7/9/200

36

Figure 3: Accuracy of the filter Bi-factor model with DEJ 7100

7000

6900

Log-likelihood

6800

6700

6600

6500 0

200

400

600

800

1000

1200

1400

1600

1800

2000

Numbers of particles

Log-likelihood of the bifactor model with double exponential jumps for different sizes of the filter.

Figure 4: QQ plots QQ plot, Gaussian Residuals

QQ plot, Jumps

1

0.9

0.9

0.8

0.8

0.7

0.7

Empirical Quantiles

Empirical Quantiles

1

0.6 0.5 0.4 0.3

0.6 0.5 0.4 0.3

0.2

0.2

0.1

0.1

0

0 0

0.2

0.4

0.6

0.8

1

0

Normal Quantiles

0.2

0.4

0.6

0.8

1

Expo Quantiles

First and second graphs: QQ plots of filtered Gaussian residuals and double exponential jumps, respectively. Note that filtered jumps smaller than 0.4% have been excluded from the sample in order to draw the QQ plot, because they form white noise generated by the particle filter.

37

Figure 5: Kernel functions

Double Exp. J.

400

Bivariate Kernel Univariate Kernel

Bivariate Kernel Univariate Kernel

400

φ(t)

φ(t)

300

Double Cst J.

500

200

300 200

100

-0.4

-0.1

20 0

-0.4

-0.3

-0.2

-0.1

15 10 5 0 -0.5

-0.4

-0.3

-0.2

-0.1

Maturity

-0.1

0

Biv.Kernel - Uni.Kernel

10 0

-0.4

-0.3

-0.2

-0.1

0

Maturity Double Cst J. Biv.Kernel / Uni.Kernel

10

5

0 -0.5

0

-0.2

Maturity Double Cst J.

15

Biv.Kernel / Uni.Kernel

-0.3

20

-10 -0.5

0

Maturity Double Exp. J.

20

-0.4

30

Biv.Kernel - Uni.Kernel

40

-20 -0.5

0 -0.5

0

φ(t) Biv / φ(t) Uni

φ(t) Biv -φ(t) Uni

-0.2

Maturity Double Exp. J.

60

φ(t) Biv / φ(t) Uni

-0.3

φ(t) Biv -φ(t) Uni

0 -0.5

100

-0.4

-0.3

-0.2

-0.1

0

Maturity

The right and left graphs compare the kernel functions of the single-factor and bifactor models, respectively, for double constant and double exponential jumps. The first row shows the kernel functions, but the scale makes a comparison difficult. Thus, the graphs in the second and third rows present the differences and ratios, respectively, of the single-factor and bifactor models.

38

Figure 6: Expected realized volatilities 0.246

0.27

0.244

0.265

θ 0 =6 λ0 =E(λ∞) θ 0 =8 λ0 =E(λ∞) θ 0 =10 λ0 =E(λ∞)

0.242

θ 0 =12 λ0 =E(λ∞)

0.26

0.24

E(σR ) T

E(σR ) T

0.255 0.238 λ0 =4 θ 0 =E(θ ∞)

0.236

λ0 =8 θ 0 =E(θ ∞)

0.234

0.25

0.245

λ0 =12 θ 0 =E(θ ∞) λ0 =16 θ 0 =E(θ ∞)

0.232

0.24

0.235

0.23 0.228

0.23 0

2

4

6

8

10

0

2

Maturity

4

6

8

10

Maturity

Curves of expected realized volatilities for different values. of λ0 and θ0 .

Figure 7: Filtered processes

40

20

0

t

t

200

0 Jan-06

Jan-08

Jan-10

Jan-12

Jan-14

Year 0.1

0.05

JdNt

0

-0.05

-0.1 Jan-06

Jan-08

Jan-10

Jan-12

Jan-14

Year

Filtered values of λt , θt , and jumps for the SEJD model with double exponential jumps.

39

Figure 8: Surface and smiles of implied volatilities

3 months smile 0.35

Volatility

0.8

Volatility

0.6 0.4 0.2 0.3

0.2

0.1

0

100

90

Maturity

110

0.3

0.25

0.2 90

Strike

Initial γ×3 γ×6

95

3 months smile

0.3

0.25

95

100

105

0.3

0.25

0.2 90

110

95

100

105

110

Strike 3 months smile

0.8

0.35 θ 0 =E(θ ∞) λ0 =E(λ∞)

θ 0 =E(θ ∞) λ0 =E(λ∞) θ 0 =E(θ ∞) λ0 = 3 E(λ∞)

0.6

Volatility

Volatility

110

Initial η × 1.5 δ×5

Strike 1 month smile

θ 0 =3 E(θ ∞) λ0 =E(λ∞)

0.4

0.2 90

105

0.35 Initial α × 0.5 β × 0.5

Volatility

Volatility

0.35

0.2 90

100

Strike 3 months smile

95

100

105

110

Strike

θ 0 =E(θ ∞) λ0 = 3 E(λ∞)

0.3

θ 0 =3 E(θ ∞) λ0 =E(λ∞)

0.25

0.2 90

95

100

105

110

Strike

Left upper graph: surface of implied volatilities for European call options on the S&P 500. The other graphs analyse the sensitivity of the volatility smile to changes in the parameters, and to the initial values λ0 and θ0 .

40

8

Tables Table 1: Descriptive Statistics for Returns Data

Mean daily return Standard daily deviation Skewness Kurtosis Jarque–Bera p-value Lillie test p-value Durbin–Watson p-value

Value 0.02% 1.30% -0.33 13.51 1e-3 1e-3 0e-3

The mean, standard deviation, skewness, kurtosis, and statistics of normality and serial dependence for the continuously compounded returns of the S&P 500, from September 2005 to October 2015.

Table 2: Filtered log-likelihood

Log. Lik. AIC

D.Exp. J. 6999 -13979

Log. Lik. AIC LRT p-value

D.Exp. J. 6989 -13962 20 0.00%

Bifactor model Pos. Exp. J. Neg. Exp. J. 6193 6388 -12367 -12756 Single-factor model Pos. Exp. J. Neg. Exp. J. 6193 6384 -12371 -12752 0 8 100% 1.83%

D. Ber. J. 6950 -13879 D. Ber. J. 6944 -13872 12 0.25%

Log-likelihoods, AIC filtered by the particle filter for the double- and single-factor models, and p-values of the log-likelihood ratio test.

41

Table 3: Parameter Estimation

µ σ α η β γ δ p ρ+ ρ− E(J + ) E(J − )

D.Exp. J. 0.05 0.12 18.78 381.80 1.77 5.07 8.37 0.37 30.47 -33.90 3.28% -2.95%

µ σ α η θ p ρ+ ρ− E(J + ) E(J − )

D.Exp. J. 0.05 0.12 14.71 337.08 6.44 0.37 30.47 -33.90 3.28% -2.95%

Bifactor model Pos. Exp. J. Neg. Exp. J. 0.05 0.05 0.18 0.16 13.19 9.91 291.40 220.17 1.62 1.36 1.71 4.28 4.68 3.80 30.47 -33.90 3.28% 0 0 -2.95% Single-factor model Pos. Exp. J. Neg. Exp. J. 0.05 0.05 0.18 0.16 11.46 8.73 273.32 208.82 2.05 4.78 30.47 -33.90 3.28% 0 0 -2.95%

D. Ber. J. 0.05 0.12 17.95 458.06 1.44 4.07 4.07 0.37 30.47 -33.90 3.28% -2.95% D. Ber. J. 0.05 0.128 16.17 436.55 4.88 0.37 30.47 -33.90 3.28% -2.95%

The parameters of the double- and single-factor models are fitted using the POT procedure. The distributions of jumps considered are as follows: double exponential, positive exponential, negative exponential, and double Bernoulli.

42

Table 4: Asymptotic Moments

limt→∞ E(θt ) limt→∞ E(λ p t) limt→∞ V(λt )

D.Exp. J. 8.27 22.03 13.01

limt→∞ E(λ p t) limt→∞ V(λt )

D.Exp. J. 21.80 12.63

Bifactor model Pos. Exp. J. Neg. Exp. J. 2.61 5.62 9.46 16.30 8.22 8.46 Single-factor model Pos. Exp. J. Neg. Exp. J. 9.42 16.18 8.13 8.38

D. Ber. J. 6.41 29.69 12.92 D. Ber. J. 28.61 12.63

Asymptotic values of expected intensities E(θt ), E(λt ), and their standard deviations.

43