1. Introduction

The success of the geometric Brownian motion is directly related to its analytical tractability. Prices of European and most of exotic options are calculable without intensive numerical computations. However, there are many piece of evidences proving that stocks returns are slightly asymmetric and have especially heavier tails than these suggested by a Brownian motion. Furthermore, an analysis of past stocks or commodities prices contradicts the assumption of continuity, inherent to a Brownian motion. Since the eighties, many alternatives have been developed to incorporate asymmetric and/or leptokurtic features in stocks dynamics. In a first category, we find models with an infinite number of jumps, obtained e.g. by subordinating a Brownian motion with an independent increasing Lévy process. This approach has been studied by Madan and Seneta (1990), Madan et al. (1998), Heyde (2000), Barndorff-Nielsen O.E., Shephard (2001) or more recently by Hainaut (2016 b). In a second category, called jump diffusion models, the evolution of prices is driven by a diffusion process, punctuated by jumps at random interval. The two most common jump-diffusion models for stocks are Merton's model with Gaussian jumps (1976) and the double-exponential jump diffusion (DEJD) model, such as presented by Kou (2002) or Lipton (2002). In this last model, the amplitude of jumps is distributed as a doubly exponential random variable. As characteristic functions and Laplace transforms have closed form expressions, Kou and Wang(2003, 2004) priced path dependent options and obtained probabilities of hitting times. Boyarchenko and Levendorskii (2002) and Levendorskii (2004) appraised American, Barrier and Touch-and-out options for the same process, using expected present value operators. Hainaut and Le Courtois (2014), Hainaut (2016 a) studied a switching regime version of the DEJD process, for credit risk applications. Cai and Kou (2011) replaced doubly exponential distributions by mixed exponential jumps. But this model, being over-parameterized, is of limited interest for econometric applications. On another hand, jump diffusion processes are not appropriate to represent commodities. Their prices tend indeed to revert to long run equilibrium prices as illustrated in Bessembinder et al. (1995). Mean reversion is mainly induced by convenience yields. To remedy to this issue, Gibson and Schwartz (1990), Cortazar and Schwartz (1994) and Schwartz (1997) modeled commodities with an Ornstein-Uhlenbeck (OU). Recently, Jaimungal and Surkov (2011) proposed a Levy OU process for modeling energy spot prices and pricing of derivatives.

This work contributes to previous researches in several directions. Firstly, it proposes parsimonious models with and without mean reversion, for stocks and commodities, capable to fit highly leptokurtic processes. To achieve this goal, the return is modeled by a diffusion and a sum of compound Poisson processes. Jumps are Laplace random variables and their frequencies of occurrences are inversely proportional to their average size. This assumption is based on the observation that sparse information has a bigger impact on stocks or commodities prices than regular information. This model, called Continuous Mixed-Laplace Jump Diffusion (CMLJD) duplicates a wide variety of leptokurtic distribution. It is adjustable to time series by likelihood maximization. A second appealing feature of CMLJD is that the number of compound Poisson processes is uncountable. CMLJD is in this sense an extension to continuous time of the Mixed Exponential Jump Diffusion model. The CMLJD converges weakly to a diffusion process punctuated by single jumps, distributed as a continuous mixture of Laplace random variables. In this setting, we infer closed form expressions for the density of jumps, for characteristic functions and moments of log-returns, both for CMLJD with and without mean reversion. Approached formulas of characteristic functions are available and can be used to speed up calculations. The last contribution is empirical. To illustrate its efficiency, the CMLJD model is fitted to four stocks indices (MS World, FTSE, S & P and CAC 40), over a period of 10 years. Whereas the mean reverting CMLJD model is calibrated to four commodities (Copper, Soy Beans, Crude Oil WTI and Wheat), observed over four years. These empirical tests confirm that the CMLJD outperforms the DEJD and Brownian processes. Probability density functions and European options prices are retrieved by a discrete Fourier's transform (DFT). Finally, we study the sensitivity of options implied volatility to parameters.

The rest of the paper is organized as follows. Section 2 introduces the Continuous Mixed-Laplace Jump Diffusion (CMLJD) process and its properties. Section 3 develops the mean reverting CMLJD. Sections 4 and 5 review the DFT methods to compute the probability density functions and options prices. Finally, the section 6 presents an empirical study.

---

2. The Continuous Mixed Laplace Jump Diffusion model

This work extends the mixed double exponential jump diffusion model of Cai and Kou (2011) by considering an uncountable number of jump processes. The construction of this model proceeds with the following steps. Firstly, we present a process for asset log-returns with a finite mixture of Laplace jumps. So as to propose a parsimonious model, parameters are replaced by functions. Secondly, we find the moment generating function of this process when the number of jump processes tends to infinity and show that it converges weakly (or in distribution) toward a jump diffusion process with a single jump component.

The asset price is a process denoted by $\left(A_{t}\right)_{t\geq0}$ and is defined on a probability space $\Omega$, endowed with its natural filtration $\left(\mathcal{F}_{t}\right)_{t\geq0}$ and a probability measure $P$. $P$ is indifferently the real historical measure or a risk neutral measure used for pricing purposes. The log return of $A_{t}$ noted $X_{t}^{n_{k}}$, is such that

$A_{t} = A_{0}\exp\left(X_{t}^{n_{k}}\right),$

(1)

where $n_{k}$ is a parameter that points out the number of jump processes involved in the dynamics of log-return. We assume that $X_{t}^{n_{k}}$ is driven by the following jump-diffusion:

$dX_{t}^{n_{k}} = \mu dt+\sigma dW_{t}+\sum\limits_{k=1:\Delta k:K}dL_{t}^{k}, \ $

(2)

where $\mu\, , \, \sigma$ are respectively the return, and volatility of the Brownian motion $W_{t}$. Whereas $K$ is constant and strictly above one ($K>1)$. The $n_{k}=\frac{K}{\Delta k}$ processes $L_{t}^{k}$, are compound Poisson processes parameterized by $k$, defined as the sum of $N_{t}^{k}$ independent and identically distributed jumps noted $J^{k}$:

$L_{t}^{k}: =\sum\limits_{j=1}^{N_{t}^{k}}J_{j}^{k}\, .$

(3)

The counting processes $N_{t}^{k}$, have intensities equal to $\lambda_{k}\Delta k$ for $k=1:\Delta k:K$. The most popular distributions for jumps are either the Gaussian as in Merton (1976) or the double exponential distributions. However, as emphasized in Kou (2002) or in Kou and Wang(2003, 2004), adding a single double exponential jump process to a diffusion considerably improves the explanatory power of the model. Furthermore, the process remains analytically tractable for options pricing and fits relatively well the surface of implied volatility. From an econometric point of view, the popularity of the double exponential jump diffusion (DEJD) comes from its ability to exhibit reasonable leptokurticity and asymmetry.

Cai and Kou (2011) extend the DEJD by considering a mixture of double exponential jumps and study a dynamics similar to the one of equation (2). But the over-parameterization of this model constitutes a serious drawback. As our purpose is to extend their model with an uncountable number of jump processes to fit processes with a high kurtosis, we remedy to this problem by doing two assumptions. Firstly, jumps $J^{k}$ are Laplace random variables. The Laplace law is a double exponential distribution, with symmetric positive and negative exponential jumps. Secondly, parameters are replaced by functions of the index $k$. The process obtained by this way is parsimonious: it counts the same number of parameters as the DEJD model of Kou (2002). We lose the asymmetry of the Cai and Kou process but our model exhibits a wider range of kurtosis, which is an important driver of option prices. Furthermore, empirical investigations concluding this work emphasizes that our approach outperforms the DEJD model. On the other hand, this specification entails that the jump part in the equation (2) is a martingale. The expectation of $dX_{t}^{n_{k}}$ is equal to the drift, $\mu dt$ and we don't need to introduce a compensator for the jump processes.

More precisely, Laplace densities of $J^{k}$ depend on a parameter $\alpha_{k}$ as follows:

$\mu_{k}(x) = \frac{\alpha_{k}}{2}e^{-\alpha_{k}|x|}\quad for\:k\in[1, K].$

(4)

This is a double exponential distribution for which the probability of observing an upward or downward shock is $\frac{1}{2}$, with respective averages $\frac{1}{\alpha_{k}}$ and $-\frac{1}{\alpha_{k}}$. With a such distribution, the expected jump is null, $\mathbb{E}\left(J^{k}\right)=0$. The characteristic function of $J^{k}$ is also equal to the following quotient:

$\begin{eqnarray} M_{J^{k}}(z) & = & \mathbb{E}\left(e^{izJ^{k}}\right)\nonumber \\ & = & \frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}\quad for\:k\in[1, K]. \end{eqnarray}$

(5)

On the other hand, jump processes $L_{t}^{k}$ have a null mean

$\mathbb{E}\left(L_{t}^{k}|\mathcal{F}_{0}\right) = \lambda_{k}\Delta k\, \mathbb{E}\left(J^{k}|\mathcal{F}_{0}\right)\, t=0\quad for\:k\in[1, K]\, ,$

(6)

and a variance given by:

$\begin{eqnarray} \mathbb{V}\left(L_{t}^{k}|\mathcal{F}_{0}\right) & := & \mathbb{E}\left(\mathbb{E}\left(\left(\sum\limits_{j=1}^{N_{t}^{k}}J_{j}^{k}\right)^{2}|\mathcal{F}_{0}\vee N_{t}^{k}\right)|\mathcal{F}_{0}\right)\nonumber\\ & = & \mathbb{E}\left(N_{t}^{k}\mathbb{E}\left(\left(J^{k}\right)^{2}\right)|\mathcal{F}_{0}\right)=\lambda_{k}\Delta k\, \mathbb{E}\left(\left(J^{k}\right)^{2}\right)\, t. \end{eqnarray}$

(7)

Furthermore $L_{t}^{k}$ are martingales, given that increments of $L_{t}^{k}$ are independent and such that:

$\begin{eqnarray} \mathbb{E}\left(L_{T}^{k}\, |\, \mathcal{F}_{t}\right) & =&L_{t}^{k}+\mathbb{E}\left(L_{T}^{k}-L_{t}^{k}\, |\, \mathcal{F}_{t}\right)\nonumber\\ & =&L_{t}^{k}\, . \end{eqnarray}$

