Likelihood Evaluation of Jump-Diffusion Models Using Deterministic Nonlinear Filters

TL;DR

This paper introduces a deterministic nonlinear filtering approach to efficiently evaluate likelihood functions for jump-diffusion models with stochastic volatility and jumps, outperforming particle filters in speed and smoothness.

Contribution

The authors develop a high-dimensional deterministic filtering algorithm for complex jump-diffusion models, providing a faster and more precise alternative to particle filters.

Findings

Deterministic filtering is more accurate and faster than particle filtering.

Maximum likelihood estimates reveal increased jump intensity during the Great Recession.

Jump components significantly contribute to volatility clustering and negative returns.

Abstract

In this study, we develop a deterministic nonlinear filtering algorithm based on a high-dimensional version of Kitagawa (1987) to evaluate the likelihood function of models that allow for stochastic volatility and jumps whose arrival intensity is also stochastic. We show numerically that the deterministic filtering method is precise and much faster than the particle filter, in addition to yielding a smooth function over the parameter space. We then find the maximum likelihood estimates of various models that include stochastic volatility, jumps in the returns and variance, and also stochastic jump arrival intensity with the S&P 500 daily returns. During the Great Recession, the jump arrival intensity increases significantly and contributes to the clustering of volatility and negative returns.

Equations

$$\frac{d{S}_t}{S_{t^{-}}}=$$

$$d{V}_t=$$

$$d{\Lambda }_t=$$

$$y_t=$$

$$v_t=$$

$${\lambda }_t=$$

$${\epsilon }_t^y=\sqrt{1-{\rho }_v^2-{\rho }_{\lambda }^2}{\epsilon }_t^{\perp }+{\rho }_v{\epsilon }_t^v+{\rho }_{\lambda }{\epsilon }_t^{\lambda }$$

$$\mathcal{L}(\Theta )=$$

$$f\left({\mathbf{x}}_t|y_{1:t-1}\right)=$$

$$u\left({\mathbf{x}}_t|y_{1:t}\right)=\frac{r\left(y_t|{\mathbf{x}}_t,y_{1:t-1}\right)f\left({\mathbf{x}}_t|y_{1:t-1}\right)}{f\left(y_t|y_{1:t-1}\right)},$$

$$f\left(y_t|y_{1:t-1}\right)=$$

$$=$$

$$u\left({\mathbf{x}}_t|y_{1:t}\right)=$$

$$=$$

$$j_t^y=\sum _{i=1}^{n_t}{z}_{t,i}^y\text{and}{j}_t^v=\sum _{i=1}^{n_t}{z}_{t,i}^v,$$

$$q\left({\mathbf{x}}_t|{\mathbf{x}}_{t-1},y_{1:t-1},n_t=n\right)=$$

$$=$$

$$r\left(y_t|{\mathbf{x}}_t,{\mathbf{x}}_{t-1},y_{1:t-1},n_t=n\right)=r\left(y_t|v_t,v_{t-1},{\lambda }_t,{\lambda }_{t-1},j_t^y,j_t^v,n_t=n\right).$$

$$f\left(y_t|y_{1:t-1}\right)=$$

$$=$$

$$u\left({\mathbf{x}}_{t-1}|y_{1:t-1}\right)\mathbb{P}\left(n_t=n|{\mathbf{x}}_t,{\mathbf{x}}_{t-1}\right)d{\mathbf{x}}_{t}d{\mathbf{x}}_{t-1}$$

$$=$$

$$q\left(j_t^y|j_t^v,n_t=n\right)q\left(v_t|v_{t-1},j_t^v\right)q\left(j_t^v|n_t=n\right)q\left({\lambda }_t|{\lambda }_{t-1}\right)$$

$$u\left(v_{t-1},{\lambda }_{t-1}|y_{1:t-1}\right)d{\mathbf{x}}_{t}d{v}_{t-1}d{\lambda }_{t-1}$$

$${\iint }_{{\mathbb{R}}_{+}\times \mathbb{R}}u\left(v_{t-1},{\lambda }_{t-1},j_{t-1}^y,j_{t-1}^v|y_{1:t-1}\right)d{j}_{t-1}^{y}d{j}_{t-1}^v=u\left(v_{t-1},{\lambda }_{t-1}|y_{1:t-1}\right).$$

$$f\left(y_t|y_{1:t-1}\right)=$$

$$q\left(v_t|v_{t-1},j_t^v\right)q\left(j_t^v|n_t=n\right)q\left({\lambda }_t|{\lambda }_{t-1}\right)$$

$$u\left(v_{t-1},{\lambda }_{t-1}|y_{1:t-1}\right)d{v}_{t}d{\lambda }_{t}d{j}_t^{v}d{v}_{t-1}d{\lambda }_{t-1},$$

$${\int }_{\mathbb{R}}r\left(y_t|v_t,v_{t-1},{\lambda }_t,{\lambda }_{t-1},j_t^y,j_t^v,n_t=n\right)q\left(j_t^y|n_t=n\right)d{j}_t^y=$$

$$y_t|v_t,v_{t-1},{\lambda }_t,{\lambda }_{t-1},j_t^v,n_t=n\sim \mathcal{N}\left({\mu }_t,{\sigma }_t^2\right),$$

$${\mu }_t=$$

$${\sigma }_t^2=$$

