Parameter Estimation for Multivariate Diffusion Systems

Diffusion processes are continuous-time, continuous-space stochastic processes that have proven to be natural modelling frameworks for many real world phenomena. Over an infinitesimal interval $dt$, the evolution of a multivariate diffusion process $\phi_{t}^{*}$ is represented by the following, possibly time inhomogeneous, stochastic differential equation (SDE):

where $\phi_{t}^{*}=\left(\phi_{i}\right)_{i=1,...,m}$; $\theta=\left(\theta_{i}\right)_{i=1,...,p}$ is the parameter vector; $\mu\left(\phi_{t}^{*},t;\theta\right)=\left(\mu_{i}\right)_{i=1,...,m}$; $\sigma^{2}\left(\phi_{t}^{*},t;\theta\right)=$ $\left(\sigma_{ij}\right)_{i,j=1,...,m}$ with $\sigma^{2}\left(\phi_{t}^{*},t;\theta\right)=\sigma\left(\phi_{t}^{*},t;\theta\right)^{T}\sigma\left(\phi_{t}^{*},t;\theta\right)$ and $B_{t}$ is an $m$-dimensional vector of independent Brownian motions. The $m$-dimensional vector $\phi_{t}^{*}$ represents a set of state variables which characterizes the diffusion system through time. The assumption that the $m$ Brownian motions are independent does not lead to any loss of generality since allowance is made for off-diagonal terms within the diffusion matrix $\sigma^{2}\left(\phi_{t}^{*},t;\theta\right)$. Within this framework, the primary focus is on estimating the parameter vector $\theta$ from discretely-sampled data.

The drift vector $\mu\left(\phi_{t}^{*},t;\theta\right)$ and the diffusion matrix $\sigma^{2}\left(\phi_{t}^{*},t;\theta\right)$ characterize the evolution of $\phi_{t}^{*}$. Their individual elements are defined as:

Over an infinitesimally small interval, the diffusion system is distributed as follows:

That is, the diffusion system has a multivariate-normal distribution characterized by the drift vector and diffusion matrix over any infinitesimal interval. Since diffusion processes are Markovian, equation (2) may be used to derive the likelihood for continuously sampled diffusion paths. However with discretely sampled diffusion paths, statistical inference is considerably more challenging. This is because the distribution of the diffusion increments over discretely sampled time points is often unknown.

Instead of representing a multivariate diffusion system $\phi_{t}^{*}$ as a stochastic differential equation, one may instead focus on the Kolmogorov forward equation, which dictates the evolution of its probability density function $p(\phi_{t}^{*})$. This is given by:

This is also known as the Fokker-Planck equation. Except where required, the dependence of the drift vector and diffusion tensor on the parameter vector shall be suppressed within the notation. Since a diffusion process is Markovian, the likelihood of a diffusion system sampled at discrete time points $(t_{1},t_{2},...,t_{N})$ is given by:

Asymptotically, for large $N$, the term $p(\phi_{t_{1}}^{*})$ may be ignored whilst the transitional probability distribution $p(\phi_{t_{i}}^{*}|\phi_{t_{i-1}}^{*})$ is the solution to equation (3) at time $t_{i}$ with the boundary condition that $p(\phi_{t}^{*})$ is given by the Dirac delta function centered around $\phi_{t_{i-1}}^{*}$ at time $t_{i-1}$. The likelihood function is central to many inference procedures: it enables us to derive parameter estimates, confidence intervals and to conduct hypothesis tests. Unfortunately, except for a few special cases, equation (3) (and hence also the likelihood) is analytically intractable.

The inability to solve equation (3) is an impediment to statistical inference. This may be circumvented by attempting to match, by choice of parameters, characteristics of the sampled path with characteristics of the diffusion model. For example, one may choose to estimate the instantaneous means and variances using the corresponding sample moments of the differenced data (Gallant and Long, 1997; Ragwitz and Kantz, 2001). Alternatively one may employ Bayesian imputation to augment the observed data so that the diffusion increments are approximately normally distributed (Roberts and Stramer, 2001).

Since likelihood based methods tend to give more precise parameter estimates than method of moments estimators (Hurn et al., 2007), we instead seek an approximation to the transitional probability distribution of the diffusion process. Monte Carlo methods may be used to approximate the likelihood (Kleinhans and Friedrich, 2007; Durham and Gallant, 2002). Alternatively, equation (3) could be solved numerically. Wojtkiewicz and Bergman (2000) discretized the spatial domain and solved the partial differential equation numerically at each point on the lattice. The finite-difference method discretizes the time domain, taking advantage of the fact that over an infinitesimally small time period, the diffusion process is normally distributed (Wehner and Wolfer, 1987). Huang (2012) developed a quasi-maximum likelihood estimator that approximates the first two conditional moments using a Wagner-Platen approximation. The resulting normal distribution can be used to approximate the transitional probability.