(8)

In order to limit the degrees of freedom, $\alpha_{k}$ and $\lambda_{k}$ are parameterized with the following reasoning. Jumps are related to the flow of information. Good news or bad news, of different importance, arrive according to Poisson processes and prices change in response, according to an exponential jump. If we assume that sparse information has a bigger impact on prices than regular information, intensities $\lambda_{k},$ and average absolute values of jumps $\frac{1}{\alpha_{k}}$, respectively increase and decrease with the index $k$. The next functions for $\alpha_{k}$ and $\lambda_{k}$ satisfy these features:

$\lambda_{k} = \lambda_{0}\, k^{\beta_{1}}\quad\forall k\in[1, K]\quad\beta_{1}>0\, ,$

(9)

$\alpha_{k} = \alpha_{0}\, k^{\beta_{2}}\quad\forall k\in[1, K]\quad\beta_{2}>0\, .$

(10)

where $\lambda_{0}$, $\alpha_{0}$, $\beta_{1}$ and $\beta_{2}$ are positive constants. The positivity of $\beta_{1}$ and $\beta_{2}$ ensures that intensities $\lambda_{k}$ are inversely proportional to average jumps, $\frac{1}{\alpha_{k}}$. Before any further developments, let us define $N_{t}$, a Poisson process with an intensity $\lambda$ equal to

$\begin{eqnarray} \lambda & = & \int_{1}^{K}\lambda_{0}\, k^{\beta_{1}}dk\nonumber \\ & = & \frac{\lambda_{0}}{1+\beta_{1}}K^{\beta_{1}+1}-\frac{\lambda_{0}}{1+\beta_{1}}. \end{eqnarray}$

(11)

Let us denote by $B$ a random variable on the interval $[1, K]$ and defined by the measure $\mu^{B}(k)$ as follows:

$\mu^{B}(k)=\begin{cases} \frac{\lambda_{k}}{\lambda} & k\in[1, K]\\ 0 & k\notin[1, K] \end{cases}.$

(12)

Let $J$ be a random variable distributed as a continuous mixture of Laplace random variables:

$J = \int_{1}^{K}J^{k}d\mathbf{\delta}_{\{B=k\}}.$

(13)

Then using nested conditional expectations, the characteristic function of $J$ is such that:

$\begin{eqnarray} \mathbb{E}\left(e^{izJ}\right) & = & \mathbb{E}\left(\mathbb{E}\left(e^{izJ^{B}}\right)\, |\, B\right)\nonumber \\ & = & \int_{1}^{K}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}d\mu^{B}(k)\nonumber \\ & = & \int_{1}^{K}\frac{\lambda_{k}}{\lambda}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk\, . \end{eqnarray}$

(14)

When the number of jumps components tends to infinity, $n_{k}\to\infty$, the process $X_{t}^{n_{k}}$ converges in distribution (weak convergence) toward $X_{t}$ that is a jump diffusion process. As stated in the next proposition, the jump component of $X_{t}$ is a unique compound Poisson process, with jumps distributed as a finite mixture of Laplace's random variables.

Proposition 2.1. $X_{t}^{n_{k}}$ converges in distribution toward $X_{t}$, $X_{t}^{n_{k}}\stackrel{d}{\to}X_{t}$, which is a process defined by:

$dX_{t} = \mu dt+\sigma dW_{t}+dL_{t}\, ,$

(15)

where $L_{t}:=\sum_{k=1}^{N_{t}}J$ is a compound Poisson process. $J$ is defined by equation (13) and $N_{t}$ is a point process with an intensity given by equation (11). The characteristic function of $X_{t}$ is equal to:

$\begin{eqnarray} M_{X_{t}}(z) & = & \mathbb{E}\left(e^{izX_{t}}\, |\, \mathcal{F}_{0}\right)\nonumber \\ & = & \exp\left(t\left(\mu iz-\frac{1}{2}\sigma^{2}z^{2}-\lambda\left(1-\int_{\text{1}}^{K}\frac{\lambda_{k}}{\lambda}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk\right)\right)\right)\, . \end{eqnarray}$

(16)

Proof. To prove this result, we show that characteristic functions of jump processes (2) converge to the one of equation (15). The $L_{t}^{k}$ have a characteristic function equal to (for a proof see e.g. Kaas et al. 2008, page 43),

$M_{L_{t}^{k}}(z)=\mathbb{E}(e^{ziL_{t}^{k}})=m_{N_{t}^{k}}(\ln M_{j^{k}}(z))\qquad for\:k=1:\, \Delta k\, :K.$

(17)

where $m_{N_{t}^{k}}(h)$ is the moment generating function of $N_{t}^{k}$, $m_{N_{t}^{k}}(h)=\mathbb{E}\left(e^{hN_{t}^{k}}\right)=e^{-(\lambda_{k}\Delta k)t(1-e^{h})}$ and $M_{j^{k}}(z)$ is the characteristic function of $J^{k}$, such as defined by equation (5). Then, $M_{L_{t}^{k}}(z)$ is equal to:

$M_{L_{t}^{k}}(z) = \exp\left(-\left(\lambda_{k}\Delta k\right)\, t\, (1-\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}})\right)\qquad for\:k=1:\, \Delta k\, :K.$

(18)

Given that $L_{t}^{k}$'s are independent, the sum of all jumps components till time $t$, $L_{t}:=\lim_{\Delta k\to0}\sum_{k=1:\Delta k:K}L_{t}^{k}$, has the following characteristic function:

$\begin{eqnarray} M_{L_{t}}(z) & = & \lim\limits_{\Delta k\to0}\mathbb{E}(e^{zi\sum\limits_{k=1:\Delta k:K}L_{t}^{k}})\nonumber\\ & = & \lim\limits_{\Delta k\to0}\prod\limits_{k=1:\Delta k:K}\mathbb{E}(e^{ziL_{t}^{k}})\, . \end{eqnarray}$

(19)

Wherein, the product in this limit is equal to:

$\prod\limits_{k=1:\Delta k:K}\mathbb{E}(e^{ziL_{t}^{k}}) = \exp\left(-\sum\limits_{k=1:\Delta k:K}\left(\lambda_{k}\Delta k\right)t(1-\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}})\right).$

(20)

If $\lambda_{\Delta k}=\sum_{k=1:\Delta k:K}\left(\lambda_{k}\Delta k\right)$, this characteristic function becomes:

$\prod\limits_{k=1:\Delta k:K}\mathbb{E}(e^{ziL_{t}^{k}})=\exp\left(-\lambda_{\Delta k}t(1-\sum\limits_{k=1:\Delta k:K}\frac{\lambda_{k}}{\lambda_{\Delta k}}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}\Delta k)\right)\nonumber \\ \qquad=\exp\left(-\lambda_{\Delta k}t\left(1-\mathbb{E}\left(\exp\left(iz\sum\limits_{k=1:\Delta k:K}I_{\{B_{\Delta k}=k\}}J^{k}\right)\, |\, B_{\Delta k}\right)\right)\right),$

(21)

where $B_{\Delta k}$ denotes here a random variable defined on the interval $[1, K]$ by a discrete measure $\mu^{B_{\Delta k}}(k)$:

$\mu^{B_{\Delta k}}(k)=\begin{cases} \frac{\lambda_{k}\Delta k}{\lambda_{\Delta k}} & \qquad k=1:\Delta k:K\\ 0 & \qquad else \end{cases}.$

(22)

$\prod_{k=1:\Delta k:K}\mathbb{E}(e^{ziL_{t}^{k, \Delta k}})$ is then the characteristic function of a jump process, of intensity $\lambda_{\Delta k}$, with jumps distributed as a finite mixture of Laplace random variables. Taking the limit of (21) when $\Delta k\to0$, and according to the definition of $\lambda$, we get that

$\begin{eqnarray} \lim\limits_{\Delta k\to0}\sum\limits_{k=1:\Delta k:K}\left(\lambda_{k}\Delta k\right) & = & \int_{1}^{K}\lambda_{k}dk\nonumber\\ & = & \lambda\:. \end{eqnarray}$

(23)

On another hand, we have that

$\lim\limits_{\Delta k\to0}\sum\limits_{k=1:\Delta k:K}\frac{\lambda_{k}}{\lambda_{\Delta k}}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}\Delta k = \int_{1}^{K}\frac{\lambda_{k}}{\lambda}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk.$

(24)

The characteristic function of $L_{t}$ is then equal to:

$\begin{eqnarray} M_{L_{t}}(z) & = & \exp\left(-\lambda t(1-\int_{1}^{K}\frac{\lambda_{k}}{\lambda}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk\right)\label{eq:MStProof}\\ & = & \exp\left(-\lambda t\left(1-\mathbb{E}\left(\exp\left(i\, z\, \int_{1}^{K}J^{k}d\mathbf{\delta}_{\{B=k\}}\right)\, |\, B\right)\right)\right). \end{eqnarray}$

(25)

As there is an unequivocal correspondence between the moment generating function of a random variable and its probability density function, we have proven that

$\lim\limits_{\Delta k\to0}P\left(\sum\limits_{k=1:\Delta k:K}L_{t}^{k}\leq x\right)=P\left(L_{t}\leq x\right)\, .$

(26)

$X_{t}^{n_{k}}$ converges then weakly or in distribution toward $X_{t}$.

This proposition reveals an interesting feature of our model and shared with all Mixed Exponential Models. Whatever the number of jumps components, the dynamics of the asset return always converges in a weak sense toward a jump diffusion model, with a single compound Poisson process for which jumps are distributed as a mixture of distributions.

The next proposition presents a closed form expression for the density of the mixture of Laplace jumps.

Proposition 2.2. Let us define a constant $\gamma$ by:

$\gamma = 1/\beta_{2}\left(1+\beta_{1}+\beta_{2}\right),$

(27)

the probability density function of jumps $J$ defined in equation (13) is given by the following expression:

$\mu^{J}(x)=\frac{\lambda_{0}}{\lambda}\, \frac{\alpha_{0}}{2}\frac{1}{\beta_{2}}\frac{1}{\left(\alpha_{0}|x|\right)^{\gamma}}\left(\Gamma\left(\gamma, \alpha_{0}\, |x|\right)-\Gamma\left(\gamma, \alpha_{0}K^{\beta_{2}}\, |x|\right)\right)\, ,$

(28)

where $\Gamma\left(a, x\right)$ is the incomplete Gamma function defined as:

$\Gamma\left(a, x\right) = \int_{x}^{+\infty}u^{a-1}e^{-u}du\, .$

(29)

Proof. By construction, the probability density function of jumps is equal to

$\begin{eqnarray} \mu^{J}(x) & = & \int_{1}^{K}\frac{\alpha_{k}}{2}e^{-\alpha_{k}|x|}\mu^{B}(dk)\nonumber \\ & = & \frac{1}{\lambda}\int_{1}^{K}\lambda_{k}\frac{\alpha_{k}}{2}e^{-\alpha_{k}|x|}dk\nonumber \\ & = & \frac{\lambda_{0}}{\lambda}\, \frac{\alpha_{0}}{2}\int_{1}^{K}k^{\left(\beta_{1}+\beta_{2}\right)}e^{-\alpha_{0}\, k^{\beta_{2}}|x|}dk\, . \end{eqnarray}$

(30)

Substituting $k'=\alpha_{0}\, k^{\beta_{2}}|x|$ to the integration variable $k$ leads to the following relations

$k = \frac{1}{\left(\alpha_{0}|x|\right)^{1/\beta_{2}}}\left(k'\right)^{\frac{1}{\beta_{2}}}\, ,$

(31)

$dk = \frac{1}{\beta_{2}\left(\alpha_{0}|x|\right)^{1/\beta_{2}}}\left(k'\right)^{\frac{1}{\beta_{2}}-1}dk'\, ,$

(32)

and to the next expression for the density:

$\mu^{J}(x) = \frac{\lambda_{0}}{\lambda}\, \frac{\alpha_{0}}{2}\frac{1}{\beta_{2}}\frac{1}{\left(\alpha_{0}|x|\right)^{1/\beta_{2}\left(1+\beta_{1}+\beta_{2}\right)}}\int_{\alpha_{0}\, |x|}^{\alpha_{0}K^{\beta_{2}}\, |x|}\left(k'\right)^{\left(\frac{1}{\beta_{2}}\left(\beta_{1}+1\right)+1\right)-1}e^{-k'}dk'.$

(33)

The integral in this last equation is the difference between two incomplete Gamma functions, such as defined by equation (29).

Remark that $X_{t}$, being a jump diffusion process, belongs to the family of Lévy processes. Its infinitesimal generator is then equal to

$\left(\mathcal{L}u\right)(x) = \mu\frac{\partial}{\partial x}u(x)+\frac{\sigma^{2}}{2}\frac{\partial}{\partial x^{2}}u(x)+\lambda\int_{-\infty}^{+\infty}\left(u(x+y)-u(x)\right)\mu^{J}(y)dy\, ,$

(34)

for any function $u(x)$ that is twice continuously differentiable and where $\mu^{J}(.)$ is given by (28). This generator is the key used later, to build the Feynman-Kac equation, satisfied by option prices. This equation is solved numerically by inverting the Fourier transform of option prices. But this approach requires to know the characteristic function of $X_{t}$. The next proposition provides us this important result:

Proposition 2.3. The characteristic function $M_{X_{t}}(z)=\mathbb{E}\left(e^{izX_{t}}\right)$ of the CMLJD process $X_{t}$ is given by

$\begin{eqnarray} M_{X_{t}}(z) & := & \exp\left(t\psi(z)\right)\label{eq:characfuncXeq1}\\ & = & \exp\left(t\left(i\, \mu\, z-\frac{1}{2}\sigma^{2}z^{2}+\int_{1}^{K}\lambda_{k}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk-\lambda\right)\right)\, , \end{eqnarray}$

(35)

and if we denote by $\theta=2\beta_{2}+\beta_{1}+1$, the integral in (35) is equal to

$\int_{1}^{K}\lambda_{k}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk=\nonumber \\ \qquad\lambda_{0}\left(\frac{\alpha_{0}}{z}\right)^{2}\frac{1}{\theta}\left[K^{\theta}G\left(\frac{\theta}{2\beta_{2}}, -\left(\frac{\alpha_{0}}{z}\right)^{2}K^{2\beta_{2}}\right)-G\left(\frac{\theta}{2\beta_{2}}, -\left(\frac{\alpha_{0}}{z}\right)^{2}\right)\right]\, ,$

(36)

where $G(b, x)$ is the hypergeometric function:

$\begin{eqnarray} G(b, x) & = & _{2}F_{1}(1, b, b+1, x)\nonumber \\ & = & b\int_{0}^{1}\frac{u^{b-1}}{1-ux}du\, . \end{eqnarray}$

(37)

Proof. As mentioned in the proof of proposition 2.1, the characteristic function of $S_{t}=\sum_{i=1}^{N_{t}}J_{i}$ is

$M_{S_{t}}(z)=e^{-\lambda t(1-M_{J}(z))}\, ,$

(38)

where $M_{J}(z)$ is provided by equation (14). Then $M_{X_{t}}(z)=\mathbb{E}\left(e^{izX_{t}}\right)$ is equal to expression (35). The integral present in equation (35) is split as follows:

$\begin{eqnarray} \int_{1}^{K}\lambda_{k}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z{}^{2}}dk & = & \lambda_{0}\int_{1}^{K}\, \frac{k^{2\beta_{2}+\beta_{1}}}{k^{2\beta_{2}}+\left(\frac{z}{\alpha_{0}}\right)^{2}}dk\nonumber\\ & = & \lambda_{0}\left(\int_{0}^{K}\, \frac{k^{2\beta_{2}+\beta_{1}}}{k^{2\beta_{2}}+\left(\frac{z}{\alpha_{0}}\right)^{2}}dk-\int_{0}^{1}\, \frac{k^{2\beta_{2}+\beta_{1}}}{k^{2\beta_{2}}+\left(\frac{z}{\alpha_{0}}\right)^{2}}dk\right). \end{eqnarray}$

(39)

To calculate the integral $\int_{0}^{s}\, \frac{k^{2\beta_{2}+\beta_{1}}}{k^{2\beta_{2}}+\left(\frac{z}{\alpha_{0}}\right)^{2}}dk$ with $s=K$ or $s=1$, the next change of variable is done: $k=u^{\frac{1}{2\beta_{2}}}s$. As $dk=\frac{1}{2\beta_{2}}\, s\, u^{\frac{1}{2\beta_{2}}-1}du$, we infer that:

$\begin{eqnarray} \int_{0}^{s}\, \frac{k^{2\beta_{2}+\beta_{1}}}{k^{2\beta_{2}}+\left(\frac{z}{\alpha_{0}}\right)^{2}}dk & = & \int_{0}^{s}\, \frac{u^{\frac{1}{2\beta_{2}}\left(2\beta_{2}+\beta_{1}\right)}s^{\left(2\beta_{2}+\beta_{1}\right)}}{us^{2\beta_{2}}+\left(\frac{z}{\alpha_{0}}\right)^{2}}\frac{1}{2\beta_{2}}\, s\, u^{\frac{1}{2\beta_{2}}-1}du\nonumber\\ & = & \frac{1}{2\beta_{2}}s^{\left(2\beta_{2}+\beta_{1}+1\right)}\left(\frac{\alpha_{0}}{z}\right)^{2}\int_{0}^{s}\, \frac{u^{\frac{1}{2\beta_{2}}\left(\beta_{1}+1\right)}}{1-u\left(-s^{2\beta_{2}}\left(\frac{\alpha_{0}}{z}\right)^{2}\right)}du\nonumber\\ & = & \frac{1}{\theta}s^{\theta}\left(\frac{\alpha_{0}}{z}\right)^{2}G(\frac{\theta}{2\beta_{2}}, -\left(\frac{\alpha_{0}}{z}\right)^{2}s^{2\beta_{2}})\, . \end{eqnarray}$

(40)

Given that the hypergeometric function, $_{2}F_{1}(a, b, c, x)$, is defined by

$_{2}F_{1}(a, b, c, x) = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_{0}^{1}u^{b-1}(1-u)^{c-b-1}(1-ux)^{-a}du\, ,$

(41)

we infer that $G(b, x)=_{2}F_{1}(1, b, b+1, x)$ and conclude.

The moments of $X_{t}$ are obtained by differentiating the characteristic function, as stated in the next proposition. The skewness is null as by construction the distribution of $X_{t}$ is symmetric. But the kurtosis is always above 3. $X_{t}$ has then fatter tails than a normal distribution.

Proposition 2.4. The mean, variance, skewness and kurtosis of $X_{t}$ are respectively given by:

$\mathbb{E}(X_{t}) = \mu\, t\, ,$

(42)

$\mathbb{V}(X_{t}) = t\left(\sigma^{2}+2\frac{\lambda_{0}}{\alpha_{0}^{2}}\frac{1}{\beta_{1}-2\beta_{2}+1}\left(K^{\beta_{1}-2\beta_{2}+1}-1\right)\right)\, ,$

(43)

$\mathbb{S}(X_{t})=0\, ,$

(44)

$\mathbb{K}(X_{t}) =3+\frac{\frac{24\lambda_{0}}{\alpha_{0}^{4}\left(\beta_{1}-4\beta_{2}+1\right)}\left(K^{\beta_{1}-4\beta_{2}+1}-1\right)}{t\left(\sigma^{2}+2\frac{\lambda_{0}}{\alpha_{0}^{2}}\frac{1}{\beta_{1}-2\beta_{2}+1}\left(K^{\beta_{1}-2\beta_{2}+1}-1\right)\right)^{2}}\, .$

(45)

Proof. The moments of $X_{t}$ are obtained by deriving the characteristic function with respect to $z$,

$\mathbb{E}(X_{t}^{k}) = \left.\frac{\partial^{k}}{\partial z^{k}}M_{X_{t}}(-iz)\right|_{z=0}\, .$

(46)

In particular,

$\mathbb{E}(X_{t}) = \mu\, t\, ,$

(47)

$\mathbb{E}(X_{t}^{2}) = \left(\mu\, t\right)^{2}+t\left(\sigma^{2}+\int_{1}^{K}2\, \frac{\lambda_{k}}{\alpha_{k}^{2}}\, dk\right)\, ,$

(48)

$\mathbb{E}(X_{t}^{3}) = \left(\mu\, t\right)^{3}+3t^{2}\mu\left(\sigma^{2}+\int_{1}^{K}2\, \frac{\lambda_{k}}{\alpha_{k}^{2}}\, dk\right)\, ,$

(49)

$\mathbb{E}(X_{t}^{4}) = \left(\mu\, t\right)^{4}+6t^{3}\mu^{2}\left(\sigma^{2}+\int_{1}^{K}2\frac{\lambda_{k}}{\alpha_{k}^{2}}\, dk\right)+\nonumber\\ \quad \quad 3t^{2}\left(\sigma^{2}+\int_{1}^{K}2\frac{\lambda_{k}}{\alpha_{k}^{2}}\, dk\right)^{2}+t\left(\int_{1}^{K}24\frac{\lambda_{k}}{\alpha_{k}^{4}}\, dk\right)\, .$

(50)

The skewness and kurtosis are inferred from following relations:

$\mathbb{S}(X_{t}) = \frac{\mathbb{E}(X_{t}^{3})-3\mathbb{E}(X_{t})\mathbb{V}(X_{t})-\mathbb{E}(X_{t})^{3}}{\mathbb{V}(X_{t})^{\frac{3}{2}}}\, ,$

(51)

$\mathbb{K}(X_{t}) =\frac{1}{\left(\mathbb{V}(X_{t})\right)^{2}}\left(\mathbb{E}(X_{t}^{4})-4\mathbb{E}(X_{t})\mathbb{E}(X_{t}^{3})+6\mathbb{E}(X_{t})^{2}\mathbb{E}(X_{t}^{2})-3\mathbb{E}(X_{t})^{4}\right)\, .$

(52)

A helpful feature of the hypergeometric function for numerical purposes, is that it can be rewritten as an infinite sum. In this case, the characteristic exponent admits the following alternative expression:

Corollary 2.5. The characteristic exponent $\psi(z)$ is equal to the sum:

$\psi(z) = i\, \mu\, z-\frac{1}{2}\sigma^{2}z^{2}-\lambda+\sum\limits_{j=0}^{\infty}\frac{\lambda_{0}}{(\theta+2\, j\, \beta_{2})}(-1)^{j}\left(K^{2j\beta_{2}+\theta}-1\right)\left(\frac{\alpha_{0}}{z}\right)^{2(j+1)}$

(53)

where $\theta=2\beta_{2}+\beta_{1}+1$.

Proof. $_{2}F_{1}(a, b, c, x)$ has the property to be equivalent to the infinite series:

$\begin{eqnarray} _{2}F_{1}(a, b, c, x) & = & 1+\frac{a\, b}{1!\, c}x+\frac{a(a+1)\, b(b+1)}{2!\, c(c+1)}x^{2}\nonumber \\ & & +\frac{a(a+1)(a+2)\, b(b+1)(b+2)}{3!\, c(c+1)(c+2)}x^{3}+\ldots \end{eqnarray}$

(54)

This feature allows us to develop $G(b, x)$ as follows:

$G(b, x) = 1+\frac{b}{(b+1)}x+\frac{\, b}{\, (b+2)}x^{2}+\frac{b}{(b+3)}x^{3}+\ldots$

(55)

and the difference present in the characteristic exponent, becomes

$K^{\theta}G\left(\frac{\theta}{2\beta_{2}}, -\left(\frac{\alpha_{0}}{z}\right)^{2}K^{2\beta_{2}}\right)-G\left(\frac{\theta}{2\beta_{2}}, -\left(\frac{\alpha_{0}}{z}\right)^{2}\right)\nonumber\\ \qquad=\left(K^{\theta}-1\right)+\sum\limits_{j=1}^{\infty}\frac{\theta}{(\theta+2\, j\, \beta_{2})}\left(K^{\theta}\left(-\left(\frac{\alpha_{0}}{z}\right)^{2}K^{2\beta_{2}}\right)^{j}-\left(-\left(\frac{\alpha_{0}}{z}\right)^{2}\right)^{j}\right)\nonumber\\ \qquad=\sum\limits_{j=0}^{\infty}\frac{\theta}{(\theta+2\, j\, \beta_{2})}(-1)^{j}\alpha_{0}^{2j}\left(K^{2j\beta_{2}+\theta}-1\right)\frac{1}{z^{2j}}$

(56)

---

3. The mean reverting CMLJD model

As mentioned in the introduction, a simple jump diffusion is not appropriate to represent commodities as their prices revert to long run equilibrium prices. To insert this feature in assets dynamics, the following mean reversion mechanism is studied

$X_{t} = \varphi(t)+Y_{t}\, ,$

(57)

where $\varphi(t)$ is a function of time defined by

$\varphi(t) = b\;(1-e^{-at})\, .$

(58)

$b$ is the constant mean reversion level whereas $a$ is the speed of mean reversion. On the other hand, $Y_{t}$ is non Gaussian Ornstein-Uhlenbeck process with $Y_{0}=X_{0}$, driven by the next SDE:

$dY_{t} = -aY_{t}dt+dZ_{t}.$

(59)

where $dZ_{t}$ is a Lévy process, sum of a Brownian component and of a jump process:

$dZ_{t}: =\sigma dW_{t}+dL_{t}\, .$

(60)

As previously, $L_{t}=\sum_{j=1}^{N_{t}}J_{j}$ where $N_{t}$ is a Poisson process and $J$ is a continuous mixture of Laplace's law. $L_{t}$ is the limit in the weak sense of the sum of processes $L_{t}^{k}$ when $\Delta_{k}$ tends to zero. In this setting, Applying the Lévy Ito formula to $e^{at}Y_{t}$ leads to the following expression for $Y_{t}$,

$Y_{t}=Y_{s}e^{-a(t-s)}+\int_{s}^{t}e^{-a(t-u)}dZ_{u}\, .$

(61)

The statistical distribution of $Y_{s}$ is unknown but may be inferred from its characteristic function in numerical applications. The asset value at time $t$, conditionally to the available information at time $s$ is given by:

$A_{t} = A_{s}\exp\left(\varphi(t)-\varphi(s)+Y_{s}e^{-a(t-s)}+\int_{s}^{t}e^{-a(t-u)}dZ_{u}\right)\, .$

(62)

Given that $Y_{0}=X_{0}$, the characteristic function of the asset return is:

$\begin{eqnarray} M_{X_{t}}(z) & = & \mathbb{E}\left(e^{izX_{t}}|\mathcal{F}_{0}\right)\nonumber\\ & = & e^{iz\left(\varphi(t)+X_{0}e^{-at}\right)}\mathbb{E}\left(e^{\int_{0}^{t}ize^{-a(t-u)}dZ_{u}}|\mathcal{F}_{0}\right) \end{eqnarray}$

(63)

The expectation is valued by the following result, proposed by Eberlein and Raible (1999):

Proposition 3.1. Let $Z_{t}$ be a Lévy process having a cumulant transform defined as follows

$\tilde{\psi}(u) = \log\mathbb{E}(\exp(uZ_{1}))\, ,$

(64)

and let f : $\mathbb{R}_{+}\to\mathbb{C}$ be a complex valued left continuous function such that $|Re(f)|\leq M$ then

$\mathbb{E}\left(\exp\left(\int_{0}^{t}f(u)\, dZ_{\theta}\right)|\mathcal{F}_{0}\right) = \exp\left(\int_{0}^{t}\, \tilde{\psi}(f(u))\, du\right)\, .$

(65)

In particular, if $Z_{t}$ is a mixed Laplace process, its cumulant transform is equal to:

$\tilde{\psi}(u) = \frac{1}{2}\sigma^{2}u^{2}+\int_{1}^{K}\lambda_{k}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}-u^{2}}dk-\lambda\, ,$