$${\epsilon }_t^v=\frac{v_t-v_{t-1}-\kappa (\theta -v_{t-1})h-j_t^v}{\sigma \sqrt{v_{t-1}h}}\text{and}{\epsilon }_t^{\lambda }=\frac{{\lambda }_t-{\lambda }_{t-1}-\chi (\theta -{\lambda }_{t-1})h}{\xi \sqrt{{\lambda }_{t-1}h}}.$$

$$n_t|$$

$$v_t|$$

$$j_t^v|$$

$${\lambda }_t|$$

$$u\left({\mathbf{x}}_t|y_{1:t}\right)=u\left(v_t,{\lambda }_t|y_{1:t}\right)$$

Full text

Keywords: Discrete nonlinear filtering; Particle filter; Sequential importance resampling; Stochastic volatility; Maximum likelihood estimation.

1   Introduction and Motivation

Over the last thirty years, jump-diffusion models (with stochastic volatility) have become increasingly popular for their ability to replicate important stylized facts such as heavy tails, no autocorrelation in the returns, volatility clustering, etc. (Cont, 2007). Yet, parameter estimation of these models is cumbersome as stochastic volatility and jumps are latent—or not directly observed. Therefore, the first studies that looked into that problem focused mostly on stochastic volatility models using the (generalized) method of moments, quasi-, simulated and approximate maximum likelihood, and Bayesian methods.333The first study looking into parameter estimation of stochastic volatility models is Taylor (1986) which was based upon the method of moments. Generalizations of the method of moments were proposed by Melino and Turnbull (1990), Andersen and Sørensen (1996), Pan (2002). Quasi-maximum likelihood approaches were suggested by Nelson (1988) and Harvey et al. (1994); Danielsson (1994) and Brandt and Santa-Clara (2002) used simulated maximum likelihood whereas approximate maximum likelihood were introduced by Aït-Sahalia (2002), Bates (2006) and Aït-Sahalia and Kimmel (2007). Finally, another family of approaches based upon the Bayesian paradigm was used in this literature, namely Shephard (1993), Jacquier et al. (1994), Johannes et al. (1999) and Eraker (2001).

Most recent progress in parameter estimation of jump-diffusion models has been accomplished with numerical recursive prediction-update algorithms, namely sequential Monte Carlo (SMC) methods and discrete nonlinear filters (DNF). Such methods provide the posterior distribution of latent state variables conditional on current and past observations in a recursive manner.444A classic example of prediction-update algorithms is the Kalman (1960) filter which provides the exact posterior distribution of the state variables in a Gaussian and linear framework. As a by-product of the method, one can obtain the likelihood function. On the one hand, sequential Monte Carlo (SMC) methods—so-called particle filters—as introduced by Gordon et al. (1993), have been used to compute and maximize the likelihood function (Johannes et al., 2009; Christoffersen et al., 2010; Pitt et al., 2014; Bégin et al., 2019; Bardgett et al., 2019). Although the SMC method is very flexible, it is Monte Carlo-based and thus computationally intensive. Moreover, the likelihood function is random and not smooth for a finite set of particles, which can be a significant issue for frequentist-based inference.555Malik and Pitt (2011) introduced a resampling algorithm that smoothes the likelihood function but according to Creal (2012), this only works when the state dimension is one, which is typically the case when only stochastic volatility is included whereas jumps are excluded.

On the other hand, one can approximate the prediction density as well as the likelihood function with numerical integration schemes; this is the framework put forward by Kitagawa (1987). As such, the dimension of the integration problem increases with the number of state variables which explains why many authors solely focused on stochastic volatility models (Fridman and Harris, 1998; Watanabe, 1999; Bartolucci and De Luca, 2001; Clements et al., 2006a; Langrock et al., 2012) whereas jumps were ignored.

In this article, we present a DNF-based method for the likelihood evaluation of stochastic volatility models that include return jumps, variance jumps, as well as stochastic jump arrival intensity. Those include well-known models such as the stochastic volatility (SV) model, stochastic volatility with return jumps (SVYJ) model, and the stochastic volatility with return and variance correlated jumps (SVCJ) model, which are nested cases of the general model. We then assess the accuracy of the proposed methodology and compare it with the most common SMC method—the sequential importance resampling (SIR) of Gordon et al. (1993). Finally, we apply the method to a typical financial time series: the S&P 500 index.

Our main contributions and findings are as follows. We first develop a prediction-update algorithm for the SVCJSI model—and its nested models. The proposed approach reduces the dimensionality of the resulting integrals when compared to a naive application of Kitagawa (1987). As a result, the proposed DNF method is both accurate and faster than the SIR. And because the DNF yields a smooth likelihood function in the parameter space, it is ideal for numerical maximization. We thus compute the maximum likelihood estimates of the SV, SVYJ, SVCJ and SVCJSI666SVCJ with stochastic jump arrival intensity models using S&P 500 daily returns. We find that during the Great Recession, the jump arrival intensity increases significantly and contributes to the clustering of volatility and negative returns.

The article is organized as follows. Section 2 presents the econometric framework. The DNF is explained in Section 3. Section 4 discusses the accuracy and the reliability of the DNF method. The speed-accuracy trade-off is assessed in Section 5. An empirical application is presented in Section 6. Finally, some concluding remarks are provided in Section 7.

2   Framework