Another possibility is to approximate the transitional probability distribution by a closed form analytic function; for example, a Hermite polynomial expansion (Ait-Sahalia, 2002). This method was shown by Ait-Sahalia to be superior to many of the competing methods - both in terms of the accuracy of the transitional distribution approximation as well as the speed of the algorithm. Stramer et al. (2010) created an MCMC procedure based on the Hermite approximation which could allow for measurement errors in the diffusion process.

It must be stressed that the Hermite approximation is only applicable for reducible diffusion processes. A diffusion process $X$ is reducible if there exists a one-to-one transformation $Y=h(X,\theta)$ such that the covariance function of $Y$ is the identity matrix. Though all univariate diffusion processes are reducible, only some multivariate diffusions share this property. Ait-Sahalia (2008) extended the method to irreducible, multivariate diffusions, but not only is the procedure more difficult to implement, there is also a reduction in the accuracy of the closed-form approximation to the transitional density. Furthermore, the Hermite approximation does not in general integrate to one. Indeed, for parameter values far from the maximum likelihood estimates, the normalizer can be very far from one. This often prevents convergence when applying the Hermite approximation within an MCMC setting (Stramer et al., 2010). Consequently, without modification, the resulting MCMC credibility intervals often suffer from considerable undercoverage.

It is proposed that the transitional probability may rather be estimated by a saddlepoint approximation (Daniels, 1954). The saddlepoint approximation is an algebraic expression based on a random variable’s cumulant generation function (CGF). In cases where the first few moments of a random variable is known but the corresponding probability density is difficult to obtain, the saddlepoint approximation to the density can be calculated. The tails of a saddlepoint approximation are more accurate than those of a Edgeworth-expansion (Barndorff-Nielsen and Klüppelberg, 1999). Saddlepoint methods have already been used to approximate the transition densities of diffusions (Ait-Sahalia and Yu, 2006; Preston and Wood, 2012). Preston and Wood apply the saddlepoint approximation to the CGF of a truncated small-time sample-path expansion whilst Ait-Sahalia and Yu (2006) apply the saddlepoint to a truncated expansion, in small-time, of the characteristic function of the transition density.

The saddlepoint approximation will be applied to a truncated expansion of the CGF of $\phi_{t}^{*}$ with respect to the CGF parameters. In contrast to the aforementioned small-time methods, the focus is instead on predicting the evolution of the diffusion system’s cumulants through time, which may be achieved with the cumulant truncation procedure (Whittle, 1957; Gillespie and Renshaw, 2007). The cumulant truncation procedure may also be applied to multivariate diffusion processes (Varughese and Fatti, 2008; Varughese, 2011). These predicted cumulants can be subsequently substituted within a saddlepoint to approximate the transitional densities. After fitting a diffusion model, one may subsequently test for model misspecification of the SDE (Zhang et al., 2012).

In Section 2, the saddlepoint approximation to the transitional probability distribution of the diffusion process is introduced. In Section 3, an MCMC algorithm which uses this approximation for parameter estimation is introduced. Since Ait-Sahalia (2002) demonstrated that the Hermite approximation is superior to a number of alternative approximations, we compare the saddlepoint and the Hermite approximation in Section 4. The two methods are tested on univariate diffusion processes. The Heston model (an irreducible, multivariate diffusion process) is fitted in Section 5 to the S&P 500 and VIX indices. This is followed by a discussion of the various results and some conclusions being drawn.

In this section, a parameter estimation algorithm for a diffusion system is proposed. The algorithm uses the saddlepoint approximation introduced in Section 2 to recreate the transitional probability distribution. This is subsequently used to approximate the likelihood function of the diffusion system. The approximate likelihoods enable us to explore the parameter space of the diffusion system with a modified MCMC procedure. The algorithm is described below.

1. 1.

   Apply the cumulant truncation procedure to the diffusion system thus deriving a system of ordinary differential equations that describe the cumulants’ evolution. (Note these differential equations will depend on the model parameters $\theta$.)

2. 2.

   Choose an arbitrary, starting set of parameter values $\theta_{0}$. Set $\theta_{old}=\theta_{0}$.

3. 3.

   Propose a jump from the old set of parameter values $\theta_{old}$ to a new set of parameter values, $\theta_{new}=\theta_{old}+\Delta\theta$ using a suitably chosen proposal distribution $q(\theta_{new}|\theta_{old})$.