(66)

and

$f(u)=ize^{-a(t-u)}\, .$

(67)

Proposition 3.2. The characteristic function of $X_{t}$ is equal to

$\begin{eqnarray} M_{X_{t}}(z) & := & \exp\left(\psi(t, z)\right)\\ & = & \exp\left(iz\left(\varphi(t)+X_{0}e^{-at}\right)+\int_{0}^{t}\, \tilde{\psi}(ize^{-a(t-u)})\, du\right)\, , \end{eqnarray}$

(68)

where the integral $\int_{0}^{t}\, \tilde{\psi}(ize^{-a(t-u)})\, du$ is given by the next expression:

$\int_{0}^{t}\, \tilde{\psi}(ize^{-a(t-u)})\, du=-\frac{1}{4a}\sigma^{2}z^{2}\left(1-e^{-2at}\right)+\label{eq:intGammaPhiAnalytic}\\ \qquad\frac{\lambda_{0}}{2a}\frac{2\beta_{2\, }K^{\beta_{1}+1}}{\left(\beta_{1}+1\right)^{2}}\left(G\left(\frac{\beta_{1}+1}{2\beta_{2}}, -\frac{\alpha_{0}^{2}}{z^{2}}K^{2\beta_{2}}e^{2at}\right)-G\left(\frac{\beta_{1}+1}{2\beta_{2}}, -\frac{\alpha_{0}^{2}}{z^{2}}K^{2\beta_{2}}\right)\right)\nonumber \\ \qquad-\frac{\lambda_{0}}{2a}\frac{2\beta_{2}}{\left(\beta_{1}+1\right)^{2}}\left(G\left(\frac{\beta_{1}+1}{2\beta_{2}}, -\frac{\alpha_{0}^{2}}{z^{2}}e^{2at}\right)-G\left(\frac{\beta_{1}+1}{2\beta_{2}}, -\frac{\alpha_{0}^{2}}{z^{2}}\right)\right)\nonumber \\ \qquad+\frac{\lambda_{0}}{2a}\frac{K^{\beta_{1}+1}}{\left(\beta_{1}+1\right)}\ln\left(\frac{\alpha_{0}^{2}K^{2\beta_{2}}+z^{2}e^{-2at}}{\alpha_{0}^{2}K^{2\beta_{2}}+z^{2}}\right)-\frac{\lambda_{0}}{2a}\frac{1}{\left(\beta_{1}+1\right)}\ln\left(\frac{\alpha_{0}^{2}+z^{2}e^{-2at}}{\alpha_{0}^{2}+z^{2}}\right)\, ,$

(69)

and where $G(b, x)$ is again the hypergeometric function:

$\begin{eqnarray} G(b, x) & = & _{2}F_{1}(1, b, b+1, x)\nonumber \\ & = & b\int_{0}^{1}\frac{u^{b-1}}{1-ux}du\, . \end{eqnarray}$

(70)

Proof. A direct calculation leads to the following development:

$\begin{eqnarray} \int_{0}^{t}\, \tilde{\psi}(ize^{-a(t-u)})\, du & = & -\frac{1}{4a}\sigma^{2}z^{2}\left(1-e^{-2at}\right)+\label{eq:IntphiTilde-1}\\ & & \int_{1}^{K}\lambda_{k}\int_{0}^{t}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z^{2}e^{-2a(t-u)}}du\, dk-\lambda t\, , \end{eqnarray}$

(71)

and the integral in the second term is equal to

$\begin{eqnarray} \int_{0}^{t}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z^{2}e^{-2a(t-u)}}du & = & \left[u-\frac{1}{2a}\ln\left(1+\left(\frac{z^{2}}{\alpha_{k}^{2}}e^{-2at}\right)e^{2au}\right)\right]_{u=0}^{u=t}\nonumber\\ & = & t+\frac{1}{2a}\ln\left(\frac{\alpha_{k}^{2}+z^{2}e^{-2at}}{\alpha_{k}^{2}+z^{2}}\right). \end{eqnarray}$

(72)

Then the integral in equation (71) becomes:

$\begin{eqnarray} \int_{1}^{K}\lambda_{k}\int_{0}^{t}\frac{\alpha_{k}^{2}}{\alpha_{k}^{2}+z^{2}e^{-2a(t-u)}}du\, dk & = & \lambda t+\frac{1}{2a}\int_{1}^{K}\lambda_{k}\ln\left(\frac{\alpha_{k}^{2}}{z^{2}}+e^{-2at}\right)\, dk\nonumber \\ & & -\frac{1}{2a}\int_{1}^{K}\lambda_{k}\ln\left(\frac{\alpha_{k}^{2}}{z^{2}}+1\right)\, dk\, . \end{eqnarray}$

(73)

For any constant $\kappa$, several changes of variables similar to these done in proposition 2.3, lead to the following expression for the integral:

$\int_{1}^{K}\lambda_{k}\ln\left(\frac{\alpha_{k}^{2}}{z^{2}}+\kappa\right)\, dk=\lambda_{0}\, \int_{1}^{K}k^{\beta_{1}}\ln\left(\frac{\alpha_{0}^{2}}{z^{2}}\, k^{2\beta_{2}}+\kappa\right)\, dk=\lambda_{0}\frac{1}{\left(\beta_{1}+1\right)^{2}}\times\nonumber \\ \left[k^{\beta_{1}+1}\left(2\beta_{2\, }G\left(\frac{\beta_{1}+1}{2\beta_{2}}, -\frac{\alpha_{0}^{2}k^{2\beta_{2}}}{z^{2}\kappa}\right)+\left(\beta_{1}+1\right)\ln\left(\frac{\alpha_{0}^{2}}{z^{2}}k^{2\beta_{2}}+\kappa\right)-2\beta_{2}\right)\right]_{k=1}^{k=K}.$

(74)

Combining equations (71) (73) and (74) allows us to conclude the proof.

Notice that the dynamics of $X_{t}$ can be described by the following SDE

$dX_{t} = a\left(b-X_{t}\right)dt+\sigma dW_{t}+dL_{t},$

(75)

We deduce from this last relation, its infinitesimal generator is equal to

$\begin{eqnarray} \left(\mathcal{L}u\right)(x) & = & a\left(b-x\right)\frac{\partial}{\partial x}u(x)+\frac{\sigma^{2}}{2}\frac{\partial}{\partial x^{2}}u(x)+\label{eq:InfinitesimalGen-1}\\ & & \lambda\int_{-\infty}^{+\infty}\left(u(x+y)-u(x)\right)\mu^{J}(y)dy\, , \end{eqnarray}$

(76)

for any function $u(x)$ that is twice continuously derivable and where $\mu^{J}(.)$ is given by (28). This generator is used for pricing purposes in appendix A.

The moments of $X_{t}$ are then obtained by differentiating the characteristic function. Again, the skewness is null and the kurtosis is always above 3.

Proposition 3.3. The mean, variance, skewness and kurtosis of $X_{t}$ are respectively given by:

$\mathbb{E}(X_{t}) = b\;(1-e^{-at})+X_{0}e^{-at}\, ,$

(77)

$\mathbb{V}(X_{t}) = \frac{1}{2a}(1-e^{-2at})\left(\sigma^{2}+2\frac{\lambda_{0}}{\alpha_{0}^{2}}\frac{1}{\beta_{1}-2\beta_{2}+1}\left(K^{\beta_{1}-2\beta_{2}+1}-1\right)\right)\, ,$

(78)

$\mathbb{S}(X_{t})=0\, ,$

(79)

$\mathbb{K}(X_{t}) =3+\frac{24\frac{1}{4a}\left(1-e^{-4at}\right)\left(\frac{\lambda_{0}}{\alpha_{0}^{4}}\frac{1}{\beta_{1}-4\beta_{2}+1}\left(K^{\beta_{1}-4\beta_{2}+1}-1\right)\right)}{\left(\frac{1}{2a}(1-e^{-2at})\left(\sigma^{2}+2\frac{\lambda_{0}}{\alpha_{0}^{2}}\frac{1}{\beta_{1}-2\beta_{2}+1}\left(K^{\beta_{1}-2\beta_{2}+1}-1\right)\right)\right)^{2}}\, .$

(80)

Proof. The moments of $X_{t}$ are obtained by differentiating the characteristic function with respect to $z$,

$\mathbb{E}(X_{t}^{k}) = \left.\frac{\partial^{k}}{\partial z^{k}}M_{X_{t}}(-iz)\right|_{z=0}\, .$

(81)

In particular, if $g(t)$ denotes the following function,

$g(t) = \sigma^{2}\int_{0}^{t}\, e^{-2a(t-u)}\, du+2\int_{0}^{t}\int_{1}^{K}\frac{\lambda_{k}}{\alpha_{k}^{2}}e^{-2a(t-u)}dk\, du\, ,$

(82)

the non centered moments of $X_{t}$ are

$\mathbb{E}(X_{t}) = b\;(1-e^{-at})+X_{0}e^{-at}\, ,$

(83)

$\mathbb{E}(X_{t}^{2}) = \left(b\;(1-e^{-at})+X_{0}e^{-at}\right)^{2}+g(t)\, ,$

(84)

$\mathbb{E}(X_{t}^{3}) = \left(b\;(1-e^{-at})+X_{0}e^{-at}\right)^{3}+3\left(b\;(1-e^{-at})+X_{0}e^{-at}\right)g(t)\, ,$

(85)

$\mathbb{E}(X_{t}^{4}) = \left(b\;(1-e^{-at})+X_{0}e^{-at}\right)^{4}+6\left(b\;(1-e^{-at})+X_{0}e^{-at}\right)^{2}g(t)\nonumber\\ \quad \quad +3g(t)^{2}+24\int_{0}^{t}\int_{1}^{K}\frac{\lambda_{k}}{\alpha_{k}^{4}}e^{-4a(t-u)}dkdu\, .$

(86)

If the hypergeometric function is developed as an infinite serie, the following result speeds up the numerical calculation of the moment generating function:

Corollary 3.4. The integral $\int_{0}^{t}\, \tilde{\psi}(ize^{-a(t-u)})\, du$ is equal to the following sum:

$\int_{0}^{t}\, \tilde{\psi}(ize^{-a(t-u)})\, du=-\frac{1}{4a}\sigma^{2}z^{2}\left(1-e^{-2at}\right)+\label{eq:intGammaPhiAnalytic-1-1-1}\\ \qquad\frac{\lambda_{0}}{2a}\frac{2\beta_{2}}{\left(\beta_{1}+1\right)^{2}}\left(\sum\limits_{j=1}^{\infty}\frac{1}{1+\frac{2\beta_{2}}{\beta_{1}+1}j}\left(K^{2j\beta_{2}+\beta_{1}+1}-1\right)\left(-\frac{\alpha_{0}^{2}}{z^{2}}\right)^{j}\left(e^{2at}-1\right)^{j}\right)\nonumber \\ \qquad+\frac{\lambda_{0}}{2a}\frac{K^{\beta_{1}+1}}{\left(\beta_{1}+1\right)}\ln\left(\frac{\alpha_{0}^{2}K^{2\beta_{2}}+z^{2}e^{-2at}}{\alpha_{0}^{2}K^{2\beta_{2}}+z^{2}}\right)-\frac{\lambda_{0}}{2a}\frac{1}{\left(\beta_{1}+1\right)}\ln\left(\frac{\alpha_{0}^{2}+z^{2}e^{-2at}}{\alpha_{0}^{2}+z^{2}}\right)\, .$

(87)

---

4. Calculation of the probability density function

The calculation of characteristic exponents can be done numerically by discretizing the integral form of $G(b, x)$ or with the development in equation (53), truncated to a finite number of terms. Both approaches may be used in numerical applications to retrieve the density function of CMLJD processes with and without mean reversion, by a discrete Fourier transform. The next proposition presents this methodology:

Proposition 4.1. Let $N$ be the number of steps used in the Discrete Fourier Transform (DFT) and $\Delta_{x}=\frac{2x_{max}}{N-1}$, be the step of discretization. Let us denote $\delta_{j}=\frac{1}{2}1_{\{j=1\}}+1_{\{j\neq1\}}$, $\Delta_{z}=\frac{2\pi}{N\, \Delta_{x}}$ and $z_{j}=\left(j-1\right)\Delta_{z}$. If $X_{t}$ is CMLJD, the values of $f_{X_{t}}(.)$ at points $x_{k}=-\frac{N}{2}\Delta_{x}+(k-1)\Delta_{x}$ are approached by

$f_{X_{t}}(x_{k})=\frac{2}{N\, \Delta_{x}}\sum\limits_{j=1}^{N}\delta_{j}\left(e^{t\psi(z_{j})}(-1)^{j-1}\right)e^{-i\frac{2\pi}{N}\left(j-1\right)(k-1)}\, ,$

(88)

where $\psi(z)$ is defined by equation (35). If $X_{t}$ is CMLJD with mean reversion, the function $f_{X_{t}}(.)$ is approached by

$f_{X_{t}}(x_{k}) = \frac{2}{N\, \Delta_{x}}\sum\limits_{j=1}^{N}\delta_{j}\left(e^{\psi(t, z_{j})}(-1)^{j-1}\right)e^{-i\frac{2\pi}{N}\left(j-1\right)(k-1)}\, ,$

(89)

where $\psi(t, z)$ is defined by equation (35). This last relation can be computed by any DFT procedure.

Proof. By definition, the characteristic function is the inverse Fourier transform of the density

$M_{X_{t}}(z) =\int_{-\infty}^{+\infty}f_{X_{t}}(x)e^{izx}dx\nonumber\\ \quad :=2\pi\mathcal{F}^{-1}[f_{X_{t}}(x)](z)\, .$

(90)

If $X_{t}$ is CMLJD, the density is retrieved by calculating the Fourier transform of $M_{X_{t}}(z)$ as follows

$\begin{eqnarray} f_{X_{t}}(x) & = & \frac{1}{2\pi}\mathcal{F}[e^{t\psi(z)}](x)\nonumber\\ & = & \frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{t\psi(z)}e^{-ixz}dz\nonumber\\ & = & \frac{1}{\pi}\int_{0}^{+\infty}e^{t\psi(z)}e^{-ixz}dz\, . \end{eqnarray}$

(91)

The last equality comes from the fact that $\psi(z)$ and $\psi(-z)$ are complex conjugate. Equation (88) is deduced by approaching this last integral with the trapezoid rule $\int_{a}^{b}h(x)dx=\frac{h(a)+h(b)}{2}\Delta_{x}+\sum_{k=1}^{N-1}h(a+k\, \Delta_{x})\Delta_{x}$. Equation (89) is proven in the same way.

The CMLJD process without and with mean reversion are respectively identified by 7 parameters ($\mu$, $\sigma$, $\alpha_{0}$, $\lambda_{0}$, $\beta_{1}$, $\beta_{2}$, $K$) and 8 parameters ($a$, $b$, $\sigma$, $\alpha_{0}$, $\lambda_{0}$, $\beta_{1}$, $\beta_{2}$, $K$). Both processes cannot be fitted to time series by the method of moments matching, without setting a priori some parameters. An alternative consists to calibrate the process by log-likelihood maximization. We adopt this approach in numerical illustrations to fit CMLJD processes. The matlab code implementing the equation (88) is provided in appendix B and may be used to evaluate the expression (89). The characteristic exponents are computed by direct integration, with equation (2.3) and (3.2). The matlab code implementing these operations is also reported in appendix B.

---

5. Options pricing

Firstly, we consider that log-returns are ruled by a CMLJD process without mean reversion. The pricing of financial securities is done under a risk neutral measure. Under this measure, the discounted price process is a martingale and the expected return is equal to the risk free rate, $r$ to avoid any arbitrage. $X_{t}$ is then defined by parameters $(r, \sigma, \alpha_{0}, \beta_{1}, \lambda_{0}, \beta_{2}, K)$. The most common methods used for pricing derivatives are based on Fourier and Inverse Fourier transforms. We denote them respectively by:

$\mathcal{F}[f](\omega) = \int_{-\infty}^{+\infty}f(x)e^{-i\omega x}dx\;, \;\mathcal{F}^{-1}[f](x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}f(\omega)e^{i\omega x}d\omega.$

(92)

The Fourier transform maps spatial derivatives $\frac{\partial}{\partial x}$ into multiplications in the frequency domain. As shown in the next result, this feature allows us to price any European derivatives.

Proposition 5.1. In the CMLJD model, the price at time $t$ and when $X_{t}=x$, of an European derivatives delivering a payoff $V(T, X_{T}$) at maturity $T,$ is given by

$\left[V(t, x)\right]=\mathcal{F}^{-1}\left[\mathcal{F}\left[V(T, x)\right](\omega)e^{\left(\psi(\omega)-r\right)(T-t)}\right](x).$

(93)

Proposition 5.2. In the CMLJD model with mean reversion, the price at time $t$ and when $Y_{t}=y$, of an European derivatives delivering a payoff $V(T, Y_{T}$) at maturity $T,$ is given by

$\left[V(t, y)\right]=\mathcal{F}^{-1}\left[\mathcal{F}\left[V(T, y)\right](e^{a(T-t)}\omega)e^{\int_{0}^{T-t}\tilde{\psi}(ie^{a(T-t-u)}\omega)du-(r-a)(T-t)}\right](y).$

(94)

where $\int_{0}^{T-t}\tilde{\psi}(ie^{a(T-t-u)}\omega)du$ is provided by equation (69).

The proofs are standard in the literature and are reproduced in appendix A for information. Jackson et al. (2008) an Jaimungal and Surkov (2011) have used iteratively a similar procedure to price American or barrier options for other Lévy processes, with or without mean reversion. In practice, integrals in equations (93) or (5.2) are discretized, and the price is obtained by Discrete Fourier Transforms as stated in next propositions.

Proposition 5.3. Let $N$ and $\Delta_{x}=\frac{2x_{max}}{N-1}$ be respectively the number of steps used in the Discrete Fourier Transforms (DFT) and the step of discretization. Let us define $\Delta_{\omega}=\frac{2\pi}{N\, \Delta_{x}}$ and $\delta_{k}=\frac{1}{2}1_{\{k=1\}}+1_{\{k\neq1\}}$. $V(t, x_{j})$ in equations (93), at points $x_{j}=-\frac{N}{2}\Delta_{x}+(j-1)\Delta_{x}$, for $j=1...N$, is approached by the following DFTs sum:

$V(t, x_{j}) \approx \frac{2}{N}\sum\limits_{k=1}^{N}\left(\delta_{k}e^{\left(\psi(\omega_{k})-r\right)(T-t)}\sum\limits_{m=1}^{N}V(T, x_{m})e^{-i\frac{2\pi}{N}(k-1)(m-1)}\right)e^{i\frac{2\pi}{N}(k-1)(j-1)}\, ,$

(95)

where $\omega_{k}=(k-1)\Delta_{\omega}$.

Proof. The Fourier transform is approached by the following sum

$\begin{eqnarray} \mathcal{F}[V(T, x)](\omega_{k}) & = & \int_{-\infty}^{+\infty}V(T, x)e^{-i\omega_{k}x}dx\nonumber\\ & \approx & \Delta_{x}e^{-i(k-1)\Delta_{\omega}x_{min}}\sum\limits_{m=1}^{N}V(T, x_{m})e^{-i(k-1)(m-1)\Delta_{\omega}\Delta_{x}}\, . \end{eqnarray}$

(96)

As $\Delta_{\omega}\Delta_{x}=\frac{2\pi}{N}$ and $x_{min}=-\frac{N}{2}\Delta_{x}$, this last expression becomes

$\mathcal{F}[V(T, x)](\omega_{k}) \approx \left(\Delta_{x}e^{i(k-1)\pi}\right)\sum\limits_{m=1}^{N}V(T, x_{m})e^{-i\frac{2\pi}{N}(k-1)(m-1)}\, .$

(97)

Let us denote $g(\omega)=\mathcal{F}\left[V(T, x)\right](\omega)e^{\left(\psi(\omega)-r\right)(T-t)}$ then

$\begin{eqnarray} V(t, x) & = & \mathcal{F}^{-1}[g(\omega)](x)\nonumber\\ & = & \frac{1}{\pi}\int_{0}^{+\infty}g(\omega)e^{i\omega x}d\omega\, . \end{eqnarray}$

(98)

The function being known at point $\omega_{k}=(k-1)\Delta_{\omega}$, approaching this last integral with the trapezoid rule, leads to:

$\begin{eqnarray} V(t, x_{j}) & \approx & \frac{1}{\pi}\sum\limits_{k=1}^{N}\delta_{k}g(\omega_{k})e^{i\omega_{k}x_{j}}\Delta_{\omega}\nonumber\\ & \approx & \frac{\Delta_{\omega}}{\pi}\sum\limits_{k=1}^{N}\delta_{k}\left(g(\omega_{k})e^{-i(k-1)\pi}\right)e^{i\frac{2\pi}{N}(k-1)(j-1)}\, , \end{eqnarray}$

(99)

which is equivalent to (95).

Proposition 5.4. Let $N$ and $\Delta_{y}=\frac{2y_{max}}{N-1}$ be respectively the number of steps used in the Discrete Fourier Transforms (DFT) and the step of discretization. Let us define $\Delta_{\omega}=\frac{2\pi}{N\, \Delta_{y}}$ and $\delta_{k}=\frac{1}{2}1_{\{k=1\}}+1_{\{k\neq1\}}$. $V(t, y_{j})$ in equations (5.2), at points $y_{j}=-\frac{N}{2}\Delta_{x}+(j-1)\Delta_{y}$, for $j=1...N$, is approached by the following DFTs sum:

$\begin{eqnarray} V(T, y_{j}) & = & \frac{2}{N}\sum\limits_{k=1}^{N}\delta_{k}e^{\int_{0}^{T-t}\tilde{\psi}(ie^{a(T-t-u)}\omega_{k})du-r(T-t)}\left(\sum\limits_{m=1}^{N}\left[V(T, e^{-a(T-t)}y_{m})\right]e^{-i\frac{2\pi}{N}(k-1)(m-1)}\right)\nonumber\\ & & \times e^{i\frac{2\pi}{N}(k-1)(j-1)}\, , \end{eqnarray}$

(100)

where $\omega_{k}=(k-1)\Delta_{\omega}$.

Proof. By a change of variable $y'=e^{a(T-t)}y$ and if $\Delta_{y'}=\Delta_{y}$, the Fourier transform of the terminal is approached by a DFT as follows

$\begin{eqnarray} \mathcal{F}[V(T, y)](e^{a(T-t)}\omega_{k}) & = & \int_{-\infty}^{+\infty}V(T, y)e^{-ie^{a(T-t)}\omega_{k}y}dy\nonumber\\ & = & \int_{-\infty}^{+\infty}\left(V(T, e^{-a(T-t)}y')e^{-a(T-t)}\right)e^{-i\omega_{k}y'}dy'\nonumber\\ & \approx & \Delta_{y'}e^{-i(k-1)\Delta_{\omega}y'_{min}}\sum\limits_{m=1}^{N}\left[V(T, e^{-a(T-t)}y'_{m})e^{-a(T-t)}\right]e^{-i(k-1)(m-1)\Delta_{\omega}\Delta_{y'}}\, . \end{eqnarray}$

(101)

As $\Delta_{\omega}\Delta_{y'}=\frac{2\pi}{N}$ and $y'_{min}=-\frac{N}{2}\Delta_{y'}$, this last expression becomes

$\mathcal{F}[V(T, y)](e^{a(T-t)}\omega_{k}) \approx \Delta_{y'}e^{i(k-1)\pi}\sum\limits_{m=1}^{N}\left[V(T, e^{-a(T-t)}y'_{m})e^{-a(T-t)}\right]e^{-i\frac{2\pi}{N}(k-1)(m-1)}\, .$

(102)

Let us denote $g(\omega)=\mathcal{F}\left[V(T, y)\right](e^{a(T-t)}\omega)e^{\int_{0}^{T-t}\tilde{\psi}(ie^{a(T-t-u)}\omega)du-(r-a)(T-t)}$ then

$\begin{eqnarray} V(t, y) & = & \mathcal{F}^{-1}[g(\omega)](y)\nonumber\\ & = & \frac{1}{\pi}\int_{0}^{+\infty}g(\omega)e^{i\omega y}d\omega\, . \end{eqnarray}$

(103)

The function being known at point $\omega_{k}=(k-1)\Delta_{\omega}$, if we approach this last integral with the trapezoidal rule, we infer that:

$\begin{eqnarray} V(t, y_{j}) & \approx & \frac{1}{\pi}\sum\limits_{k=1}^{N}g(\omega_{k})e^{i\omega_{k}y_{j}}\Delta_{\omega}\nonumber\\ & \approx & \frac{\Delta_{\omega}}{\pi}\sum\limits_{k=1}^{N}\left(g(\omega_{k})e^{-i(k-1)\pi}\right)e^{i\frac{2\pi}{N}(k-1)(j-1)}\, , \end{eqnarray}$

(104)

and we can conclude.

---

6. Numerical applications

Firstly, we fit the CMLJD to daily log-returns of four stocks indices, over the period June 2006 to June 2016 (2600 observations). The chosen indices are the Morgan Stanley World stocks indice, the FTSE 100, the S & P 500 and the CAC 40. Given that $\mathbb{E}\left(dX_{t}\right)=\mu dt$, we set the drift to the average of log-returns. The calibration of other parameters is done by log-likelihood maximization. The density is retrieved numerically by inverting the characteristic function of the CMLJD process. The number of steps for the DFT is set to $N=2^{18}$, the minimum and maximum daily returns are respectively equal to $-\frac{N}{2}\Delta_{x}=-0.30$ and $\frac{N}{2}\Delta_{x}=0.30$. The characteristic exponent is computed numerically with equation (35). The CMLJD is also compared to the DEJD model of Kou (2002) that postulates the following dynamics for the log-return:

$dX_{t}^{DEJD}=\mu dt+\sigma dW_{t}+dL_{t}^{DEJ},$

(105)

where $L_{t}^{DEJ}:=\sum_{j=1}^{N_{t}^{DEJ}}J_{j}^{DEJ}\, .$ In this last expression, $N_{t}^{DEJ}$ is a Poisson process with a constant intensity $\lambda^{DEJ}$ and $J_{j}^{DEJ}$ are distributed according to a double exponential law with a density:

$\mu_{DEJ}(x) = p\lambda^{+}e^{-\lambda^{+}x}1_{\{x\geq0\}}-(1-p)\lambda^{-}e^{-\lambda^{-}y}1_{\{x<0\}}$

(106)

where $p$ and $\lambda^{+}$ are positive constants and $\lambda^{-}$ is a negative constant. They represent the probability of observing respectively upward and downward exponential jumps. The expectation of $J$ is set to zero by setting $\lambda^{-}=-\frac{p}{1-p}\lambda^{+}$. This ensures the identifiability of the model and that $L_{t}^{DEJ}$ is a martingale. The characteristic exponent of this process is equal to

$\psi^{DEJ}(z) := i\mu z-\frac{1}{2}\sigma^{2}z^{2}+p\frac{\lambda^{+}}{\lambda^{+}-z\, i}-(1-p)\frac{\lambda^{-}}{z\, i-\lambda^{-}}\, .$

(107)

As the CMLJD, the DEJD does not admit analytical expression for its probability density function. The same DFT algorithm is used to compute it. Fitted parameters, standard errors^* and log-likelihoods are reported in table 1. The CMLJD consistently outperform the Brownian motion and the DEJD model as underlined by the comparison of log-likelihoods (lines "LogLik.", "DEJD LogLik." and "B.M. LogLik." in exhibit 1). Parameters are well behaved and consistent in so much that they exhibit stability and some form of erratic variations among the different index. The Brownian volatility $\sigma$ is between 5.81% and 12.22%: the lowest Brownian volatility being obtained for the most diversified stocks indice (MS world). The $\lambda_{0}$'s that measure the average number of jumps per year are in the range $[22.80\, , \, 27.72]$. At least for the MS world, FTSE 100 and S & P 500, $K$ and $\alpha_{0}$ are comparable.

^*The standard error of estimation for a parameter $\theta\in\left(\mu, \sigma, \lambda_{0}, \beta_{1}, \beta_{2}, \alpha_{0}, K\right)$ is computed numerically as the square root of the asymptotic Fisher information:

$Std.\, err.\, (\theta) = \sqrt{-\left(\frac{\partial^{2}\ln L(\theta)}{\partial\theta^{2}}\right)^{-1}}.$

The table 2 compares the empirical moments of daily returns with these of CMLJD processes. These figures clearly confirm that CMLJD's exhibit a high kurtosis, comparable to the one displayed by observations. By construction, the skewness is null as the distribution of jump is symmetric.

To check that the CMLJD process duplicates smiles of implied volatilities similar to these observed in financial markets, European call prices on the S & P 500 are computed with fitted parameters of table 1. Next, implied volatilities are retrieved by inverting the Black & Scholes formula. The volatility surface obtained by this method is shown in the left graph of exhibit 1. The shape of this surface is coherent and realistic. For short term maturities, the smile of volatilities is well visible and a minimum is attained for "At The Money" option (moneyness =100). The curvature of the smile is inversely proportional to time to maturity. The implied volatility of 1 year and 1 month ATM options are respectively 21.09% and 23.72%.

Figure 2 illustrates the sensitivity of S & P 500 implied volatilities to modification of parameters. The time to maturity is 1 month. Keeping all other parameters equal to these of table 1, increasing $\lambda_{0}$ or $K$ shift up the volatility surface, without important modification of the curvature. Whereas raising $\beta_{2}$ or $\alpha_{0}$ moves down the curve and slightly flattens the smile of volatilities.

Tests of the mean reverting CMLJD are carried out on four times series of commodities: Copper (LME), Soy Beans (yellow soybeans, Chicago), Crude Oil WTI spot and Wheat (USDA 2 soft red winter wheat, Chicago) from December 2011 to November 2015 (1000 observations). It is well known that commodities exhibit cointegration in prices. Like equity prices they are also exposed to sudden jumps of price. However, unlike equities, commodities tend to revert to long run equilibrium prices. Parameters are fitted by maximization of the log-likelihood. However the calculation of this log-likelihood requires a pre-treatment of data to reduce the computation time. Firstly, we calculate the log-return process, $X_{j}=\frac{P_{j}}{P_{0}}$ where $P_{j}$ is the price of commodity on the $j^{th}$ day. The process $Y_{j}$ is next obtained by subtracting from $X_{j}$ the function $\varphi(t_{j})$ as defined by equation (58),

$Y_{j}=X_{j}-\varphi(t_{j})\, .$

(108)

If $\Delta_{t}$ represents one day of trading, according to equation (61), the process $V_{j}$ is distributed as follows

$V_{j}=Y_{j}-Y_{j-1}e^{-a\Delta_{t}}\sim\int_{0}^{\Delta_{t}}e^{-a(t-u)}dZ_{u}\, .$

(109)

The distribution of $V_{j}$ is finally retrieved by inverting numerically the characteristic function of $\int_{0}^{\Delta_{t}}e^{-a(t-u)}dZ_{u}$, that is reported in proposition 3.2. The log-likelihood is computed with this distribution. Results of the calibration procedure are shown in table 3. Soy Beans and Copper prices display respectively the lowest and the highest speed of mean reversion. Furthermore, log-likelihoods of the mean reverting CMLJD are significantly higher than these obtained with and without a mean reverting Brownian motions (lines "M.R. B.M. LogLik." and"B.M. LogLik." in the exhibit 3).

As for stocks indices, fitted parameters exhibit stability among the different times series. The frequency of jumps $\lambda_{0}$ is in the range $[58.51\, , \, 65.34]$. $\beta_{1}$ is small (around 4%). $\beta_{2}$ for Crude Oil and Wheat are similar and close to 30% whereas $\beta_{2}$ of Soy beans and Copper are higher. $K$ varies a lot: from 3.58 for Crude Oil to 25 for Copper. The parameter $\alpha_{0}$ is in a corridor {[}73.50, 90.64{]}. Table 4 compares the moments of daily returns with these of models. Again, the CMLJD demonstrates its capacity to fit the variance and leptokurticity of time series.

According to the proposition 3.3, the expectation of $X_{t}$ converges to $b$ when $t$ tends to $\infty$. Then, the estimated $b$ should be close to the average of $X_{t}$ for the considered period. The first column of table 4 confirms this intuition. On the other hand, the parameter $b$ is only involved in the mean of the process and does not influence the variance and kurtosis. Then, if we calibrate the model with data observed on another time window, $b$ will be different if and only if we observe a significant modification of the average of $X_{t}$.

To assess the impact of the mean reversion on implied volatilities, European call options on Wheat are priced with fitted parameters. The surface of implied volatilities obtained by inverting the Black & Scholes formula is displayed in the right graph of exhibit 1. Contrary to equities, we don't observe a smile of volatilities, with a local minimum for ATM options. Instead, volatilities are convex and inversely proportional to the moneyness. Given that a similar trend is not observed for the S & P500 volatility surface, this behavior may be associated to the mean reversion of returns.

Figure 3 illustrates the sensitivity of Wheat implied volatilities to parameters $a$ and $b$. Two maturities are considered: 3 and 6 months. Keeping all other parameters equal to these of table 3, increasing the speed of mean reversion $a$ raises the steepness of the volatility surface, without any strong modification of the curvature. On the contrary, raising the level of mean reversion $b$, shifts up the curve in an asymmetrical manner.

---

7. Conclusion

This article proposes two parsimonious models for stocks or commodities, driven by a diffusion and a mixture of Laplace jumps processes. These dynamics aim to replicate time series for which increments have a distribution with a high kurtosis. Despite the somewhat lengthy expressions of characteristic functions, the CMLJD remains analytically tractable and its density can be retrieved numerically by a discrete Fourier transform. Its ability to duplicate leptokurtic processes makes it eligible for option pricing or risk management.

As underlined by the empirical study, the CMLJD processes fit stocks and commodities better than a Brownian motion or a DEJD. Furthermore, the parameters estimated by log-likelihood maximization show stability and consistency over a variety of assets. On the other side, CMLJD processes generate realistic surfaces of implied volatilities. In the CMLJD model with mean reversion, we observe a strong asymmetry in this surface, due to the reversion of prices.

There exist many possibilities for further researches. Firstly, the main criticism about the CMLJD is that it does not capture the asymmetry of returns. It is possible to remedy to this problem by considering double exponential jumps. However, this increases the number of parameters and requires additional parameterization. Secondly, we only consider power functions for $\alpha_{k}$ and $\lambda_{k}$. It is probably possible to find another type of functions for which we can obtain a closed form expression for the limit of the characteristic function of $X_{t}$. Finally, it would be interesting to develop a multivariate extension of this model.

---

Appendix A

The Fourier transform maps spatial derivatives $\frac{\partial}{\partial x}$ into multiplications in the frequency domain. For any differentiable function $f$, we have then:

$\mathcal{F}\left[\frac{\partial^{n}}{\partial x^{n}}f\right](t, \omega) = \left(i\omega\right)^{n}\mathcal{F}[f](t, \omega)$

(110)

and

$\mathcal{F}\left[x\frac{\partial}{\partial x}f\right](t, \omega) = -\mathcal{F}[f](t, \omega)-\frac{\partial}{\partial\omega}\mathcal{F}[f](t, \omega)$

(111)

If $X_{t}$ is a CMLJD process, the price of any derivatives is arbitrage free if and only it is solution of the Feynman-Kac equation:

$\frac{\partial}{\partial t}V+\mathcal{L}V = rV$

(112)

where $\mathcal{L}$ is the infinitesimal generator of $X_{t}$, such as introduced in equation (34). If we transport the equation (112) in the frequency domain, we infer that

$\begin{eqnarray} \mathcal{F}[\mathcal{L}V](\omega) & = & \left(i\, r\, \omega-\frac{1}{2}\omega^{2}\sigma^{2}+\int_{\mathbb{R}}\left(e^{i\omega z}-1\right)\mu^{J}(dz)\right)\mathcal{F}[V](\omega)\nonumber\\ & = & \psi(\omega)\mathcal{F}[V](\omega) \end{eqnarray}$

(113)

and that the Feynman-Kac equation becomes an ODE:

$\frac{\partial}{\partial t}\mathcal{F}[V](\omega)+\left(\psi(\omega)-r\right)\mathcal{F}[V](\omega)=0.$

(114)

$\mathcal{F}[V](\omega)$ is solution of this ODE and the option price is obtained by inversion, as stated in proposition 5.1.

If $X_{t}$ is a mean reverting CMLJD process, the Fourier's transform of the Feynman Kac equation is now

$\frac{\partial}{\partial t}\mathcal{F}[V]+a\omega\frac{\partial}{\partial\omega}\mathcal{F}[V]+a\mathcal{F}[V]-\frac{\sigma^{2}}{2}\omega^{2}\mathcal{F}[V]+\lambda\int_{\mathbb{R}}\left(e^{i\omega z}-1\right)\mu^{J}(dz)\mathcal{F}[V]=r\mathcal{F}[V]$

(115)

Change of variables $\gamma=e^{-a(T-t)}\omega$ or $\omega=e^{a(T-t)}\gamma$, then

$\begin{eqnarray} \frac{\partial}{\partial t}\mathcal{F}[V](\omega) & = & \frac{\partial}{\partial t}\mathcal{F}[V](e^{a(T-t)}\gamma)\nonumber\\ & = & -ae^{a(T-t)}\gamma\frac{\partial}{\partial\gamma}\mathcal{F}[V](\gamma)+\frac{\partial}{\partial t}\mathcal{F}[V](\gamma) \end{eqnarray}$

(116)

and the Feynman Kac Equation becomes

$\frac{\partial}{\partial t}\mathcal{F}[V](\gamma)+\left(a-r-\frac{\sigma^{2}}{2}\left(e^{a(T-t)}\gamma\right)^{2}+\lambda\int_{\mathbb{R}}\left(e^{ie^{a(T-t)}\gamma z}-1\right)\mu^{J}(dz)\right)\mathcal{F}[V](\gamma)=0$

(117)

then

$\mathcal{F}[V](t, \gamma) = \exp\left(-(r-a)(T-t)+\int_{0}^{T-t}\tilde{\psi}(ie^{a(T-t-u)}\gamma)du\right)\mathcal{F}[V](T, \gamma)$

(118)

and this proves proposition 5.2.

---

Appendix B

The next Matlab function implements the algorithm of proposition 4.1 that computes the probability density function of the CMLJD process.

function [x, fXt]=DensityXt(T, mu, sig, lam0, b1, alp0, b2, K, N, x_max)

% T: time horizon

% mu, sig, lam0, b1, alp0, b2, K : parameters defining the CMLJD

% x_max : maximum value for the domain of Xt

% N : number of DFT steps

     % definition of the grid

     x_min = -x_max;

     dx = (x_max-x_min)/(N-1);

     x = [ x_min:dx:x_max, ]';

     dz = 2*pi/(N*dx);

     z = [ 0:dz:(N-1)*dz]';

     % optimization the speed of the DFT

     fftw('planner', 'estimate');

     % evaluation of characteristic exponent and function

     psi = CharacteristicExponent(z, mu, sig, lam0, alp0, b1, b2, K, T);

     J = [1:N]';

     phi = (-1).^(J-1).*exp(psi.*T); % characteristic function

     phi(1) = 0.5* phi(1);

     % inversion by DFT of the characteristic function

     fXtraw = 2/((N)*dx)*real(fft(phi));

     % parallel shift to avoid negative value for the pdf

     fXtraw = fXtraw -min(fXtraw);

     adj = sum(fXtraw*dx);

     % ensure that the integral of the pdf is 1

     fXt = fXtraw./adj;

end

The following Matlab function evaluates by direct integration, the characteristic exponent of the CMLJD process, such as presented in equation (35).

function [psi] = CharacteristicExponent(z, mu, sig, lam0, alp0, b1, b2, K, T)

% This function calculates the Characteristic function by numerical

% integration of the function a_k^2/(a_k^2+z^2), such as introduced

% in proposition 2.1

     lam = lam0/(1+b1)*K^(b1+1)-lam0/(1+b1);

     fun = @(k) k.^(2*b2+b1)./(k.^(2*b2)+ (z./alp0).^2);

     q = integral(fun, 1, K, 'ArrayValued', true);

     psi = (T*(mu*i.*z -0.5*sig^2.*z.^2-lam+lam0*q));

end

This last function computes the integral $\int_{0}^{t}\tilde{\psi}\left(ize^{-a(t-u)}\right)du$, that is present in the characteristic exponent of the mean reverting CMLJD, equation (68).

function [logM] = CharacExponentMeanRevert(z, t, sig, lam0, alp0, b1, b2, K, a, b)

     lam = lam0/(1+b1)*K^(b1+1)-lam0/(1+b1);

     fun1 = @(k, u) lam0*k^b1*(alp0*k^b2)^2./((alp0*k^b2)^2-u.^2);

     fun2 = @(u) 0.5*sig^2*u.^2+integral(@(k) fun1(k, u), 1, K, 'ArrayValued', ...

         true)-lam;

     fun3 = @(u) fun2(1i.*z.*exp(-a.*(t-u)));

     fun4 = @(z) 1i.*z.*b*(1-exp(-a*t))+...

         integral(fun3, 1e-5, t, 'ArrayValued', true);

     logM = fun4(z);

end

---

Acknowledgment

I thank for its support the Chair "Data Analytics and Models for insurance" of BNP Paribas Cardiff, hosted by ISFA (Université Claude Bernard, Lyon France).

---

Conflict of Interest

All authors declare no conflict of interest in this paper.

---