This section describes the continuous- and discrete-time frameworks. We then specifically introduce the four models that are used throughout this paper, namely the SV, SVYJ, SVCJ, and the full model with stochastic jump arrival intensity, known as SVCJSI.

2.1   Continuous-Time Framework

We fix a filtered probability space $(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ and a filtration $\mathbb{F}=\{\mathcal{F}_{t}\,:\,t\geq 0\}$ satisfying the usual conditions. Let $S_{t}$ be the time-$t$ price of a security, $V_{t}$ the time-$t$ instantaneous variance and $\Lambda_{t}$ the time-$t$ jump intensity. Under the objective measure $\mathbb{P}$, their dynamics are given by the following equations:

[TABLE]

where $W^{S}$, $W^{V}$ and $W^{\Lambda}$ are three Brownian motions such that $d\langle W^{S},W^{V}\rangle_{t}=\rho_{v}\,dt$, $d\langle W^{S},W^{\Lambda}\rangle_{t}=\rho_{\lambda}\,dt$ and $d\langle W^{V},W^{\Lambda}\rangle_{t}=0\,dt$. Moreover, $\{N_{t}\}_{\{t\geq 0\}}$ is a Cox process (doubly stochastic Poisson process) with jump arrival intensity $\Lambda_{t}$ and jump sizes to be further discussed. Note that only $S_{t}$ is observed whereas $V_{t}$ and $\Lambda_{t}$ are unobserved or latent.

The first stochastic differential equation (SDE) resembles that proposed by Merton (1976) as it includes both a jump process as well as a diffusive component. Yet, in opposition to Merton, our diffusive component allows for stochastic volatility—modelled by the second SDE. The stochastic variance SDE is given by the square-root process—similar to that used by Heston (1993)—to which variance jumps are added. Parameter $\kappa$ denotes the speed of mean reversion, $\theta$ is the unconditional mean variance, and $\sigma$ is the so-called variance of the variance parameter. Finally, the last SDE models the stochastic jump arrival intensity; parameter $\chi$ is related to the speed of mean reversion, $\omega$ is the long run level of the stochastic jump intensity and $\xi$ is the variance of the intensity parameter.

We assume that each of the return jumps $\{Z_{n}^{S}\}_{\{n\in\mathbb{N}\}}$ are normally distributed with a mean of $\alpha+\rho_{z}Z_{n}^{V}$ and a standard deviation of $\delta$, i.e., $Z_{n}^{S}\sim\mathcal{N}(\alpha+\rho_{z}Z_{n}^{V},\delta^{2})$ and that the variance jumps $\{Z_{n}^{V}\}_{\{n\in\mathbb{N}\}}$ are exponentially distributed with a mean of $\nu$, i.e., $Z_{n}^{V}\sim\text{Exp}(\nu)$. Similar return jump innovations have been used by Bates (1996), Bakshi et al. (1997), Duffie et al. (2000), Pan (2002), Eraker et al. (2003), and Johannes et al. (2009) in continuous-time models whereas Maheu and McCurdy (2004), Christoffersen et al. (2012) and Ornthanalai (2014) (among others) have introduced normally distributed jumps in discrete-time models. Evidence for the presence of positive jumps in the variance process has been found by numerous authors, namely Bates (2000), Duffie et al. (2000), Pan (2002), Eraker et al. (2003) and Todorov and Tauchen (2011).

Finally, the term $\Lambda_{t}\overline{\alpha}$ is the jump compensator that makes sure that $\mathbb{E}^{\mathbb{P}}\left[S_{T}\,\middle|\,\mathcal{F}_{t}\right]=S_{t}e^{\,\mu(T-t)}$, with $\mu$ being the annual expected rate of return on the stock in our setting. We thus have $\overline{\alpha}={\exp\left(\alpha+\frac{1}{2}\delta^{2}\right)}/{\left(1-\rho_{z}\nu\right)}-1$.

2.2   Discrete-Time Framework

As we typically observe most continuous-time processes only in discrete-time, we thus need to discretize our continuous-time dynamics. Setting $h$ as the length of a time interval (say e.g. a month, a week or a calendar/trading day), then we let $y_{t}\equiv\log(S_{th}/S_{(t-1)h})$ be the time-$t$ (continuously compounded ex-dividend) return of a security. As such, the information set available to the econometrician $\mathbb{G}=\{\mathcal{G}_{t}\,:\,t\in\mathbb{N}\}$ is coarser than that of the model, i.e., $\mathcal{G}_{t}\subseteq\mathcal{F}_{t}$ since $\mathcal{G}_{t}=\sigma\left(\{y_{s}\}_{s=1}^{t}\right)$. Similarly, we define $v_{t}$ and $\lambda_{t}$ as the time-$th$ discretized variance and jump intensity processes, respectively. Again, as in the continuous-time version of the model, $v_{t}$ and $\lambda_{t}$ are latent.

The Euler-Maruyama discretization scheme applied to the continuous-time version of the model yields

[TABLE]

where $n_{t}=N_{th}-N_{(t-1)h}$ is given by a Poisson random variable with parameter $\lambda_{t-1}h$. Moreover, the standardized innovations are normally distributed, i.e., $\varepsilon_{t}^{i}\sim\mathcal{N}(0,1)$ for $i\in\{\perp,v,\lambda\}$ with

[TABLE]

obtained from the Cholesky decomposition. Finally, return jumps are also normally distributed, i.e. $z^{y}_{t,i}\sim\mathcal{N}\left(\alpha+\rho z_{t,i}^{v},\,\delta^{2}\right)$, whereas the distribution of variance jumps is still exponential, i.e., $z^{v}_{t,i}\sim\text{Exp}(\nu)$.

2.3   Nested Models

This general specification embeds numerous stochastic volatility-type models. The four specific models considered in this study are detailed below:

• –

SV: Stochastic volatility model with no jumps, obtained by letting $\chi=\omega=\xi=0$. It is similar to the model proposed by Heston (1993) in the sense that it incorporates a leverage effect. Specifications alike were also investigated by Melino and Turnbull (1990), Danielsson and Richard (1993), Shephard (1993), and Jacquier et al. (1994), among others.

• –

SVYJ: Stochastic volatility model with jumps in return only, obtained by letting $\chi=1$ and $\xi=\nu=0$. It resembles the Bates’s (1996) model. In econometrics, Pitt et al. (2014) proposed the use of jumps in the return’s dynamics.

• –

SVCJ: Stochastic volatility with simultaneous and correlated jumps in return and variance, obtained by letting $\chi=1$ and $\xi=0$. The specification was first used by Duffie et al. (2000).

• –

SVCJSI: Stochastic volatility with simultaneous and correlated jumps in return and variance along with stochastic jump arrival intensity. This is the most complicated and general specification used in this study.

3   Discrete Nonlinear Filtering

Because the volatility, jumps and the jump intensity factors are latent variables in the model, it is difficult to evaluate the likelihood function $\mathcal{L}$ in a direct manner. As such, we propose in this section a deterministic DNF scheme that uses a recursive prediction-update algorithm in the spirit of Kitagawa (1987). This methodology is based on a discretization of the state-space continuous latent variables, similar to the method proposed by Fridman and Harris (1998), Bartolucci and De Luca (2001), Clements et al. (2006b), and Langrock et al. (2012). It is important to note that such a methodology has not been applied to the likelihood evaluation of SVYJ, SVCJ and SVCJSI models, which is a key contribution of this paper.

Let $\mathbf{x}_{t}\in\mathcal{X}$ be the latent state variables at time $t$ where $\mathbf{x}_{t}$ is either a scalar or a vector. Moreover, let us define by ${y}_{t},t=1,2,\dots,T$ the time-$t$ observations—i.e., returns in this study. Series of latent factors and observations are denoted by $\mathbf{x}_{0:t}\equiv\{\mathbf{x}_{s}\}_{s=0}^{t}$ and $y_{1:t}\equiv\{y_{s}\}_{s=1}^{t}$. For a parameter set $\Theta$, the likelihood function of $\Theta$ is

[TABLE]

which involves a multidimensional integration problem.

3.1   General Methodology

The one-step ahead prediction of the latent states’ distribution (conditional on the past returns) is given by

[TABLE]

where $q$ is the transition probability distribution of $\mathbf{x}_{t}$ given $\mathbf{x}_{t-1}$ and the past returns. Once a new return is available, the posterior distribution of the state variables at time $t$, conditional on $\mathcal{G}_{t}$, may now be obtained as

[TABLE]

where $r$ is the measurement probability density—or the conditional time-$t$ return density given the state space variables $\mathbf{x}_{t}$. This is also known as the update step of the algorithm. The denominator of Equation (2) is the time-$t$ likelihood contribution of return $y_{t}$ conditional on the past returns, i.e.,

[TABLE]

Given the model formulation of Section 2, $\mathbf{x}_{t}$ depends upon $\mathbf{x}_{t-1}$. Hence, the time-$t$ likelihood contribution can be rewritten as

[TABLE]

where $u\left(\mathbf{x}_{t-1}\,\middle|\,y_{1:t-1}\right)$ is similar to Equation (2). Then, the posterior density of Equation (2) can also be obtained recursively using the same idea:

[TABLE]

Therefore, to obtain the likelihood function of Equation (1), one needs to find the posterior distribution of the latent states and the likelihood contribution at all times. Algorithm 1 summarizes the various steps of the general DNF.

3.2   Likelihood Function of the SVCJSI Model

We seek to derive the likelihood function of the SVCJSI model using the DNF method.777Being nested models, the SV, SVYJ and SVCJ models could also be estimated with the same method by reducing the problem’s dimension. Specifically, for the framework of Section 2, we have that \mathbf{x}_{t}=\big{[}\,v_{t}\quad\lambda_{t}\quad j_{t}^{y}\quad j_{t}^{\vphantom{y}v}\,\big{]}, where

[TABLE]

as knowing every individual jump sizes over a day is futile—only the aggregated jump is important in the discretization. Therefore, $\mathcal{X}=\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}\times\mathbb{R}_{+}$ and the latent state vector is four-dimensional for the SVCJSI model.