4. 4.

   Calculate the likelihoods $L(\theta_{old})$ and $L(\theta_{new})$. $L(\theta)$ is calculated as follows:

   1. i.

      Set the likelihood, $L(\theta)=1$.

   2. ii.

      Set $i$ to 1.

   3. iii.

      Use the system of ordinary differential equations from Step 1 together with the data values observed at time $t_{i}$, $\phi_{t_{i}}^{*}=(\phi_{1},\phi_{2},...,\phi_{m})$ to predict the values of the cumulants at time $t_{i+1}$.

   4. iv.

      Given the data $\phi_{t_{i}}^{*}$ at time $t_{i}$, the transitional probability distribution at time $t_{i+1}$ can be approximated by a saddlepoint approximation $f_{n}(\phi_{t_{i+1}}^{*}|\phi_{t_{i}}^{*})$ after substituting for the cumulants derived from step iii.

   5. v.

      Set the likelihood $L(\theta)=L(\theta)\times f_{n}(\phi_{t_{i+1}}^{*}|\phi_{t_{i}}^{*})$. This is an implementation of equation (4).

   6. vi.

      Set $i=i+1$ and if $i<N$ (where $N$ is the number of data points) go to step iii.

5. 5.

   Calculate the acceptance ratio $R$. This is given by:

   $$R=\frac{L(\theta_{new})}{L(\theta_{old})}\frac{\pi(\theta_{new})}{\pi(\theta_{old})}\frac{q(\theta_{old}|\theta_{new})}{q(\theta_{new}|\theta_{old})},$$

   where $\pi(\theta)$ is the prior density. We then accept the proposed parameter values $\theta_{new}$ with probability $\min(1,R)$. That is, set $\theta_{new}=\theta_{old}$ with probability $1-\min(1,R)$. Otherwise leave $\theta_{new}$ unchanged.

6. 6.

   Go back to step 3.

After a suitably chosen burn-in period, the MCMC chains sample from the posterior distribution of the diffusion parameters. Hence, in addition to using the medians of the chains as estimates of the diffusion parameters, the $\alpha/2$-th and $(1-\alpha/2)$-th percentile may be used to construct $(100-\alpha)$% credibility intervals.

Since the Hermite approximation (Ait-Sahalia, 2002) is a popular method for estimating the parameters of a diffusion system, we compare its performance against the saddlepoint. In order to judge the relative accuracies of the saddlepoint and Hermite procedures, the true transitional probability distribution is required as a baseline. Unfortunately, the diffusion models for which the transitional distribution is known are few. Amongst, diffusion processes with analytical, non-normal transitional distributions, the Cox-Ingersoll-Ross (CIR) process and Geometric Brownian motion are arguably the most well known. Neither the saddlepoint nor the Hermite approximation in general integrate to one and the approximations are not normalized when performing the comparisons in this section.

Within this section, the CIR process is analyzed didactically: the analysis is meant to illustrate the general derivation of the saddlepoint approximation to a transitional distribution. This is followed by a study of the relative accuracy, for the CIR process, of the Hermite and saddlepoint approximations to the transitional distribution. The section concludes by comparing MCMC implementations for the Hermite and saddlepoint procedures for both the CIR process as well as Geometric Brownian motion.

We study the S&P 500 and it’s relation to the VIX index (a popular measure of implied volatility) over the period December 2009 to November 2010. Figure 4 shows the two indices over this period. The Heston model is fitted to the dataset. It uses proxies for market volatility - such as the VIX index - to account for market movements. The key feature of this model is that the volatility of the asset price movements are assumed to be stochastic.

where $S_{t}$ represents the asset price, $V_{t}$ is the underlying volatility that drives the asset movements and $r$, $\delta$, $\theta$, $\rho$ and $\sigma$ are parameters. The Heston model is a multivariate, irreducible diffusion processes for which many of the competing methods would not be applicable.

We wish to run the saddle MCMC in order to obtain both estimates and credibility intervals for the parameters for the Heston model. To run the saddle MCMC, the diffusion tensor $\sigma_{ij}(S_{t},V_{t};\theta)$ is needed. This is given by:

By applying equations (8) and (10) to the above diffusion process, we obtain the following partial differential equation for $K(\nu_{1},\nu_{2},t)=\ln\mbox{E}[\exp\{\nu_{1}S_{t}+\nu_{2}V_{t}\}]$:

A third order cumulant truncation is performed to yield a system of ordinary differential equations for the cumulants. These may be subsequently used to approximate the likelihood of the Heston model. The modified MCMC (described in section 3) may then be run on the S&P and scaled VIX dataset.

