Bayesian Approach for Parameter Estimation of Continuous-Time Stochastic Volatility Models using Fourier Transform Methods

TL;DR

This paper introduces a Bayesian method utilizing Fourier transform-based volatility estimation to infer parameters of continuous-time stochastic volatility models, demonstrating high accuracy for Heston models but limited success for more complex models.

Contribution

It develops a novel two-stage Bayesian estimation procedure combining Fourier-based volatility estimation with prior-informed likelihood construction.

Findings

Highly effective for Heston model parameter estimation

Limited success with exponential-Ornstein-Uhlenbeck model

Demonstrates the utility of Fourier methods in Bayesian inference

Abstract

We propose a two stage procedure for the estimation of the parameters of a fairly general, continuous-time stochastic volatility. An important ingredient of the proposed method is the Cuchiero-Teichmann volatility estimator, which is based on Fourier transforms and provides a continuous time estimate of the latent process. This estimate is then used to construct an approximate likelihood for the parameters of interest, whose restrictions are taken into account through prior distributions. The procedure is shown to be highly successful for constructing the posterior distribution of the parameters of a Heston model, while limited success is achieved when applied to the highly parametrized exponential-Ornstein-Uhlenbeck.

Figures

Tables

Table 1. Table 1 : The function ϕ italic-ϕ \phi for the equi-volatility and Exp-OU models.

Model	Specification	$\varphi (v;\theta )$
Equi-Volatility	$f=g$	$\mu -\frac{\kappa \rho m}{\xi }+\frac{\kappa \rho }{\xi }v$
Exp-OU	$f\left(v\right)=e^v$ and $g\equiv 1$	$\mu -\frac{\rho \kappa m}{\xi }{e}^v+\frac{\rho \kappa }{\xi }v{e}^v$

Table 2. Table 2 : Parameters’ Description and Values

Parameter	Description	Value
$T$	Time horizon	1.0
$n$	Number of observations	$2^{19}$
$N$	C-T Frequency	$2^9$
$X_0$	initial log-price	0.0
$\mu$	return rate	0.0

Table 3. Table 3 : Priors and Posteriors distributions

Parameter	Priors
$\rho$	$U[-1,1]$
$\xi$	$U[0,5]$
$\kappa$	$U[0,100]$
$m$	$U[0,1]$
$\mu$	$U[-1,1]$

Equations

$$X_t=x_0+{\int }_0^t{\mu }_{s}ds+{\int }_0^t{\sigma }_{s}d{W}_s,0\le t\le T,$$

$$V_t=v_0+{\int }_0^t{\mu }_s^{V}ds+{\int }_0^t{\sigma }_s^{V}d{W}_s^V,0\le t\le T,$$

$$\displaystyle\left\{\begin{array}[]{l}dX_{t}=\displaystyle\mu dt+f(V_{t})dW_{t},\\ \\ dV_{t}=\kappa(m-V_{t})dt+\xi g(V_{t})dW_{t}^{V},\\ \\ W_{t}^{V}=\rho W_{t}+\sqrt{1-\rho^{2}}W_{t}^{\perp},\end{array}\right.$$

$$X_t-X_s|(V_u{)}_{u\in [0,T]}\sim N\left({\int }_s^t\varphi (V_u;\theta )du+\frac{\rho }{\xi }{\int }_s^t\frac{f(V_u)}{g(V_u)}d{V}_u,(1-{\rho }^2){\int }_s^t{f}^2(V_u)du\right)$$

$$\varphi (v;\theta )$$

$${\int }_s^{t}f(V_u)d{W}_u=\rho {\int }_s^{t}f(V_u)d{W}_u^V+\sqrt{1-{\rho }^2}{\int }_s^{t}f(V_u)d{W}_u^{\perp },$$

$${\int }_s^{t}f(V_u)d{W}_u^{\perp }|(V_u{)}_{u\in [s,t]}\sim N\left(0,{\int }_s^t{f}^2(V_u)du\right).$$

$${\int }_s^{t}f(V_u)d{W}_u^V=\frac{1}{\xi }{\int }_s^t\frac{f(V_u)}{g(V_u)}d{V}_u-\frac{1}{\xi }{\int }_s^t\frac{f(V_u)}{g(V_u)}\kappa (m-V_u)du$$

$$d{X}_t={\mu }_{t}dt+f(V_t)d{W}_t,$$

$${\hat{V}}_t^{n,N}=f^{-1}\left(\sqrt{{\rho }_h^{-1}\left({\widehat{{\rho }_h(V)}}_t^{n,N}\right)}\right),$$

$${\widehat{{\rho }_h(V)}}_t^{n,N}$$

$$G(X,h,k{)}_T^n$$

$$\Delta X_k|(V_u{)}_{u\in [0,T]}\stackrel{ind}{\sim }N\left(M_k(\theta ),(1-{\rho }^2)F_k^2\right)$$

$$M_k(\theta )={\int }_{t_{k-1}}^{t_k}\varphi (V_u;\theta )du+\frac{\rho }{\xi }{\int }_{t_{k-1}}^{t_k}\frac{f(V_u)}{g(V_u)}d{V}_u.$$

$$\displaystyle\ell(\theta;\,\Delta\mathbf{x}\,|\,(V_{u})_{u\in[0,T]})=\big{(}2\pi(1-\rho^{2})\big{)}^{-K/2}\left(\prod_{k=1}^{K}F_{k}^{-1}\right)\exp\left\{-\frac{1}{2(1-\rho^{2})}\sum_{k=1}^{K}\left(\frac{\Delta x_{k}-M_{k}(\theta)}{F_{k}}\right)^{2}\right\}.$$

$$\displaystyle\widehat{\ell}_{n,N}(\theta;\,\Delta\mathbf{x})=\big{(}2\pi(1-\rho^{2})\big{)}^{-K/2}\left(\prod_{k=1}^{K}\widehat{F}^{-1}_{k}\right)\exp\left\{-\frac{1}{2(1-\rho^{2})^{2}}\sum_{k=1}^{K}\left(\frac{\Delta x_{k}-\widehat{M_{k}}(\theta)}{\widehat{F}_{k}}\right)^{2}\right\},$$

$${\hat{F}}_k^2=f^2({\hat{V}}_{t_{k-1}})\Delta t\text{ and }{\hat{M}}_k(\theta )=\varphi ({\hat{V}}_{t_{k-1}};\theta )\Delta t+\frac{\rho }{\xi }\frac{f({\hat{V}}_{t_{k-1}})}{g({\hat{V}}_{t_{k-1}})}\Delta {\hat{V}}_k.$$

$${ϵ}_k=\frac{\Delta X_k-M_k(\theta )}{\sqrt{1-{\rho }^2}{F}_k}\text{mbox}and{\widehat{ϵ}}_k=\frac{\Delta X_k-{\hat{M}}_k(\theta )}{\sqrt{1-{\rho }^2}{\hat{F}}_k}$$

$${\int }_{t_{k-1}}^{t_k}\varphi (V_u;\theta )du=\left(\mu -\frac{\kappa \rho m}{\xi }\right)(t_k-t_{k-1})+\frac{\kappa \rho }{\xi }{\int }_{t_{k-1}}^{t_k}{V}_{u}du$$

Full text

1 Introduction

For a given filtered probability space $(\Omega,\mathcal{A},(\mathcal{A}_{t})_{t\in[0,T]},\mathbb{P})$ with a Brownian motion $W$, a process $X$ is called standard Itô process if it allows the integral representation

[TABLE]

where $\mu$ (drift) and $\sigma>0$ (diffusion coefficient) are adapted measurable processes that satisfy certain conditions to ensure existence of integrals (see, for example, [15]). The diffusion coefficient of an Itô processes, known as volatility in financial models, has been studied intensely due to its immense significance in applications. If the volatility effectively depends on $\omega\in\Omega$ via another standard Itô process $V_{t}$ with

[TABLE]

where $W^{V}$ is a Brownian motion, we say that the process (1.1) is equipped with stochastic volatility. From the theory of Itô processes it follows that $\int_{0}^{t}\sigma^{2}_{s}ds$ (the integrated realized variance) is equal to quadratic variation of $X$ in the interval $[0,t]$, and can be estimated from a discretely observed path of $X$ as the sum of squares of increments on $[0,t]$. Then $\sigma^{2}$ as a function of $t$ (instantaneous variance) can be recovered from the estimate by differentiation. However, the implementation with high frequency samples from real time series of asset prices showed certain drawbacks due to the market noise, as explained and analyzed in Zhou [19] and Zhang et al. [18], see also Aït-Sahalia et al. [1] for a large-scale simulation study of the integrated variance estimator from [18]. Malliavin and Mancino ([12], [13]) offered a non-parametric method for a direct estimation of instantaneous variance, based on Fouried transform. Using similar ideas, Cuchiero and Teichmann [2] proposed a robust estimator, which allows processes with jumps. Both Malliavin-Mancino and Cuchiero-Teichmann estimators are applicable in a multidimensional setup. There are several other instantaneous variance estimators (see [2] and [13] for references).

Diffusion processes usually involve some parameters, which, in practice, would have to be estimated. There is a voluminous literature on Bayesian inference in this context (for example, [4], [5], [9], [16]). In this paper we consider a general family of continuous-time stochastic volatility models with five parameters. With particular specifications, this family covers several well known models, such as the Heston and exponential-Ornstein-Uhlenbeck (exp-OU) models. We propose a two-stage method for parameter estimation. At the first stage, we recover the realized volatility process using Cucheiro-Teichmann estimator, and at the second stage, conditionally on the recovered volatility, we estimate the parameters using Bayesian technique.

2 A general stochastic volatility model: properties and volatility estimation

In this section we present the model and its properties, together with the one-dimensional version of the Cucheiro-Teichmann procedure [2]. The Bayesian parameter estimation will be discussed in the Section 3.

2.1 The model and its properties

Let $W$ and $W^{\perp}$ be independent standard Brownian motions on the fixed probability space $(\Omega,\mathcal{A},\mathbb{P})$. Consider the following stochastic volatility model in the differential form, as a system of stochastic differential equations (SDEs):

[TABLE]

where $\rho\in[-1,1]$. This parameter reflects a correlation between $W$ and $W^{V}$. Other parameters of the model under consideration and their ranges are $\mu\in\mathbb{R}$, $m\in\mathbb{R}$, $\kappa>0$, $\xi>0$. The usual names in the financial literature are exhibited in Table 2.

For example, $f(v)=g(v)=\sqrt{v}$ gives us the Heston model, see [8], and $f(v)=e^{v},g(v)=1$, the exponential-Ornstein Uhlenbeck (exp-OU) model, see [14]. Different models (that consider different drifts for the volatility) could be considered with little modification of the technique presented below. A class of models with $f=Cg$, where $C>0$ is a specified known constant, is frequently used in applications. We call these models equi-volatility models; the Heston model is one example of such model. The Inverse Gamma model of [10] is also a member of this class.

The additional necessary assumptions on the model are stated as follows.

Assumption 2.1**.**

$(i)$ $f$ and $g$ are positive functions on the support of the volatility process $V$, and the function $f$ is strictly monotone; $(ii)$ The SDE (2.6) has a unique strong solution, with a positive volatility process $V$.

The following theorem states that the increments $X_{t}-X_{s}$, given the volatility path, are normally distributed with mean and variance that can be explicitly computed. This result is a paramount in our Bayesian estimation procedure described in Section 3.

Theorem 2.2**.**

Under Assumption 2.1, given the volatility path $(V_{u})_{u\in[0,T]}$, the increments $X_{t}-X_{s}$, $0\leq s<t\leq T$, are independent over disjoint intervals and normally distributed with mean and variance as below:

[TABLE]

where $\theta=(\mu,\kappa,m,\rho,\xi)$ and

[TABLE]

Proof.

From (2.6) we find that

[TABLE]

where $W^{\perp}$ is a Brownian motion independent of $W^{V}$. Then it follows that

[TABLE]

Moreover, from the SDE describing the dynamics of the process $V$, we find that

[TABLE]

and this yields the result. ∎

In Table 1, we specify the function $\phi$ of two particular models.

2.2 Pathwise covariance estimation of Cuchiero and Teichmann

The estimation procedure for the parameters of the SV model discussed in the previous subsection is built on the following estimation method of the hidden volatility process, which has been proposed and shown to be consistent in [2]. We will now describe this method. More generally, we assume that $X$ follows the dynamics

[TABLE]

where $W$ is an one-dimensional Brownian motion, $\mu$ is a locally bounded process and $V$ is a continuous stochastic process. Fix a time horizon $T>0$ and define $s_{m}^{n}=m/n$, for $m=0,\ldots,\lfloor{nT}\rfloor$. It is assumed that we observe the process $X$ at times $\{s_{m}^{n}\}_{m=0}^{\lfloor{nT}\rfloor}$.

Given a continuous function $h$ with at most polynomial growth, the Cuchiero-Teichmann estimator of the instantaneous variance is given by

[TABLE]

where

[TABLE]

and the function $x\mapsto\rho_{h}(x)$ is defined as $\rho_{h}(x)=\mathbb{E}[h(Z)]$, for $Z\sim N(0,x)$.

Example 2.3**.**

For instance, we may choose $h(x)=\cos(x)$, which gives us $\rho_{h}(x)=e^{-\frac{1}{2}x}$ and $\rho_{h}^{-1}(x)=-2\log x$.

The consistency of the estimator is guaranteed if $\displaystyle\lim_{n,N\to+\infty}n/N^{\gamma}>0$, for some $\gamma>1$. One might take $N=\lfloor{\sqrt{n}}\rfloor$.

Even though the volatility estimator from (2.9) is defined for continuous time $t$, one could evaluate the process $\widehat{V}^{n,N}$ at times $t_{k}=kT/(2N+1)$, $k=0,\ldots,2N+1$. We choose $n$ and $N$ in such a way that every $t_{k}$ is one of the $s_{m}^{n}$, i.e. we may consider both $X$ and $\widehat{V}^{n,N}$ at times $t_{k}$.

3 Bayesian estimation procedures for the parameters for stochastic volatility models

In this section, we will describe the estimation procedure based on Theorem 2.2 and the volatility estimation of Cuchiero-Teichmann presented in Section 2.2.

For a sequence of equally spaced time steps $(t_{k})_{k=1,\ldots,K}$, with $\Delta t=t_{k}-t_{k-1}$, we define the increments of the log-price process as $\Delta X_{k}=X_{t_{k}}-X_{t_{k-1}}$. Theorem 2.2 states that, conditional on the entire path of true volatility process $V$, the increments are independent and normally distributed:

[TABLE]

with $F_{k}^{2}=\int_{t_{k-1}}^{t_{k}}f^{2}(V_{u})du$ and

[TABLE]

Then, under model (3.1), the likelihood for the parameter $\theta$, conditional on the volatility process, is given by

[TABLE]

As the process $V$ is non-observable, in order to be able to perform parametric inference for $\theta$, we use an approximate likelihood approach (see, e.g., [3] for several simulation-based approaches in the context of Bayesian inference). For fixed $n$ and $N$, let $\widehat{\ell}_{n,N}$ denote the approximated likelihood, defined as

[TABLE]

where the process $V$ is replaced by its Cuchiero-Teichmann estimate, $\widehat{V}$. In (3.3) all the integrals are approximated by quadrature, providing the following definitions:

[TABLE]

Note that, for the sake of notational simplicity, we are dropping the superscript $n,N$ from the volatility estimate and assuming that $K=2N+1$.

Remark 3.1**.**

For fixed $(t_{k})_{k=1,\ldots,K}$, under some additional regularity assumptions, one could show the convergence of the approximated likelihood $\widehat{\ell}_{n,N}$ to $\ell$, as $n$ and $N$ go to infinity satisfying $\displaystyle\lim_{n,N\to+\infty}n/N^{\gamma}>0$, for some $\gamma>1$. Although the rigorous verification of this convergence is outside the scope of this letter, in Figure 2 we motivate this result by comparing the quantiles of two different normal distributions:

[TABLE]

Note that on both distributions the parameters are set to their true values and on the first model we condition on the true volatility path.

3.1 Parameter identifiability

Firstly, the variance of the increments of $X$ depends only on the unknown parameter $\rho$ through its square. This means that absolute value of $\rho$ is identifiable. See A for a more throughout discussion on the identifiability issue for parameter estimation. Secondly, the presence of a Brownian motion in $V$ makes the increments of $F(V)$ and the integrals of $\phi(V;\theta)$ fundamentally different, i.e. the increments of $F(V)$ have the rough behavior of the Brownian motion and the integrals of $\phi(V;\theta)$ are of bounded variation, which gives us a strong indication that $\beta$ is identifiable.

Therefore, since $\xi>0$, we have $\mbox{sign}(\rho)=\mbox{sign}(\beta)$, which implies $\rho=|\rho|\mbox{sign}(\beta)$ and $\xi=\rho/\beta$ are identifiable in the estimation of the general statistical model (3.1). In order to study the identifiability of the other parameters $\mu,\kappa$ and $m$, one would need to specify the volatility functions $f$ and $g$.

In the equi-volatility models, using the specification of $\phi$ given in Table 1, the first integral in (3.2) involves a constant and an integral of $V$:

[TABLE]

Since these two terms are different functions of the data, we may conclude that we can identify the coefficients of the time increment and the integral of $V$, which implies that $\kappa$ is identifiable, but $\mu$ and $m$ are not. In the exp-OU model, following a similar reasoning, we conclude that $\mu$, $m$ and $\kappa$ are identifiable.

Remark 3.2** (A drawback of frequentist inference).**

As the model in (3.3) is a simple linear regression it is clear that a pure likelihood-based estimation procedure could possibly estimate the variance term $1-\rho^{2}$ by a negative value. In order to avoid inestimability, the Bayesian procedure proposed in the next section assigns the uniform prior on the interval $[-1,1]$ to $\rho$.

3.2 The estimation procedure and implementation on selected models

In order to construct the posterior distribution of the parameters of interest, the (approximated) likelihood in (3.3) is combined with appropriated prior distributions. Samples from the posterior distribution are generated using the following procedure:

from the observed data $X_{s_{1}},\ldots,X_{s_{n}}$, estimate $V_{t_{1}},\ldots,V_{t_{N}}$ using the Cuchiero-Teichmann procedure described in Section 2.2; 2. 2.

assign independent prior distributions to $\xi$ and $\rho$ respecting the restrictions that $\xi>0$ and $\rho\in[-1,1]$;

• –

in our procedure, we choose uninformative priors on (sensible) finite intervals, see Table 3; 3. 3.

using $(X_{t_{1}},\widehat{V}_{t_{1}}),\ldots,(X_{t_{N}},\widehat{V}_{t_{N}})$ and Equation (3.3), generate samples from the posterior distributions of $(\xi,\rho)$;

1. (a)

Under the equi-volatility models, we may also generate samples from the posterior distribution of $\kappa$, but we cannot separate the effects of $\mu$ and $m$, see Table (1). 2. (b)

Under the Exp-OU model, we are able to generate samples from the posterior distribution of all the other parameters: $\kappa,m$ and $\mu$.

Step 3 of the above procedure is performed with the aid of a Hamiltonian Monte Carlo algorithm, implemented through R-Stan [17].

4 Numerical exercise

The numerical experiment will consider simulated data from the Heston model (an example of equi-volatility model) and the exp-OU model (where $f$ and $g$ are different). The function $h$ for both models is defined as in the Example 2.3. The parameters (of the model and the numerical procedure) are described in Table 2 and the prior distributions in Table 3.

In Figure 1 the true volatility process (dark blue) and its Cuchiero-Teichmann estimate (light blue) are presented, both for the (a) Heston and the (b) Exp-OU models. It can be seen that the estimate is able to closely follow the true unobserved paths for both models, with the estimate for the Heston model been (at least visually) more precise. As the estimate is offline, in the sense that it its computed for a batch of data and needs to be recomputed when new observation arrive, errors at the beginning/end of the estimation window should not be notoriously different. It should also be noticed that the estimates are reliable at both low and high (absolute) volatility regimes, see Figure 1. One important aspect, though, is the apparent lower volatility of the estimated paths, i.e., the reconstructed functions appear to be smoother than the original ones.

As hinted by Figure 1, Figure 3 shows that the parameters of the Heston model are, indeed, better estimated when compared to the exp-OU model. For both models we present all two dimensional posterior distributions in the lower triangular plots. At these plots we also include (in blue) the real parameter values used in the data simulation. In the main diagonal plots we present the histogram of the identifiable parameters, where the blue solid line represents the real value and the dashed line its marginal posterior mean. Both for the Heston and the exp-OU models the correlation parameter $\rho$ is estimated remarkably well, with high posterior probability been assigned to its correct sign (negative on both cases). While the two remaining parameters of the Heston model are also very well estimated, the same does not hold true for the other exp-OU parameters, as most of them have uniform marginal posterior distributions.

It should be stressed that these results are consistently observed for these models, independent of the particular sample path generated, which leads us to believe that the proposed method can be a competitive alternative to the state-of-the art algorithms used for inference in continuous time stochastic volatility models.

Acknowledgemens. We express our gratitude to Josef Teichmann for sharing with us the code for calculation of the estimator that we used in this paper. Milan Merkle acknowledges the support by grants III 44006 and 174024 from Ministry of Education, Science and Technological Development of Republic of Serbia. Yuri F. Saporito acknowledges the support by grant 210.168/2017 from Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro.

Appendix A Identifiability

Given a stochastic model with a parameter $\theta\in\Theta$, let $\hat{\theta}$ be an estimator for $\theta$ based on a sample of size $n$. Let $\hat{\theta}\;|\;\theta_{0}$ denote the random element $\hat{\theta}$ under condition that the true value of the parameter $\theta$ equals $\theta_{0}$. It is said that $\theta$ is non-identifiable if for any $K\in\mathbb{N}$ there is a sample of size $n\geq K$ and two distinct values $\theta_{1},\theta_{2}\in\Theta$ such that the random elements $\hat{\theta}\;|\;\theta_{1}$ and $\hat{\theta}\;|\;\theta_{2}$ have the same probability distribution. This can be one among many possible ways to define the phenomenon that has been discussed through statistical literature since long ago. There is a considerable number of papers on the topic, especially in the econometric context, see [11] or [7] in the framework of Bayesian methods, and many others.

The parameter $\theta$ is said to be identifiable if it is not non-identifiable. Note that identifiability property is relative to the given estimator. As shown in [6], the identifiability is a weaker property than consistency: if $\hat{\theta}$ is consistent estimator of $\theta$, then $\theta$ is identifiable; however $\theta$ can be identifiable even if there does not exist a consistent estimator.