Naively computing Equations (3) and (4) might be cumbersome given the high-dimensionality of the state variables. Fortunately, some simplifications can be applied to reduce the dimension of the problem. First, by conditioning on the number of jumps at time $t$, we can simplify the transition density in the following way:

[TABLE]

as $v_{t}$ and $\lambda_{t}$ depend upon $v_{t-1}$ and $\lambda_{t-1}$, respectively. Yet, $j_{t}^{(y)}$ and $j_{t}^{(v)}$ do not depend on the previous value of the latent state; in fact, return and variance jump sizes are (iid) transient latent states. By conditioning again on the number of jumps at time $t$, the observation density is given by:

[TABLE]

Then, following the rationale outlined at the beginning of Section 3, the time-$t$ likelihood contribution of Equation (3) can be rewritten as

[TABLE]

because

[TABLE]

Conditioning on the number of jumps also allows us to further simplify Equation (5) so that

[TABLE]

as

[TABLE]

which has a normal density. Specifically, we have that

[TABLE]

where

[TABLE]

[TABLE]

We can obtain the integrand of Equation (6) in closed form because

[TABLE]

and the corresponding densities are known in closed form as well. Note that $\mathcal{P}(m)$ is a Poisson distribution with mean $m$ and $\Gamma(n,\nu)$ is a gamma distribution with mean $n\nu$.

Finally, one can get the time-$t$ posterior density of the latent states in a similar fashion by computing the following integral

[TABLE]

3.3   Numerical Implementation

To be able to numerically compute the integrals above, we need to define a discretization of the state space—so-called likely sequences. The interval boundaries used to discretize the state space of the variance factor $v$, the jump intensity factor $\lambda$ and the variance jumps $j$, are chosen over the following ranges:

[TABLE]

where $\mathcal{E}^{V}_{\infty}$ and $\mathcal{V}^{V}_{\infty}$ are the long-run expected value and variance of the stationary process of the variance factor, respectively. Moreover, $\mathcal{E}^{\Lambda}_{\infty}$, $\mathcal{V}^{\Lambda}_{\infty}$, $\mathcal{E}^{J}_{\infty}$ and $\mathcal{V}^{J}_{\infty}$ are defined analogously for the jump intensity factor and for the variance jumps, respectively. The intermediate points are determined by the following equations:

[TABLE]

making the variance (jump intensity) likely sequences uniformly distributed in the volatility (square root of the jump intensity) domain.

The functions $\delta^{V}_{N}$, $\delta^{\Lambda}_{M}$ and $\delta^{J}_{K}$ assure that the DNF converges and shall therefore respect two sets of conditions:

[TABLE]

The first set of conditions makes sure that we cover the domain as the number of nodes increases. On the other hand, the second set of conditions ensures that the partition becomes finer and finer as the number of nodes increases. Obviously, a myriad of such functions exists. In this study, we use $\delta_{N}^{V}=3+\log(N)$, $\delta_{M}^{\Lambda}=3+\log(M)$ and $\delta_{K}^{J}=3+\log(K)$, although other functions could be used (as long as they satisfy the two sets of conditions).

Moreover, in the spirit of Langrock et al. (2012), a special integration rule will be used in this study: broadly speaking, if $g_{1}$ and $g_{2}$ are two functions, then

[TABLE]

where $c$ is a representative point in $[a,b)$ (e.g., the midpoint). Given that we are integrating products of probability density functions (pdfs), this integral rule will allow us to rewrite the integral of Equation (7) as $g_{1}(c)\left(G_{2}(b)-G_{2}(a)\right)$, where $G_{2}$ is the corresponding cumulative distribution function of the pdf $g_{2}$.888Other quadratures have been proposed in the literature. For instance, Fridman and Harris (1998) use of Gauss-Legendre quadrature. Bartolucci and De Luca (2001) and Clements et al. (2006a) use a simple numerical integral scheme based on the midpoint $c$, i.e.,

$\int_{a}^{b}f_{1}(x)g_{2}(x)\,dx\approx(b-a)g_{1}(c)g_{2}(c).$

Then, to compute numerical integrals, it is often more convenient to define intervals. Therefore, let

[TABLE]

be three different sets of intervals that cover $[0,\infty)$, where $v^{(0)}=-v^{(1)}$, $\lambda^{(0)}=-\lambda^{(1)}$, $j^{(0)}=-j^{(1)}$, $v^{(N+1)}=\infty$, $\lambda^{(M+1)}=\infty$, and $j^{(K+1)}=\infty$.

Finally, the time-$t$ likelihood contribution is approximated by

[TABLE]

where $R$ is the truncation level of the Poisson random variable. The time-$t$ posterior density $u\left(v_{t},\lambda_{t}\,\middle|\,y_{1:t}\right)$ is also estimated in a similar way:

[TABLE]

4   Accuracy and Reliability

We now aim to assess the accuracy and reliability of the DNF. Specifically, we first verify that it converges to the correct likelihood value by comparing the log-likelihood obtained with the DNF to that of the SIR assuming that a large computing budget is available. This experiment will tell us how the DNF performs in a best-case scenario.

As explained in Malik and Pitt (2011), the SIR-based likelihood function is unbiased (see Section A of the Supplementary Material for more details on the SIR method). Therefore, when a large computing budget is used, the SIR-based likelihood estimate should be very close to its true value with very little sampling error. Hence, it can be used as a reliable benchmark to assess the accuracy and reliability of the DNF.

For the SIR method, we use an increasing number of particles as a function of the model complexity. Specifically, we use 100,000 particles for the SV model, 250,000 particles for the SVYJ model and 1,000,000 particles for both the SVCJ and SVCJSI models. Then, we set the number of nodes in the variance (variance jump) grid to $N=200$ ($K=100$) for the SV, SVYJ and SVCJ models.999According to our results, a reliable rule of thumb is to set the number of nodes of the variance jump grid, $K$, to be two and three times fewer than the number of points of the variance. As for the SVCJSI model, we set $N=M=50$ and $K=25$ due to memory size limits.101010We used a graphics processing unit (GPU) to compute the likelihood function of the SVCJSI model with the DNF; see Section 5.1.

4.1   Random Series

In this first test, we compare the likelihood of 1,000 one-year daily series taken at random using parameters that are also randomly generated. The aim of this test is to determine the accuracy of the DNF for a wide range of possible paths and parameter sets.

Therefore, we generate parameter sets using the following bounds:

[TABLE]

These values are reasonable bounds for the parameters; they span multiple potential and realistic cases. Then, based on the simulated parameters, we generate one-year paths and we compute the log-likelihood using both methods. Specifically, we compare both estimates using the absolute percentage error (APE):

[TABLE]

Table 4.1 reports seven quantiles of the APE distribution for the four models considered. Broadly speaking, the medians hover between 0.01% and 0.05% for the four models which is very small. The right tail of the APE distributions is also very thin with 99.5${}^{\text{th}}$ quantiles that are between 0.40% and 0.85%.

Figure 1 supplements Table 4.1 by reporting scatter plots of SIR- and DNF-based log-likelihoods. Most of these are aligned on the diagonal meaning that both methods are yielding very much alike estimates when using a large computing budget. Based upon this first test, we can comfortably conclude that the DNF is very accurate, at least within the range of parameters provided above and for a large computing budget.

4.2   S&P 500 Index

Instead of using randomly generated series that might not resemble those used in empirical finance, we now focus on the S&P 500 returns as these would be representative of typical asset returns. Specifically, we use daily S&P 500 index returns (excluding dividends), from January 1990 to September 2018, obtained from the Bloomberg terminal. For this specific time series, we then randomly generate a set of parameters and compute the resulting likelihood with both the SIR and DNF methods.

Table 4.2 mimics the results of Table 4.1, only this time we consider the S&P 500 returns over the last three decades. The results are similar to those obtained with random series: in fact, the median APE is rather small, ranging between 0.02% and 0.10%. Again, even the worst cases are still very decent, with the 99.5${}^{\text{th}}$ quantile being below 1% for the four models. Therefore, for a typical financial time series and a large computing budget, the DNF is very accurate.

5   Precision and Computing Times

As shown in the previous section, the DNF is capable of providing very accurate log-likelihood estimates when a large computing budget is available. Nonetheless, it is of great interest to also assess its performance when a smaller budget is considered because most inference methods—both frequentist and Bayesian—require a considerable number of likelihood function evaluations. The overarching goal of this section is to analyze the precision of the DNF compared with computing times required to attain such precision.

We first discuss the use of central processing unit (CPU) and graphics processing unit (GPU) for the problem at hand. We then compare the DNF and SIR in terms of computing times and precision. Finally, we analyze the speed and accuracy trade-off of the DNF by determining the grid size or computing time necessary to achieve a given level of precision (and vice versa).

5.1   Motivation for CPU and GPU Computing

The computing times involved in this study are greatly influenced by the number of state variables considered, which has an impact on the dimension of the sum of Equation (8). Indeed, Equation (8) involves a sum over two, three, four or six dimensions for the SV, SVYJ, SVCJ and SVCJSI models, respectively. Obviously, such calculations can become cumbersome when $N$, $M$, $K$, $R$ and $T$ are large.

It is easy to see why a GPU can be very useful. For instance, if $N=50$, $M=50$, $K=25$ and $R=2$, then the size of the six-dimensional matrix containing all the elements that need to be summed up is

[TABLE]

Moreover, the number of operations needed to obtain one likelihood contribution is rather large and it increases exponentially, e.g., $5\times 50\times 50\times 50\times 50\times 25\times 2=1{,}562{,}500{,}000$ multiplications, and $50\times 50\times 50\times 50\times 25\times 2=312{,}500{,}000$ additions for the example above.

CUDA-enabled GPUs are extremely efficient when it comes to calculating sums. In fact, for very large matrices, the time needed to compute sums is manifold lower when using a GPU instead of a CPU. Figure 2 shows the computing time for summing all the entries of an $n\times n$ matrix while using a GPU (black) or a CPU (grey). For instance, the GPU-based calculation is eight times faster for a matrix of size 20,000 $\times$ 20,000—which is a significant decrease.

But for GPU computing to cut overall computing times, it is important to take into account the additional time required to move data between the CPU and GPU, known as the overhead. For array sizes involved with deterministic jump arrival intensity models, our experiments showed that a conventional CPU is still faster overall.

Therefore and unless stated otherwise, the SIR and DNF methods have been coded in Matlab 2018b on a single thread with a typical CPU. However, the DNF method applied to the SVCJSI model has been implemented with an NVIDIA QUADRO P6000 GPU, also using Matlab 2018b.111111The computing times are obtained by running Matlab programs on a desktop computer with two 2.1 GHz Intel Xeon E5-2620 v4 and 48 GB RAM. The NVIDIA QUADRO P6000 we use has 3840 CUDA parallel-processing cores and a GPU memory of 24 GB.

5.2   Deterministic Jump Arrival Intensity: SV, SVYJ and SVCJ Models

Since models with deterministic and stochastic jump arrival intensity have been implemented differently, we analyze computing times and precision for both families of models separately.

5.2.1   Computing Times: Comparison Between DNF and SIR

The goal of this first test is to analyze the accuracy of SIR and DNF for different fixed and finite computing budgets. For the SIR method, the computing budget is determined by the number of particles whereas for the DNF method, the grid size parameters $N$ and $K$ determine the computing budget. For all computations, we use the same daily S&P 500 return series over 5 years (about 1250 observations) evaluated with the parameters that produced the highest likelihood in Section 4.2 as a proxy for maximum likelihood parameters.

Then, for a given computing budget, we compute the likelihood value and record the computing time. Moreover, for the SIR, we repeat the exercise 50 times to approximate the distribution of the log-likelihood values. This entire exercise is then rerun for several computing budgets.

Figure 3 reports this comparison. We show boxplots (in grey) of the distribution of the SIR-based log-likelihood values and we compare them to values obtained with the DNF method (black circles). The true value (dashed line) is given by the average obtained with SIR using a very large computing budget (this number is virtually the same as that obtained with the DNF method for a large budget as well).

For the SV model, the DNF method yields values that are very close to the true log-likelihood value—and this, after only 100 milliseconds (ms). For such small computing budget, the SIR does not produce adequate results: it is biased and shows important variations. For a budget of one second per likelihood valuation, the SIR still gives very noisy log-likelihood estimates even though this value is asymptotically consistent: the whiskers are quite long, but the median is aligned with the true value.

The DNF provides similar results for the SVYJ model: after 600 ms, the deterministic method produces very steady and stable estimates of the log-likelihood function—consistent with the true value obtained with a large computing budget. The SIR-based log-likelihood, on the other hand, is still volatile even with a budget of two seconds.

The results for the SVCJ model tell a similar story: after a few seconds, the DNF is precise and stable and quickly converges to the true value whereas the SIR-based log-likelihood still shows significant sampling variations as the whiskers are still clearly visible to the naked eye.

5.2.2   Speed and Accuracy Trade-Off

We know the DNF method performs very well compared to the SIR providing accurate results much quicker than the SIR. In this second test, we would like to determine the precision gained (or lost) if we have a limited computing budget, as defined by either the grid size or computing time.

In this experiment, we generate 100 sets of parameters and we evaluate the likelihood function for varying number of nodes. We use the number of variance nodes $N$ ranging from 25 to 200 with steps of 5. The number of variance jump nodes is consistent with the rule of thumb explained above, i.e., $N/2.5$. The number of nodes associated with the jumps, $R$, is set to 2. We still use daily S&P 500 returns over 5 years for about 1250 observations.

The likelihood values obtained for these varying number of nodes is then compared to that obtained with the SIR while using a large computing budget (see Section 4). We still use the APE as a basis of comparison and compute mean APE (MAPE) to assess the average difference between the log-likelihood values:

[TABLE]

The leftmost panels of Figure 4 report the smoothed version of the MAPE against size—here, proxied by the number of nodes in the variance grid $N$.121212We use Matlab’s fit function with a power model, i.e., $y=ax^{b}$, to fit the MAPE curve. The results are robust to other specifications. For the four models considered, a size between 50 and 60 would lead to a MAPE below 0.1%. This result is in fact consistent with the size used in Watanabe (1999), Clements et al. (2004), Clements et al. (2006a) and Langrock et al. (2012).

The rightmost panels of Figure 4 detail the relationship between the MAPE and the computing time per likelihood valuation (in seconds). Again, as explained above, the computing time tends to increase with the complexity of the model. It takes approximately 75 ms, 200 ms and 300 ms to attain an MAPE of about 0.1% with the SV, SVYJ and SVCJ models, respectively. Therefore, with a conventional CPU, one is able to compute the likelihood function in much less than a second and reach an MAPE of 0.1%.

5.3   Stochastic Jump Arrival Intensity: SVCJSI Model

We now analyze the performance of the DNF for the SVCJSI model. As opposed to the family of deterministic jump arrival intensity models, the SVCJSI model is comprised of four state variables and to accelerate computations, we implemented the DNF using a high-end GPU.

As in Section 5.2, we first seek to compare the precision and computing times of the SIR and DNF methods. Figure 5 is the same as Figure 3, but for the SVCJSI model.

For a very small computing budget (a couple of seconds), the SIR generates a very wide range of likelihood values. As we add particles and thus the required computing time, the width of whiskers decreases as expected. However, even if the computing budget is large (2 or 3 minutes per likelihood valuation), there is still a significant uncertainty that remains around the likelihood value obtained. But if we are equipped with a powerful GPU, one is capable of achieving a very precise likelihood value with the DNF within a minute.

Figure 5 also illustrates the difficulty of precisely evaluating the likelihood function for the SVCJSI model (with the technology available to most users in 2019). If the computing budget available is small (e.g., a couple of seconds), then one should turn to the SIR implemented on a CPU but the resulting likelihood value will be biased and random. If a GPU is available, the computing budget must be large enough to overcome the extra overhead: at least 30 to 60 seconds are needed to achieve a level of accuracy that is impossible to attain with the SIR with a finite computing budget.

We also analyze the MAPE as a function of the grid size or computing time. The number of variance and jump intensity nodes ranges from 25 to 50 with steps of 1, with $N=M$.131313A number of nodes of 50 is the physical size limit the GPU handles at the time of writing this article. Moreover, the number of variance jumps is set to $\lceil N/2.5\rceil$ and the number of jump nodes is set to 2 still.

The results are illustrated in Figure 6 for the SVCJSI model. Again, we observe that a grid size of about 50-60 achieves an MAPE of 0.1% (as with the SV, SVYJ and SVCJ models) but the required computing time to achieve such MAPE is about 160 seconds with GPU computing.

6   Empirical Applications

In the previous applications, we have shown the precision of the DNF in evaluating the likelihood function. For frequentist-based inference such as the maximum likelihood, the likelihood function must be smooth enough for an optimizer to find a global maximum in the parameter space. As such, the DNF has a key advantage over the SIR as its construction leads to a smooth likelihood function.

We thus present in this last section an important application of the DNF methodology: we compute the maximum likelihood estimates (MLE) of each of the four models studied in this paper. This section is certainly one of the few attempts in the literature to compute the MLE of such complex financial models.

We use the S&P 500 return series (excluding dividends) from 1990 to 2018 (see Figure 7). This time frame contains a few recessions and eras of turmoil: the early 1990s recession (subsequent to the 1990 oil price shock), the early 2000s recession (partly caused by the collapse of the speculative dot-com bubble and the 9/11 attacks) and the Great Recession (a by-product of the collapse of the US housing bubble, followed by a global financial crisis). This sample should therefore contain a time-varying volatility, return jumps, variance jumps, and potentially, a changing jump intensity. The grid sizes used for the computing of the likelihood function are such that $N=M=50$, $K=20$ and $R=2$, which are consistent with our results.

Table 6 reports the MLE for each model along with their standard error in brackets.141414Robust standard errors (in brackets) are calculated from the outer product of the gradient at the optimum parameter value. For more details on this method, see Appendix B.5 of Rémillard (2013). These estimates are comparable to those found in other studies and using different statistical methodologies (e.g., Andersen et al., 2002; Eraker et al., 2003; Christoffersen et al., 2010; Hurn et al., 2015).

Figure 8 shows the filtered annualized volatility (i.e., $\sqrt{V_{t}/252}$), the jump intensity, the return jumps and the variance jumps. Generally speaking, the volatility is rather similar between the various models, the only difference being that the SVCJ and SVCJSI models allow for sharper increases because of the volatility jumps. For instance, during the last financial crisis, the volatility goes higher with the models that permit variance jumps.

For models with constant jump intensity, the number of jumps per year is small—about three to five, on average. This figure is consistent with other studies such as Eraker et al. (2003) and Hurn et al. (2015). The filtered return jumps are rather similar for the SVYJ and SVCJ models.

On the other hand, the SVCJSI model allows for jump clustering as the jump intensity varies over time, and the jump processes are slightly different of those obtained with the SVYJ and SVCJ models. Specifically, the number of filtered jumps reduces during calm eras and increases during crises.

7   Concluding Remarks

In this paper, we developed a new estimation method for jump-diffusion models based upon the DNF method of Kitagawa (1987). The method we introduced allows for the likelihood evaluation of models that include stochastic volatility, return jumps, variance jumps, as well as stochastic jump intensity.

We then compared the performance of the DNF with the SIR. We found that the proposed DNF method is stable, accurate and faster than the SIR. We showed that a grid with 50 to 60 nodes leads to an MAPE of 0.1%, on average, for the S&P 500 index—a threshold that is reasonable. As an application of the method, we also found the maximum likelihood estimates of the SVCJSI model and its nested models with the S&P 500 daily returns. The parameters obtained are consistent with those estimated in the literature.

Finally, this method is not only helpful for maximum likelihood estimation, but also for Bayesian inference as the likelihood function is an important building block when calculating posterior distributions. For instance, it would be interesting to combine the DNF with Markov chain Monte Carlo (MCMC) methods; this is left for future research.

Appendix A Particle Filter

In this section, we briefly describe the sequential Monte Carlo method used to compute the likelihood function in this study. This implementation of the particle filter closely follows the work of Gordon et al. (1993); specifically, we use the bootstrap filter. Algorithm 2 describes the various steps of the general SIR method.

Appendix B Study of the DNF Bias

The efficacy of the DNF method for estimating the parameters of jump-diffusion models using maximum likelihood estimation is now assessed. In this experiment, we compare the parameter estimates obtained by using DNF with their respective true values. The primary objective of this simulation exercise is to assess whether the DNF leads to systematic biases in the estimation parameters.

Specifically, for a set of likely parameters—consistent with most studies cited in the introduction of this study—we generate 100 ten-year return paths, i.e., $T=2,520$. Then, we use the DNF method coupled with a numerical optimizer to find the maximum likelihood estimates.

Table B reports the average value, the bias as well as the root mean square error for each parameter and each model. In summary, the average is very close to the true value for most of the parameters, leading to small biases. The only exception is parameter $\kappa$: the bias associated with this parameter hovers between 0.25 and 0.75. Yet, even though this bias looks large at first sight, it is not so dreadful. In fact, when expressed in terms of daily persistence of the process—as it is typically done in econometrics—the difference is rather small, e.g., $e^{-3.717h}=0.9882$ versus $e^{-3.000h}=0.9854$ for the SV model.