The resulting parameter estimates and 95% credibility intervals are shown in Table 3. One may see that there is strong negative correlation between the two Brownian motions driving $S_{t}$ and $V_{t}$. The parameters $r$ and $\delta$ have wide credibility intervals, suggesting that they are not strongly constrained by the data.

A general statistical framework for estimating the parameters of a diffusion system is proposed. It is assumed that the elements of both the mean vector and the diffusion matrix may be represented by finite order polynomials. Apart from this, the method may be readily applied to multivariate diffusion systems. Many existing parameter estimation methods are only applicable for univariate diffusion processes or, if applicable to multivariate systems, require the diffusion system to be reducible (Ait-Sahalia, 2008; Beskos et al., 2006).

The transitional distribution of a diffusion system is approximated with a saddlepoint. We show how the cumulants of the system may be predicted and subsequently used to approximate the true transitional distribution of the system. The proposed method is compared with the Hermite approximation proposed by Ait-Sahalia (2002). For both the Cox-Ingersoll Ross model and Geometric Brownian motion, the accuracy of the saddlepoint is comparable to the Hermite approximation at the true parameter values. Furthermore, since the Hermite approximation is a small-time expansion, the Hermite approximation rapidly becomes inaccurate as the time interval between successive observations increases.

The saddlepoint approximations may be easily incorporated within an MCMC algorithm. Since, in regions far from the maximum likelihood estimates, the saddlepoint approximation to the transitional density is better behaved than the Hermite approximation, the saddlepoint MCMC chains have better convergence properties. Consequently, the coverage of the credibility intervals for the saddlepoint MCMC is close to those obtained from the true MCMC. In contrast, the credibility intervals for the unmodified Hermite MCMC suffer from undercoverage. This undercoverage is especially severe when the chains are not started at the true parameter values since the chains often become stuck. This is mitigated by the modified algorithm developed by Stramer et al. (2010) which does not require the computation of the posterior’s normalizing constant. However, the modified chains are not only slower than the saddlepoint MCMC, but they take longer to converge, which caused undercoverage for the Geometric Brownian motion example studied in Section 4 since the chains had not converged by the end of the burn-in period.

The comparison between the saddlepoint and Hermite approximations was only performed for univariate diffusion processes and hence caution must be exercised in making claims for multivariate diffusion processes. However, the good coverage observed for the credibility intervals of the saddlepoint MCMC for the non-linear, multivariate example suggests that the procedure still works well for multivariate diffusion processes. Also, for irreducible multivariate processes, the Hermite procedure requires a two-fold approximation, which leads to a further drop in its accuracy (Ait-Sahalia, 2008).

The saddlepoint MCMC is used to fit the Heston model to the S&P 500 and the VIX indices over the period December 2009 to November 2010. Since the Heston model is a multivariate, irreducible diffusion process, the estimation of the model parameters presents a formidable challenge, but the saddlepoint method provides reliable estimates and credibility intervals for the model parameters. As expected, there was significant negative correlation between the two Brownian motions driving the asset prices and the volatility process.

The proposed estimation algorithm has several virtues: it is fast; applicable for a wide range of diffusion systems and gives parameter estimates close to the true values. Furthermore, the corresponding credibility intervals for these estimates have coverage close to their advertised values. This makes the algorithm suitable for estimating the parameters of many diffusion systems.

The MGF of a diffusion system behaves according to the following partial differential equation:

Proof: Consider the multivariate diffusion system $\phi_{t}^{*}$ with instantaneous mean vector $\mu(\phi_{t}^{*},t)$ and diffusion matrix $\sigma(\phi_{t}^{*},t)$. We denote the transition probability by:

Let $\Lambda^{T}=(\upsilon_{1},...,\upsilon_{m})$ be the parameters of the MGF and let $\Lambda^{*}=(\partial/\partial\upsilon_{1},...,\partial/\partial\upsilon_{m})$. The MGF, $M(\Lambda,t)=E\left[\exp[\upsilon_{1}\phi_{1}(t)+...+\upsilon_{m}\phi_{m}(t)]\right]$ obeys the following differential equation (Barbour, 1972):

For small $\Delta x$, we have:

For diffusion processes the higher order terms are negligible. We thus only require the joint and marginal transitional rates of the diffusion system to characterize the evolution of the MGF. We denote the joint transition probabilities for $\phi_{i}$ and $\phi_{j}$ by:

The marginal transition rate is denoted by:

Note that the instantaneous mean $\mu_{i}(\phi_{t}^{*},t)$ and covariances $\sigma_{i}(\phi_{t}^{*},t)$ can be represented in terms of the joint and marginal transitional rates:

By substituting (A.3) within (A.2), we obtain: