Realized volatility and parametric estimation of Heston SDEs

TL;DR

This paper develops and analyzes observable moment-based estimators for Heston stochastic differential equations, providing explicit convergence rates and practical guidelines for parameter estimation from stock price data.

Contribution

The paper introduces new observable moment estimators for Heston SDE parameters, with explicit convergence bounds and practical recommendations for data sampling.

Findings

Explicit $L^q$ convergence rates for realized volatilities.

Guidelines for selecting sampling frequency and window size.

Numerical simulations demonstrating estimator performance.

Abstract

We present a detailed analysis of \emph{observable} moments based parameter estimators for the Heston SDEs jointly driving the rate of returns $R_t$ and the squared volatilities $V_t$. Since volatilities are not directly observable, our parameter estimators are constructed from empirical moments of realized volatilities $Y_t$, which are of course observable. Realized volatilities are computed over sliding windows of size $\varepsilon$, partitioned into $J(\varepsilon)$ intervals. We establish criteria for the joint selection of $J(\varepsilon)$ and of the sub-sampling frequency of return rates data. We obtain explicit bounds for the $L^q$ speed of convergence of realized volatilities to true volatilities as $\varepsilon \to 0$. In turn, these bounds provide also $L^q$ speeds of convergence of our observable estimators for the parameters of the Heston volatility SDE. Our theoretical analysis is supplemented by extensive numerical simulations of joint Heston SDEs to investigate the actual performances of our moments based parameter estimators. Our results provide practical guidelines for adequately fitting Heston SDEs parameters to observed stock prices series.

Tables

Table 1. Table 1: Values of estimated E ^ 2 superscript ^ 𝐸 2 \hat{E}^{2} and E ^ 4 superscript ^ 𝐸 4 \hat{E}^{4} errors as defined in ( 75 ) in numerical simulations with J ​ ( ε ) = 10 𝐽 𝜀 10 J(\varepsilon)=10 , 40 40 40 , 1 / ε 1 𝜀 1/\varepsilon , and 1 / ε 2 1 superscript 𝜀 2 1/\varepsilon^{2} . 95% confidence intervals are indicated with “ ± plus-or-minus \pm ” numbers.

$\epsilon =$	0.005	0.01	0.02	0.05	0.1
${\hat{E}}^2$, $J=10$	$1.85\pm 0.13$	$1.85\pm 0.14$	$1.86\pm 0.14$	$1.88\pm 0.14$	$1.94\pm 0.14$
${\hat{E}}^4$, $J=10$	$2.98\pm 0.58$	$3.00\pm 0.6$	$3.01\pm 0.63$	$2.99\pm 0.52$	$3.09\pm 0.58$
${\hat{E}}^2$, $J=40$	$0.95\pm 0.07$	$0.96\pm 0.07$	$0.98\pm 0.07$	$1.05\pm 0.07$	$1.15\pm 0.08$
${\hat{E}}^4$, $J=40$	$1.52\pm 0.24$	$1.53\pm 0.24$	$1.56\pm 0.27$	$1.63\pm 0.22$	$1.76\pm 0.24$
${\hat{E}}^2$, $J=1/\epsilon$	$0.45\pm 0.03$	$0.64\pm 0.05$	$0.90\pm 0.07$	$1.39\pm 0.09$	$1.94\pm 0.15$
${\hat{E}}^4$, $J=1/\epsilon$	$0.71\pm 0.1$	$1.01\pm 0.15$	$1.41\pm 0.24$	$2.20\pm 0.38$	$3.09\pm 0.58$
${\hat{E}}^2$, $J=1/{\epsilon }^2$	$0.16\pm 0.008$	$0.23\pm 0.013$	$0.34\pm 0.019$	$0.57\pm 0.033$	$0.90\pm 0.056$
${\hat{E}}^4$, $J=1/{\epsilon }^2$	$0.23\pm 0.017$	$0.32\pm 0.028$	$0.48\pm 0.041$	$0.83\pm 0.074$	$1.34\pm 0.16$

Equations

$$d{R}_t=\mu dt+\sqrt{V_t}d{Z}_t,$$

$$d{R}_t$$

$$d{V}_t$$

$$\frac{\kappa \theta }{{\gamma }^2}>\frac{1}{2}.$$

$$\mathbb{E}[U|V_0=y]\equiv {\mathbb{E}}_y[U]\text{for any random variable U}.$$

$$\Vert V_t{\Vert }_q\le C\text{and}\Vert V_t-V_s{\Vert }_q\le C\sqrt{t-s},$$

$$\Vert R_t{\Vert }_q\le C\text{and}\Vert R_t-R_s{\Vert }_q\le C\sqrt{t-s}.$$

$$\mathbb{E}\left[{\left({\int }_s^t{X}_{u}d{W}_u\right)}^q\right]\le H_q\mathbb{E}\left[{\left({\int }_s^t{X}_u^{2}du\right)}^{q/2}\right],$$

$${\left|\left|{\int }_s^t{X}_{u}d{W}_u\right|\right|}_q\le h_q{\left({\left|\left|{\int }_s^t{X}_u^{2}du\right|\right|}_{q/2}\right)}^{1/2}.$$

$${\left|\left|{\int }_s^t\sqrt{V_u}d{B}_u\right|\right|}_q\le h_q{\left({\left|\left|{\int }_s^t{V}_{u}du\right|\right|}_{q/2}\right)}^{1/2}\le h_q{\left({\int }_s^t\Vert V_u{\Vert }_{q/2}du\right)}^{1/2},$$

$${\left|\left|{\int }_s^t\sqrt{V_u}d{Z}_u\right|\right|}_q\le h_q{\left({\left|\left|{\int }_s^t{V}_{u}du\right|\right|}_{q/2}\right)}^{1/2}\le h_q{\left({\int }_s^t\Vert V_u{\Vert }_{q/2}du\right)}^{1/2}.$$

$${\left|\left|{\int }_s^t\sqrt{V_u}d{B}_u\right|\right|}_q\le k_q(t-s{)}^{1/2},$$

$${\left|\left|{\int }_s^t\sqrt{V_u}d{Z}_u\right|\right|}_q\le k_q(t-s{)}^{1/2}.$$

$$V_t-V_s=(t-s)\kappa \theta -\kappa {\int }_s^t{V}_{u}du+\gamma {\int }_s^t\sqrt{V_u}d{B}_u,$$

$$R_t-R_s=(t-s)\mu +{\int }_s^t\sqrt{V_u}d{Z}_u.$$

$$\Vert V_t-V_s{\Vert }_q\le {\alpha }_q(t-s{)}^{1/2}\text{ and }\Vert R_t-R_s{\Vert }_q\le {\delta }_q(t-s{)}^{1/2}$$

$${\delta }_q=\sqrt{T}\mu +k_q\text{ and }{\alpha }_q=\sqrt{T}\kappa (\theta +c_q)+\gamma k_q.$$

$$t_n=t-\epsilon +n\epsilon /J(\epsilon )\text{ for }n=0,\dots ,J(\epsilon ).$$

$$Y_t^{\epsilon }=\frac{1}{\epsilon }\sum _{n=1}^{J(\epsilon )}(R_{t_n}-R_{t_{n-1}}{)}^2.$$

$$\underset{\epsilon \to 0}{\lim }J(\epsilon )=+\infty .$$

$$\Vert Y_t^{\epsilon }-V_t{\Vert }_q\le c\left(\frac{1}{J^{1/q}(\epsilon )}+{\epsilon }^{1/2}\right).$$

$$\Vert Y_t^{\epsilon }-V_t{\Vert }_q\le c\left(1+1/\sqrt{a}\right){\epsilon }^{1/2}.$$

$$\big{|}{\mathbb{E}}[W_{1}W_{2}\ldots W_{q}]\big{|}\leq\|W_{1}\|_{q}\|W_{2}\|_{q}\ldots\|W_{q}\|_{q}.$$

$$\Vert V_t{\Vert }_q+\Vert R_t{\Vert }_q\le C_q\text{ and }\Vert V_t-V_s{\Vert }_q+\Vert R_t-R_s{\Vert }_q\le C_q(t-s{)}^{1/2}.$$

$$(R_t-R_s{)}^2={\int }_s^t{V}_{u}du+2{\int }_s^t(R_u-R_s)\sqrt{V_u}d{Z}_u.$$

$$D(s,t)=(R_t-R_s{)}^2-{\int }_s^t{V}_{u}du$$

$$\mathbb{E}(D(s,t)|{\mathcal{G}}_s)=0.$$

$$\Vert (R_t-R_s{)}^2{\Vert }_q={\left(\Vert R_t-R_s{\Vert }_{2q}\right)}^2\le C_{2q}(t-s)$$

$$\Vert D(s,t){\Vert }_q\le \Vert (R_t-R_s{)}^2{\Vert }_q+{\int }_s^t\Vert V_u{\Vert }_{q}du\le b_q(t-s)$$

$$Y_t^{\epsilon }=\sum _{n=1}^J(R_{t_n}-R_{t_{n-1}}{)}^2.$$

$$Y_t^{\epsilon }-V_t=H(t,\epsilon )+K(t,\epsilon ).$$

Full text

We obtain explicit bounds for the $L^{q}$ speed of convergence of realized volatilities to true volatilities as $\varepsilon\to 0$. In turn, these bounds provide also $L^{q}$ speeds of convergence of our observable estimators for the parameters of the Heston volatility SDE.

Our theoretical analysis is supplemented by extensive numerical simulations of joint Heston SDEs to investigate the actual performances of our moments based parameter estimators. Our results provide practical guidelines for adequately fitting Heston SDEs parameters to observed stock prices series.

Keywords: Heston model, parameter estimation, realized volatility, indirect observability

1 Introduction

Parametric estimation of stochastic differential equations (SDEs) has been an active research area for several decades. The majority of published results focus on Direct Observability situations, where the observable data $X_{t}$ are assumed to be generated by the SDEs themselves. But in many practical situations, the SDEs driving an unobservable process $X_{t}$ are parametrized by a vector $\boldsymbol{\theta}$ which needs to be estimated from observable data $Y_{t}^{\varepsilon}$ which are only known to converge to $X_{t}$ as $\varepsilon\to 0$. We refer to these situations as Indirect Observability contexts. A crucial point is then to assess estimation errors due to the use of approximate data (see, for example, [24, 36, 35, 39, 4, 12, 33]). In our papers [7, 6, 8, 11], we analyzed asymptotic consistency of parameter estimation under indirect observability in multiple contexts. In particular, in [11] we proved the asymptotic accuracy of parameter estimators based on empirical moments of indirect approximate observations, for a wide class of unobservable stationary non-Gaussian processes $X_{t}$ with “fast” mixing properties. Here we extend and deepen results from [11] to parameter estimation for the well known Heston SDEs [31] driving jointly the rate of returns $R_{t}$ of an arbitrary asset and its squared volatility $V_{t}$. Since the volatilities are not directly observable, classical observable approximations of $V_{t}$ are provided by realized volatilities $Y_{t}^{\varepsilon}$ computed on averaging time windows $(t-\varepsilon,t)$. Such volatility approximations have been studied for instance in [32, 30, 13, 20, 21, 3, 38, 5].

In this paper focused on feasible parameter estimation for the Heston volatility SDEs, we construct observable parameter estimators from the first and second order empirical moments of the realized volatility process $Y_{t}^{\varepsilon}$, and analyze their $L^{2}$-consistency as $\varepsilon\to 0$. In order to do this, we extend the $L^{2}$ convergence of the realized volatility to $L^{q}$ (where $q$ is odd). Maximum Likelihood estimation and estimates on $L^{q}$ norms for square-root diffusions were also addressed for example in [16, 17]. While the Maximum Likelihood Estimators (MLEs) have been used in many context to estimate parameters of stochastic differential equations, including the Heston model (e.g. [10, 29]), MLEs can be quite sensitive to model errors. We expect moment estimators of low order to be much more robust with respect to small perturbations of the underlying model fitted to observable and possibly noisy data. The empirical moments of $Y_{t}^{\varepsilon}$ rely on explicit sub-sampling schemes which specify key computational parameters (e.g. number of points in the window $(t-\varepsilon,t)$, observational time-step, and total number of observations) as functions of the window size $\varepsilon$. In particular, optimal sub-sampling schemes involve explicit expressions for selecting the observational time-step $\Delta(\varepsilon)$ and number of observations $N(\varepsilon)$ for the realized volatility. We demonstrate that the optimal speed of convergence for moment estimators computed under indirect observability is $O(\varepsilon^{1/2})$.

When realized volatilities are computed over sliding windows of small duration $\varepsilon$, our target is to determine nearly optimal stockprice sub-sampling rates enabling good control of estimation errors for the parameters of the Heston SDE driving the (unobservable) volatilities. In contrast with other results on parametric estimation of the volatility Heston SDE under indirect observability (e.g. [30]), our results address estimation of both drift and diffusion parameters in the volatility SDE. Parametric estimation of the diffusion term for generic SDEs is a delicate task; several methods have been introduced [32, 18, 29, 22], but comparing performance and robustness for various estimators still remains an active research area.

Application of the general theory developed in [11] requires a substantial analytical investigation of the Heston model. For realized volatilities $Y_{t}^{\varepsilon}$ computed on sliding windows of length $\varepsilon\to 0$, we give concrete estimates for the $L^{q}$ convergence speed of $Y_{t}^{\varepsilon}$ to true volatility $V_{t}$ and we derive explicit nearly optimal sub-sampling schemes of $Y_{t}^{\varepsilon}$ for consistent estimation of empirical moments. We compute theoretical convergence speeds for our observable estimators of the Heston SDE parameters as $\varepsilon\to 0$ and we compare them to numerically evaluated convergence speeds. To this end, we perform numerical investigation of the Heston model, through extensive simulations with $\varepsilon\to 0$. Our simulations refine and confirm our theoretical convergence rates for the realized volatilities as well as for our estimators of the Heston parameters. We thus validate asymptotically optimal ranges for the number of data points used to compute each realized volatility. Our numerical results indicate that for small but realistic values of $\varepsilon$, nearly optimal convergence rates of our observable parameter estimators can still be achieved with data subsampling less frequent than the theoretically prescribed rates.

We introduce the Heston model and address the $L^{q}$-Hölder continuity for squared volatilities in section 2. We introduce the realized volatility and the notion of indirect observability in section 2.4. We prove the $L^{q}$ convergence of realized volatilities in section 3. Some analytic properties of the Heston model are discussed in sections 5, 6, and moment-based parameter estimators for the volatility process are introduced in section 7. Theorem 3 in section 7 is one of our main analytical results, since it computes $L^{2}$ convergence rates for subsampled empirical moments of realized volatilities and hence yields convergence rates for our observable parameter estimators based on realized volatilities. In section 8 we discuss pragmatic discretization rates for the averaging windows $(t-\varepsilon,t)$. In sections 9 and 10 we perform an extensive numerical investigation of the Heston model, including numerical computation of the $L^{2}$ and $L^{4}$ speeds of convergence for the realized volatility process, speed of convergence of parameter estimators, and empirical covariance estimators. Conclusions are presented in section 11.

2 The Heston stochastic volatility models

2.1 Generic stochastic volatility models

In the well known Barnsdorff-Nielsen paper [13], generic stochastic volatility models consider asset price processes $A_{t}$ such that the rate of return process $dR_{t}=dA_{t}/A_{t}$ is driven by an SDE of the form

[TABLE]

where $\mu$ is a constant, $Z_{t}$ is a standard one-dimensional Brownian motion, and the square integrable continuous process $V_{t}>0$ is called the spot variance or squared volatility of the return rate. In this paper we will focus on the classical Heston joint SDEs which are widely used examples of stochastic volatility models.

2.2 The Heston joint SDEs

Recall that in the Heston model [31] for the stochastic dynamics of asset price $A_{t}$ and squared volatility $V_{t}$, two coupled SDEs jointly drive $V_{t}$ and the rate of return $dR_{t}=dA_{t}/A_{t}$, namely

[TABLE]

Here $Z_{t}$ and $B_{t}$ are standard one dimensional Brownian motions with constant instantaneous correlation ${\mathbb{E}}(dZ_{t}dB_{t})=\rho dt$ where $-1<\rho<1$.

The autonomous volatility SDE (3) is parametrized by 3 parameters, the “long run mean” $\theta>0$ of $V_{t}$, the “reversion rate” $\kappa>0$, and $\gamma>0$. To ensure that $V_{t}$ remains almost surely positive for all $t$ provided $V_{0}$ is almost surely positive, the parameter vector $\boldsymbol{\theta}=\left[\kappa,\theta,\gamma\right]$ must verify the classical Feller condition [26]

[TABLE]

In this paper we assume that parameters in the equation for volatility satisfy the Feller condition above. The first Heston SDE (2) is parametrized by the constant “mean return rate” $\mu>0$ of the asset price, and by the correlation coefficient $\rho$. To be specific, we define $Z_{t}=\rho B_{t}+(1-\rho^{2})^{1/2}\beta_{t}$, where $B_{t}$ and $\beta_{t}$ are two independent standard Brownian motions. Let ${\cal F}_{t}$ be the increasing filtration generated by $B_{t}$. Let ${\cal G}_{t}$ be the increasing filtration generated by the pair $(Z_{t},B_{t})$, or equivalently by $(\beta_{t},B_{t})$. Then the joint SDEs (2), (3) have a unique solution $V_{t},R_{t}$ starting at any fixed $V_{0}>0$ and $R_{0}$. Moreover, $V_{t}$ is nonanticipating with respect to both ${\cal F}_{t}$, ${\cal G}_{t}$, and $R_{t}$ is non-anticipating with respect to ${\cal G}_{t}$. Both $V_{t}$ and $R_{t}$ are continuous in $t$. We will use the shortcut notation

[TABLE]

2.3 $L^{q}$-Hölder continuity for squared volatilities

We now prove that the solutions $R_{t},V_{t}$ of Heston SDEs are Hölder continuous in $L^{q}$, with Hölder exponent $1/2$.

Proposition 1**.**

Let $R_{t},V_{t}$ be the returns rate and squared volatility jointly driven by the Heston SDEs (2), (3), starting at any fixed $V_{0}=y>0$ and $R_{0}=r$. Fix any time $T>0$. Then $(T,y,r)$ and the Heston SDEs parameters determine for each $q\geq 2$ a constant $C=C(q)$ such that, for all $0\leq s\leq t\leq T$,

[TABLE]

Moreover, equations (5) and (6) still hold when $V_{t}$ is the unique stationary process driven by the Heston volatility SDE, and $R_{0}$ is fixed.

Proof.

Consider an increasing filtration $\mathcal{W}_{t}$ and let $W_{t}$ be a progressively measurable standard Brownian motion with respect to $\mathcal{W}_{t}$. Let $\mathcal{X}$ be the set of all continuous processes $\{X_{t}\}$, which are progressively measurable with respect to $\mathcal{W}_{t}$ and non anticipative with respect to $W_{t}$. A classical martingale inequality (see equation (3.1) in [19]) states that for each positive $q$, there is a universal constant $H_{q}$ such that for all $(X_{t})\in\mathcal{X}$ and all $0\leq s\leq t$,

[TABLE]

which, denoting $h_{q}=H_{q}^{1/q}$, is clearly equivalent to

[TABLE]

Fix $q\geq 2$, $T>0$. Let $R_{t},V_{t}$ be the solutions of the Heston SDEs starting at any fixed $V_{0}>0$ and $R_{0}$. Denote ${\cal F}_{t}$ and ${\cal G}_{t}$ the filtrations respectively generated by $B_{t}$ and by the pair $(Z_{t},B_{t})$. Then the solution $V_{t}$ of the volatility SDE is ${\cal F}_{t}$ measurable and non-anticipating with respect to $B_{t}$. The solution $R_{t}$ of the returns rate SDE is ${\cal G}_{t}$ measurable and non-anticipating with respect to $Z_{t}$.

Next, apply (8) first with $\mathcal{W}_{t}={\cal F}_{t}$, the $W_{u}=B_{u}$, $X_{u}=\sqrt{V_{u}}$ , and then a second time with $\mathcal{W}_{t}={\cal G}_{t}$, $W_{u}=Z_{u}$, $X_{u}=\sqrt{V_{u}}$. This yields, for all $0\leq s\leq t$, all $q\geq 2$,

[TABLE]

[TABLE]

By Proposition 2 which is proved further on, there is a constant $c_{q}$ such that $\|V_{t}\|_{q}\leq c_{q}$ for all $t$. Hence (9) and (10) imply, with $k_{q}=h_{q}\sqrt{c_{q/2}}$,

[TABLE]

[TABLE]

Integrating the Heston SDEs (2), (3) we obtain

[TABLE]

[TABLE]

Due to (11), (12), this implies for $0\leq s\leq t\leq T$,

[TABLE]

with

[TABLE]

This proves (5) for $V_{0}>0$ and $R_{0}$ fixed. A fully similar proof holds when the volatilities $V_{t}$ are stationary, with only $R_{0}$ fixed.

2.4 Realized Volatilities and actual volatilities

Daily or Intraday market data provide observed asset prices $A_{t}$ at discretized times t, and hence discretized versions of the rate of return $R_{t}$, but the spot variances $V_{t}$ cannot be directly observed or precisely derived from market data and hence are unobservable. However observable approximations of $V_{t}$ are provided for each small $\varepsilon>0$ by the Realized Volatilities $Y_{t}^{\varepsilon}$ computed as follows from the discretized rates of returns $R_{s}$.

Partition each sliding time window $W_{t}^{\varepsilon}=\left[t-\varepsilon,t\right]$ into $J(\varepsilon)$ equal intervals, we define $J(\varepsilon)+1$ time-instants

[TABLE]

The realized volatilities $Y_{t}^{\varepsilon}$ are then computed by the formula

[TABLE]

We will always assume that the partition size $J(\varepsilon)$ verifies

[TABLE]

3 $L^{q}$-approximation of Heston Volatilities by Realized Volatilities

As shown in [13], when $\varepsilon\to 0$, the realized volatilities $Y_{t}^{\varepsilon}$ must converge to $V_{t}$ in $L^{2}$. For the Heston SDEs we now extend this result to convergence in $L^{q}$ for all $q\geq 2$, with estimates of the $L^{q}$ speeds of convergence.

Theorem 1**.**

Fix any starting points $V_{0}>0$ and $R_{0}$ for the returns rate $R_{t}$ and squared volatilities $V_{t}$ driven by Heston SDEs (2), (3). Fix $T>0$ and any even integer $q\geq 2$. Then there is a constant $c(q)$ such that, for any $\varepsilon>0$ and the choice of partition sizes $J(\varepsilon)$, and any $t<T$, the realized volatilities $Y_{t}^{\varepsilon}$ defined by (16) verify

[TABLE]

Hence, when $\varepsilon\to 0$ and $J(\varepsilon)\to\infty$, the $Y_{t}^{\varepsilon}$ converge to $V_{t}$ in $L^{q}$, uniformly over $0\leq t\leq T$. Moreover if one imposes $J(\varepsilon)>(a/\varepsilon)^{q/2}$ for some fixed $a>0$, one has then

[TABLE]

Proof.

Recall a well known lemma, easily proved by recurrence.

Lemma 1**.**

For each integer $q\geq 1$ and any random variables $W_{1},\ldots,W_{q}$ in $L^{q}$, one has

[TABLE]

Fix $T>0,V_{0}>0,R_{0}$. By (5) and (6), for each $q\geq 1$, there is a constant $C_{q}$ such that for $0\leq s\leq t\leq T$ one has

[TABLE]

Assume first that parameter $\mu=0$ in the Heston SDE (2) for the returns. Then $R_{t}-R_{s}=\int_{s}^{t}\sqrt{V_{u}}dZu$, and by Ito formula,

[TABLE]

Hence, the variables $D(s,t)$ defined by

[TABLE]

must verify

[TABLE]

Moreover (19) gives

[TABLE]

and hence (20) yields

[TABLE]

with $b_{q}=C_{2q}^{2}+C_{q}$. For each $\varepsilon>0$, select the partition size $J(\varepsilon)$, and partition the sliding window $[t-\varepsilon,t]$ by the time points $t_{n}=t-\varepsilon+n\varepsilon/J$, with $n=0,\ldots,J$. Recall formula (16) for the realized volatilities $Y_{t}^{\varepsilon}$

[TABLE]

To study $Y_{t}^{\varepsilon}-V_{t}$ we introduce the decomposition

[TABLE]

where the terms $H(t,\varepsilon)$ and $K(t,\varepsilon)$ are defined as

[TABLE]

We can rewrite $K(t,\varepsilon)$ as

[TABLE]

and this implies

[TABLE]

Equation (19) gives the bound $\|V_{t}-V_{u}\|_{q}\leq C_{q}(t-u)^{1/2}$ for $u\leq t\leq T$, and hence

[TABLE]

We now study $H(t,\varepsilon)$. Define $U_{n}$ for $n=1,\ldots,J$ by

[TABLE]

Formula (25) then implies directly

[TABLE]

From (23) and (27) we obtain

[TABLE]

Define the polynomial $Q=Q(U_{1},\ldots,U_{n})$ by

[TABLE]

Since $q$ is even, we then get, due to (28),

[TABLE]

Next, we derive the bound on ${\mathbb{E}}[Q]$ with $Q$ defined in (30). To this end, we first define the set $M\equiv M(q,J)$ of all multi-integers

[TABLE]

For any $\vec{m}\in M$ denote by $Q_{\vec{m}}$ the monomial $Q_{\vec{m}}=U_{m(1)}U_{m(2)}\ldots U_{m(q)}$. Then, we can expand the polynomial $Q$ as follows

[TABLE]

Then (31) is equivalent to

[TABLE]

Due to lemma 1 and (29), we have for all $\vec{m}\in M$,

[TABLE]

The above bound is sufficient for most $\vec{m}\in M$ but needs to be refined on a specific subset of $M$. So for each $\vec{m}\in M$, let $m^{*}=\max(m(1),\ldots,m(q))$. Let $z(\vec{m})$ be the number of indices $m(k)$ such that $m(k)=m^{*}$. For $1\leq r\leq q$, call $M_{r}$ the set of $\vec{m}\in M$ such that $z(\vec{m})=r$. Then $M$ is the union of disjoint subsets $M_{r}$, $r=1,\ldots,q$. Each $\vec{m}\in M_{r}$ contains at most $q-r+1$ distinct indices $m(k)$. Hence each $M_{r}$ has cardinal $Card(M_{r})\leq J^{q-r+1}$. We now consider two cases separately (i) $r\geq 2$ and (ii) $r=1$.

For $r\geq 2$ we can have a an upper bound $Card(M_{r})\leq J^{q-1}$, and therefore

[TABLE]

This yields, in view of (33),

[TABLE]

We now consider $M_{1}$ separately. Fix any $\vec{m}\in M_{1}$. The maximum $m^{*}=\max(m(1)\ldots m(q))$ is then reached by a single index $i^{*}$ such that only $m(i^{*})=m^{*}$ and $m(i)<m^{*}$ for $i\neq i^{*}$. This implies $2\leq m^{*}\leq J$ since $q\geq 2$. We can then re-order $\vec{m}$ as a multi-index $\vec{\nu}$ verifying

[TABLE]

Let $s=t_{j-1}$ and $t=t_{j}$. For $1\leq k\leq q-1$ all the $U_{\nu(k)}$ are ${\cal G}_{s}$-measurable, so that

[TABLE]

Due to (21) and definition (27), we have

[TABLE]

and, therefore, ${\mathbb{E}}[Q_{\vec{m}}|{\cal G}_{s}]=0$ for $\vec{m}\in M_{1}$.

Thus, one has

[TABLE]

and equation (34) implies

[TABLE]

Equations (31) and (35) then yield the bound

[TABLE]

Combining equation (24) with the two bounds (36) and (26), we conclude that for all even $q\geq 2$, all $t\leq T$, and all $\varepsilon>0$, one has

[TABLE]

Hence when $\lim_{\varepsilon\to 0}J(\varepsilon)=+\infty$, we obtain the $L^{q}$ convergence

[TABLE]

and this convergence is uniform for $0\leq t\leq T$. To minimize the upper bound on the speed of convergence given by (37), one must clearly impose $J(\varepsilon)\sim 1/\varepsilon^{q/2}$, and then as $\varepsilon\to 0$ one has

[TABLE]

For $\mu\neq 0$, the inequalities on the return process presented in the proof above still hold since we’re considering a fixed time $T>0$ and do not strive to have uniform bound in $T$. Consider, for instance, (22) with $\mu\neq 0$

[TABLE]

and using (12) we obtain the bound

[TABLE]

This concludes the proof.

4 Observable Estimators for the Heston Model Parameters

4.1 Parameter Estimation from true volatility data

To fit the Heston model to asset prices data, one needs to estimate the parameters $\mu,\rho$ of the SDE (2) and the parameter vector $\boldsymbol{\theta}$ of the Heston volatility SDE (3). Since the volatility $V_{t}$ is unobservable, the key issue in parametric estimation of the Heston model is to estimate $\boldsymbol{\theta}$ (see, e.g. [9]). Consider first the ideal but unrealistic case where we are given a large finite set of $N$ true volatilities values $\mathcal{V}=\{V_{n\Delta}\}$, sub-sampled at time intervals $\Delta$. For processes driven by smoothly parametrized SDEs, many publications have studied parameter estimation from large data sets actually generated by the underlying SDEs (see for instance [1, 3, 2, 14, 25, 28, 27, 37, 34, 15], etc.). Several of these approaches rely either on Maximum Likelihood Estimators (MLEs) or on Methods of Moments.

Maximum Likelihood Estimators:

For the Heston volatility SDE, the MLEs of $\boldsymbol{\theta}$ have been thoroughly analyzed in [9], where they are explicitly computed from any finite set of true squared volatilities $\mathcal{V}$. Under minor parameter constraints and as $N\to\infty$, these MLEs were shown to be asymptotically consistent, and asymptotically normal when $\kappa\theta/\gamma^{2}>1$. Note that the impact of replacing the unobservable volatilities $V_{t}$ by the realized volatilities $Y_{t}^{\varepsilon}$ in the explicit MLE formulas of [9] is a quite technical task which we will complete in a future paper.

Moments based Estimators:

In this paper, we will focus instead on natural Moment Estimators $\hat{\boldsymbol{\theta}}$ of the Heston SDE parameter vector $\boldsymbol{\theta}$. Since the true squared volatilities $V_{t}$ are unobservable, $\hat{\boldsymbol{\theta}}$ is constructed as an explicit smooth function of the empirical mean and two lagged empirical covariances of the observable realized volatilities $Y_{t}^{\varepsilon}$.

4.2 Parameter estimation under indirect observability

The Moments Estimators approach considered in this paper falls formally within the generic Indirect Observability framework we introduced in [11]. In this framework, we analyze the observable processes $Y_{t}^{\varepsilon}$ which, as $\varepsilon\to 0$, converge in $L^{4}$ to an unobservable processes $X_{t}$ parametrized by a vector $\boldsymbol{\theta}$. In particular, in [11], after selecting a number of observables $N(\varepsilon)$ and a sub-sampling rate $\Delta(\varepsilon)$, the observable estimators $\hat{\boldsymbol{\theta}}(\varepsilon)$ of $\boldsymbol{\theta}$ are constructed as smooth functions of the empirical mean and a finite set of empirical lagged covariances of the observables $Y_{n\Delta(\varepsilon)}^{\varepsilon},1\leq n\leq N(\varepsilon)$. Under a broadly applicable set of Indirect Observability Hypotheses, which however require $X_{t}$ to be weakly stationary, we proved in [11] that one can construct observable moment estimators achieving consistency as $\varepsilon\to 0$, provided $N(\varepsilon)$ and $\Delta(\varepsilon)$ are adequately selected.

Here we focus on the following indirect observability situation:

(i) the unobservable process $X_{t}$ is the squared volatility process $V_{t}$

(ii) the observable $Y_{t}^{\varepsilon}$ converging to $V_{t}$ as $\varepsilon\to 0$ are the realized volatilities defined by the rate of returns process associated to $V_{t}$.

Note that in the present paper the volatility process $V_{t}$ starting at a deterministic $V_{0}=y>0$ is not stationary; therefore, the analytical results of the present paper have requires several quite technical enhancements of the methods previously used in [11].

5 Transition Densities for squared volatilities

5.1 Explicit transition density

Consider squared volatilities $V_{t}$ driven by the Heston volatility SDE (3) parametrized by $\boldsymbol{\theta}$. We will always assume that $V_{0}>0$ has finite moments of all orders, which is of course the case if $V_{0}$ is deterministic. For $T>0$, introduce the following short-hand notations

[TABLE]

As shown in [23], the Markov diffusion process $V_{t}$ has an explicit transition density $p_{T}(z,y)$, which we often denote $p(z,y)$ for short, given by

[TABLE]

where $I_{r}$ is the modified Bessel function of the 1st kind of order $r$. As noted in [23], for fixed $T$, the linear rescaling $V_{t}\to 2\lambda_{T}V_{t}$ transforms the transition density $p_{T}(z,y)$ into

[TABLE]

which, for each fixed $y$, is a non-central $\chi^{2}$ density with non-centrality parameter $NCP(T,y)=y\nu_{T}$ and $DFR=2r+2$ degrees of freedom. Note that $DFR={4\kappa\theta}/{\gamma^{2}}=2\theta\Lambda$.

5.2 The stationary squared volatility process $V_{t}$

Since $\nu_{T}\to 0$ and $\lambda_{T}\to\Lambda$ as $T\to\infty$, $p_{T}(z,y)$ converges pointwise, at the speed $e^{-\kappa T}$, to the unique stationary density $\psi(z)$ of the autonomous Heston volatility SDE. This stationary density is given for all $z>0$ by

[TABLE]

When the initial condition $V_{0}$ is random with density $\psi$, all $V_{t}$ have then the same density $\psi$, and the process $V_{t}$ driven by the Heston volatility SDE becomes strictly stationary. Expectations with respect to this stationary diffusion will be denoted ${\mathbb{E}}_{\psi}$.

Note, that since $\lim_{T\to\infty}NCP(T,y)=0$, the linear rescaling $z\to 2\Lambda z$ transforms the stationary density $\psi(z)$ into $({2\Lambda})^{-1}\psi(z/2\Lambda)$ which is a standard $\chi^{2}$ density with $DFR=2r+2$ degrees of freedom.

6 Conditional Moments of squared volatilities

6.1 Moments of non-central $\chi^{2}$

Let $X$ be a random variable having a non-central $\chi^{2}$ density with $DFR$ degrees of freedom and non-centrality parameter $NCP$. Then, the Laplace transform $Lap(z)={\mathbb{E}}(e^{zX})$ is, for $0\leq z<1/2$,

[TABLE]

Expanding $1/(1-2z)^{n}$ as a series in $z$, we obtain, for a fixed $DFR$,

[TABLE]

where the $\pi_{q}(NCP)$ are polynomials of degree $q$ in $NCP$, with coefficients fully determined by $DFR$ and $q$. Denote $nc\chi(q)$ and $st\chi(q)$ the respective moments of order $q$ for the non-central $\chi^{2}$ density and for the standard $\chi^{2}$ density with $DFR$ degrees of freedom. We then have the polynomial expressions

[TABLE]

For the first two moments of the non-central $\chi^{2}$ and the standard $\chi^{2}$, one has, for instance, the following well known formulas

[TABLE]

6.2 Conditional Moments of the squared volatilities

Recall that $\nu_{T}=e^{-\kappa T}$ and that $DFR=2r+2$ is determined by $\boldsymbol{\theta}$. The next proposition addresses the computation of conditional moments for $V_{t}$

[TABLE]

Proposition 2**.**

For each $q\geq 1$ and for each $y>0$, the conditional moments ${\mathbb{E}}_{y}[V_{T}^{q}]$ remain uniformly bounded for all $T\geq 0$. There is a polynomial $Q_{q}$ with coefficients depending only on $q$ and $\boldsymbol{\theta}$, such that for all $s$, $T$ and all $y>0$,

[TABLE]

As $T\to\infty$, moments $M_{q}(y,T)$ converge at the exponential speed $\nu_{T}=e^{-\kappa T}$ to finite moments $m_{q}={\mathbb{E}}_{\psi}[V_{t}^{q}]$ of the stationary diffusion $V_{t}$ driven by the Heston volatility SDE.

Proof.

The rescaling $V_{s}\to 2\lambda_{T}V_{s}$, with $\lambda_{T}$ in (38), transforms the conditional distribution of $V_{s+T}$ given that $V_{s}=y$ into a non-central $\chi^{2}$ with

[TABLE]

This rescaling implies, using the non-central $\chi^{2}$ moments (41),

[TABLE]

Define the homogeneous polynomial $H_{q}(a,b)$ of total degree $q$ by

[TABLE]

where $\pi_{q}(\cdot)$ is defined by the Laplace transform introduced in section 6.1. Therefore, coefficients of $H_{q}$ depend only on $q$ and on $\boldsymbol{\theta}$. Next, we define

[TABLE]

Since $0\leq\nu_{T}\leq 1$ is a constant given by (38), the expression (46) for $Q_{q}(y)$ proves equation (44). Equation (44) also implies the existence of a constant $C$ such that

[TABLE]

Therefore for each $V_{0}=y>0$, moments ${\mathbb{E}}_{y}[V_{T}^{q}]$ remain bounded for all $T\geq 0$. As discussed in section 5.2, rescaling by $\Lambda$ transforms the stationary density $\psi$ into a standard $\chi^{2}$ distribution, and hence stationary moments

[TABLE]

must be finite. As $T\to\infty$, both $\nu_{T}=e^{-\kappa T}$ and $NCP$ tend to [math], and $\lambda_{t}\to\Lambda$ while DFR remains constant. Hence, due to (44), (46) the $M_{q}(y,T)$ converge at exponential speed $\nu_{T}$ to

[TABLE]

6.3 Mean and Covariances of the stationary diffusion $V_{t}$

Using (42), and the appropriate rescaling of $V_{t}$ by $2\lambda_{T}$ one easily computes the first two conditional moments of the squared volatility process starting at $y>0$, namely

[TABLE]

where ${\mathbb{E}}_{y}[\cdot]$ is the conditional moment defined in (43). In particular, as $T\to\infty$, equations (47) and (48) yield the first two moments of the stationary diffusion $V_{t}$

[TABLE]

Moreover, the stationary diffusion driven by the Heston volatility SDE has mean $m_{1}=\theta$ and lagged covariances $K(u)=cov_{\psi}[V_{s}V_{s+u}]$ given by

[TABLE]

for any time lag $u\geq 0$. This yields the stationary covariances

[TABLE]

6.4 Heston SDE parameters as functions of asymptotic moments

We can now express $\boldsymbol{\theta}=(\kappa,\theta,\gamma)$ as an explicit smooth function

[TABLE]

of three moments of the stationary volatility diffusion $V_{t}$, namely its mean $m_{1}$, its variance $K(0)$, and one lagged covariance $K(u)$ for some fixed (but arbitrary) $u>0$. Equations (49) and (50) indeed imply that parameters $(\kappa,\theta,\gamma)$ can be expressed using the moments $m_{1}$, $K(0)$, and $K(u)$ as follows

[TABLE]

which defines the function $\Phi$ above.

7 Moments based observable estimators

We now use our preceding results on the Heston volatility SDEs to study a class of moment-based estimators of the Heston parameters and to discuss their consistency when the observable data are generated by the realized volatilities.

7.1 Computation of Moments Based Observable Estimators

Given the window-size $\varepsilon>0$, select a sub-sampling time interval $\Delta\equiv\Delta(\varepsilon)$ and a number of observations $N\equiv N(\varepsilon)$. Then, the realized volatility process (16) generates an observable data set of $N(\varepsilon)$ realized volatilities

[TABLE]

Next, we specify how we use these $N(\varepsilon)$ observable data to estimate any lagged covariance $K(u)$ of the stationary diffusion $V_{t}$. Denoting $\left[a\right]_{int}$ the closest integer to $a$, we approximate the lag $u$ by $U\Delta(\varepsilon)$ where

[TABLE]

Since $K(u)$ is Lipschitz in $u$, there is a constant $C\equiv C(u)$ such that

[TABLE]

Then, for any $\varepsilon$ and time lag $u$, we define observable empirical estimators of the mean $m_{1}$ and lagged covariances $K(u)$ of the stationary diffusion $V_{t}$ as follows

[TABLE]

where $U=U(u,\varepsilon)$ and $N=N(\varepsilon)$. Formulas (51) express the parameter vector $\boldsymbol{\theta}$ as an explicit $C^{1}$ function $\Phi(m_{1},K(0),K(u))$. This naturally leads to the definition of an observable parameter estimator $\hat{\boldsymbol{\theta}}(\varepsilon)$ of $\boldsymbol{\theta}$ by

[TABLE]

This definition yields the following explicit observable estimators of the Heston parameters

[TABLE]

7.2 Asymptotics for polynomial functionals of squared volatilities

Theorem 2**.**

Consider any fixed polynomial $h(x_{1},\ldots,x_{k})$ of total degree $n$ in $k$ variables $(x_{1},\ldots,x_{k})$. Let $0=u(0)<u(1)<\ldots<u(k)$ be any sequence of $k+1$ lag instants. For $T>0$, define $H$ and $H_{T}$ by

[TABLE]

Recall that $\nu_{T}=e^{-\kappa T}$. Define $w_{j}=e^{-\kappa(u(j+1)-u(j))}$ for $j=0,\ldots,k-1$. Then, there is a polynomial $POL$ in $k+2$ variables such that for all $T>0$ and all $y>0$

[TABLE]

The degree and coefficients of POL are determined by the integers $n$, $k$, the coefficients of $h$, and the vector $\boldsymbol{\theta}$. The asymptotic polynomial moments are then given by

[TABLE]

For any integer $q\geq 1$ there is a positive constant $C$, and an integer $p\geq 1$, determined only by $q$, $k$, $\boldsymbol{\theta}$ and the polynomial $h$ such that, for all positive $T$ and $y$, and all $0=u(0)<u(1)<\ldots<u(k)$

[TABLE]

In particular for $q=1$ one has

[TABLE]

Proof.

For better readability, the detailed proof is given in Appendix A.

Remarks. Equation (56) also implies that as $T\to\infty$, the random polynomial functions $H_{T}$ converge in $L^{q}$-norm to the constants ${\mathbb{E}}_{\psi}(H)$, where $L^{q}$-norms are computed under ${\mathbb{E}}_{y}$. Note that the constant $C$ introduced in the theorem does not depend on the time lags $u(0)<u(1)<\ldots<u(k)$.

7.3 Consistency of observable estimators

Since $\hat{\boldsymbol{\theta}}(\varepsilon)$ is a $C^{1}$ function of three specific empirical moments of realized volatilities, the key consistency issue is to estimate, as $\varepsilon\to 0$, the speeds of convergence of $\hat{m}^{\varepsilon}$ to $m_{1}$ and $\hat{K}^{\varepsilon}(u)$ to $K(u)$. These speeds of convergence strongly depend on the sub-sampling scheme defined by $N(\varepsilon)$ and $\Delta(\varepsilon)$. In [11], we have determined sub-sampling schemes optimizing these speeds of convergence for stationary unobservable limit processes. We now prove similar results for the non-stationary volatilities driven by the Heston SDEs.

Theorem 3**.**

Consider an asset with return rate $R_{t}$ and squared volatility $V_{t}$, jointly driven by the Heston SDEs (2), (3). Fix deterministic initial conditions $R_{0}$ and $V_{0}=y>0$. Call $P_{y}$ the probability distribution in path space of the trajectories $\{R_{t},V_{t}\}$. Realized volatilities $Y_{t}^{\varepsilon}$ are computed by formula (16) with $J(\varepsilon)\sim\varepsilon^{-2}$. The $Y_{t}^{\varepsilon}$ are sub-sampled with time step $\Delta(\varepsilon)$ to generate $N(\varepsilon)$ observations $W_{k}=Y_{k\Delta(\varepsilon)}^{\varepsilon}$. We then apply formulas (53) to compute observable empirical estimators $\hat{K}^{\varepsilon}(u)$ and $\hat{m}^{\varepsilon}$ of the asymptotic lagged covariances $K(u)-\theta^{2}=\lim_{t\to\infty}{\mathbb{E}}[V_{t}V_{t+u}]$ and mean $m_{1}=\lim_{t\to\infty}{\mathbb{E}}[V_{t}]$ of true volatilities.

Then there exists a sub-sampling scheme which guarantees that for any fixed positive $L$ and $y$ there is a constant $C=C(L,y,\boldsymbol{\theta})$ such that, for all lags $0\leq u\leq L$, one has

[TABLE]

where $L^{2}$-norms are computed with respect to $P_{y}$. Moreover, under $P_{y}$ the observable parameters estimators $\hat{\boldsymbol{\theta}}^{\varepsilon}$ given by formulas (54) converge in probability to the true Heston parameters $\boldsymbol{\theta}$ as $\varepsilon\to 0$. One has, for an adequate constant $C$,

[TABLE]

Proof.

Fix the time lag $u$ and $V_{0}=y>0$. All $L^{q}$-norms are computed under $P_{y}$. The notation “constant $C$” will designate a generic constant which can change values from one bound to another. By Theorem 1, there is a constant $c_{4}$ such that for all $t$,

[TABLE]

Therefore, there is a constant $c_{2}$ such that for all $s$ and $t$,

[TABLE]

Denote $\Delta\equiv\Delta(\varepsilon)$. The sub-sampled realized volatilities $W_{k}=Y_{k\Delta}^{\varepsilon}$ determine the observable empirical estimators of 1st and 2nd moments of volatilities, through formula (53). Since $\|V_{k\Delta}-W_{k}\|_{2}\leq c_{4}\varepsilon^{1/2}$ by (60), the definition (53) of $\hat{m}^{\varepsilon}$ gives

[TABLE]

where we used (50) which implies that

[TABLE]

We would like to point out that the term $O(1/\sqrt{N\Delta})$ in the expression above arises from the $L^{2}$ error of the empirical mean, $m_{1}$, computed form direct observations, $V_{k\Delta}$. Therefore, the estimate in (61) cannot be improved analytically.

Provided we take $N\Delta=\varepsilon^{-1}$, this proves convergence, at speed $\varepsilon^{1/2}$, of the empirical mean of realized volatilities $\hat{m}^{\varepsilon}$ to the asymptotic mean of $V_{t}$. Next, we study our estimators of lagged covariances. Let $U$ be the closest integer to $\left[u/\Delta(\varepsilon)\right]_{int}$, so that $|U\Delta-u\Delta|\leq\Delta$ . Define

[TABLE]

By subadditivity of norms, we obtain

[TABLE]

Let $\psi$ be the asymptotic density of $V_{t}$, and denote by

[TABLE]

the stationary lagged 2nd moments, which do not depend on $t$. By Theorem 2, for every fixed $a$ and $y>0$ there is a constant $C$ such that for all $T$ and all $s\leq a$ one has the bound

[TABLE]

The $L^{2}$ norm under $P_{y}$ hence verifies

[TABLE]

The above two bounds in (63) and (64) provide constants $C$ and $c=\kappa/2$ such that for all $\varepsilon$

[TABLE]

By subadditivity of $L^{2}$ norms, inequality (65) then implies,

[TABLE]

Regrouping our definitions and notations, we have

[TABLE]

This implies since $K(u)$ is Lipschitz in $u$,

[TABLE]

We have the obvious identity

[TABLE]

and hence

[TABLE]

We now use the bounds (61), (62), (66), and (67) to obtain

[TABLE]

To optimize this last bound and ensure that all terms have the same rate of convergence as $\varepsilon\to 0$, we impose the choice

[TABLE]

which is equivalent to selecting $\Delta\sim\varepsilon^{1/2}$ and $N\sim\varepsilon^{-3/2}$.

Therefore, for each fixed $V_{0}=y>0$ and for all time lags $u$ within any fixed interval $\left[0,L\right]$ there is a constant $C$ realizing the bound

[TABLE]

The $L^{2}$ convergence under $P_{y}$ of $\hat{K}^{\varepsilon}(u)$ to $K(u)$ and of $\hat{m}^{\varepsilon}$ to $m_{1}$, implies their convergence in probability under $P_{y}$. By equation (54) our estimators of Heston parameters are of the form

[TABLE]

where $\Phi$ is a $C^{1}$ function. Thus, estimators $\hat{\boldsymbol{\theta}}^{\varepsilon}$ converge in probability to $\boldsymbol{\theta}=$ $\Phi(K(0),K(u),m_{1})$ as $\varepsilon\to 0$. The $L^{2}$-speeds of convergence $\varepsilon^{1/2}$ for the 1st and 2nd moments imply, by Chebyshev inequality,

[TABLE]

with a similar inequality for $\hat{m}^{\varepsilon}$. Since $\Phi$ is $C^{1}$, these speeds of convergence in probability under $P_{y}$ imply, by the first order Taylor expansion of the function $\Phi$, the same speed of convergence in probability for the parameter estimators themselves.

8 Asymptotically Optimal Partition Sizes $J(\varepsilon)$

Denote by $T(\varepsilon)$ the total observation time available for the rate of returns process $R_{t}$ The realized volatilities $Y^{\varepsilon}_{t}$ given by formula (16) involve subdividing the sliding window $(t-\varepsilon,t)$ into $J(\varepsilon)$ intervals and averaging the corresponding $J(\varepsilon)$ squared increments of the returns rate. To compute our 1st and 2nd moments estimators from the observable process $Y^{\varepsilon}_{t}$ with $t\leq T(\varepsilon)$, we sub-sample this process at the $N(\varepsilon)$ instants $k\Delta(\varepsilon)$, with $k=1,\ldots,N(\varepsilon)$, with the obvious relation $N(\varepsilon)\Delta(\varepsilon)=T(\varepsilon)$.

Provided one uses the subsampling scheme

[TABLE]

and a partition size $J(\varepsilon)\sim 1/\varepsilon^{2}$, our current theoretical bounds can guarantee $L^{4}$ speeds of convergence $\sim\sqrt{\varepsilon}$ for $Y^{\varepsilon}_{t}-V_{t}$ and consistency in probability of our observable estimators for the parameters in the volatility Heston SDE. The theoretical choice $J(\varepsilon)\sim 1/\varepsilon^{2}$ seems overwhelmingly large for concrete fitting of Heston SDEs to actual stockprices data. Therefore, we also examine numerically a more pragmatic choice $J(\varepsilon)\sim\varepsilon^{-1}$.

Remark. The optimized sub-sampling scheme (70) necessitates a total observational time $T(\varepsilon)\sim\varepsilon^{-1}$. The associated convergence rate (58) for the 1st and 2nd moments estimators can hence be reformulated as

[TABLE]

However, to compute realized volatilities $Y_{t}^{\varepsilon}$, we need $J(\varepsilon)\sim 1/\varepsilon^{2}$ time points in each sliding window, and hence the computation of the observable moments estimators requires a total number of observational points

[TABLE]

Therefore, the convergence rates of our moments estimators can be expressed as follows in terms of the total number $n$ of observational time points for the return process as

[TABLE]

Even the more pragmatic choice $J(\varepsilon)=1/\varepsilon$ results in the scaling

[TABLE]

These upper bounds on convergence rates are suboptimal for the somewhat theoretical case when the available observational time is not a priori bounded. Indeed Hoffmann [32] indicates that when the total observational time $T\to\infty$ one should expect a more classical convergence rate $n(\varepsilon)^{-1/2}$, while for fixed finite $T$, the paper [32] suggests that the optimal convergence rate should be $n(\varepsilon)^{-1/4}$. However, results in [32] focus on approximate maximum likelihood estimators, and, therefore, cannot be directly applied to the moments based estimators which are considered in this paper.

9 Effective Speed of Convergence of realized volatilities to true volatilities

9.1 Generic stochastic volatility models versus Heston SDEs

Recall that realized volatilities $Y_{t}^{\varepsilon}$ are computed by formula (16) with partition size $J(\varepsilon)$ for the time windows $\left[t-\varepsilon,t\right]$. When the rates of return $R_{t}$ and the squared volatilities $V_{t}$ are driven by joint Heston SDEs, we have proved in Theorem 1 that for each fixed even integer $q$ and for $s$ bounded, the $L^{q}$ norms $\|Y_{s}^{\varepsilon}-V_{s}\|_{q}$ must verify, for some constant $C=C(q)$, the bounds

[TABLE]

Our numerical simulations suggest that for the ”moderate” partition size $J(\varepsilon)\sim 1/\varepsilon$, one can possibly improve (71) to yield the following convergence speeds, valid for $q=2$, $4$ and $s$ bounded,

[TABLE]

For $q=2$, this is indeed implied by (71). However, for $q=4$, our theoretical bound (71) seems to overestimate the size of the partition $J(\varepsilon)$ required to achieve the $L^{4}$ speed of convergence $\sim\sqrt{\varepsilon}$

In this paper, to validate numerically the effective $L^{2}$ and $L^{4}$ speeds of convergence of realized volatilities $Y_{t}^{\varepsilon}$ to true volatilities $V_{t}$ we have carried out the following intensive simulations for joint Heston SDEs.

9.2 Outline of our Heston SDEs simulations

We have numerically simulated the joint Heston SDEs with the following specific parameters

[TABLE]

and for 3 values $\rho=0,~{}0.3,~{}0.7$ of the correlation coefficient between the Brownian noises driving the joint Heston SDEs. The Feller condition is valid since $2\kappa\theta/\gamma^{2}=3.4$. To emulate asymptotics as $\varepsilon\to 0$, we consider the partition sizes

[TABLE]

and five values of $\varepsilon=0.1$, $0.05$, $0.02$, $0.01$, $0.005$. Simulations with the fixed partition sizes $J=10,\,40$ are presented here to illustrate that the errors in the approximation of volatility by the realized volatility will not decay as $\varepsilon\to 0$ if the partition has a fixed number of points. Numerical simulations with $J=\varepsilon^{-1},\,\varepsilon^{-2}$ are more interesting since they provide an insight into the convergence rate and selecting the optimal sub-sampling regime for parameter estimation under indirect observability.

Simulations of true volatility paths $V_{t}$ are implemented by an Euler dicretization scheme for SDEs, with time step $10^{-6}$, except for $\varepsilon=0.005$, $J(\varepsilon)=1/\varepsilon^{2}$, where the time step was $1.25\times 10^{-7}$. We perform a Monte-Carlo simulation by generating $200,000$ independent simulated paths $\{V_{t},Y_{t}^{\varepsilon}\}$. We then partition these $200,000$ paths into sub-ensembles, as discussed in section 9.4.

9.3 Snapshots of joint sample paths $\{V_{t},Y_{t}^{\varepsilon}\}$

For $\varepsilon=0.01$ and $\rho=0$, Figure 1 displays four examples of joint paths $\{V_{t},Y_{t}^{\varepsilon}\}$, where realized volatilities $Y_{s}^{\varepsilon}$ are successively computed with the four $J(\varepsilon)$ listed in (74). Clearly, the accuracy of the approximation of $V_{t}$ by $Y_{t}^{\varepsilon}$ increases drastically for larger partition sizes $J(\varepsilon)$. The smallest $J(\varepsilon)$, equal to 10, generates many quite significant inaccuracies for $Y_{t}^{\varepsilon}-V_{t}$. For $J=40$, we still note several significant inaccuracies. But for $J(\varepsilon)=10000$ the sample paths of $V_{t}$ and $Y_{t}^{\varepsilon}$ nearly coincide. Such large partition sizes are generally not feasible: for intraday stock prices sampled every minute, a partition size $J=10000$ would require an unrealistic sliding time window of about 20 trading days; for stock prices sampled every second, such partition would require a sliding window of approximately 2.7 hours. For small values of $J(\varepsilon)$, a practical remedy to eliminate large sharp peaks of $|Y_{t}^{\varepsilon}-V_{t}|$ is time smoothing of the $Y_{t}^{\varepsilon}$ either directly, or by using a weighted average in (16).

9.4 Numerical Asymptotics of $\|Y_{t}^{\varepsilon}-V_{t}\|_{q}$ as $\varepsilon\to 0$

We partition all $200,000$ simulated paths into sub-ensembles of $1000$ paths resulting in $G=200$ of such sub-ensembles. This allows us to compute confidence intervals for the numerically estimated $L^{q}$ errors between the volatility and realized volatility. In particular, we fix $T=1$, and for each sub-ensemble the empirical mean $M_{k}(q)$ of $|Y_{T}^{\varepsilon}-V_{T}|^{q}$ , $k=1,\ldots,G$ provides an estimator $\hat{e}^{q}_{k}=\left(M_{k}(q)\right)^{1/q}$ for the $L^{q}$ errors $\|Y_{T}^{\varepsilon}-V_{T}\|_{q}$. Final estimates for these $L^{q}$ errors are then given by

[TABLE]

with 95% confidence intervals $\hat{E}^{q}\pm 1.96\sigma(q)$ where

[TABLE]

Numerical results for the $L^{2}$ and $L^{4}$ convergence are presented in Table 1 and Figure 2. For constant partition sizes $J(\varepsilon)=10$, $40$ both $\hat{E}^{2}$ and $\hat{E}^{4}$ error estimates are nearly constant (independent of $\varepsilon$), as predicted by our analytical bound (71). Moreover, $\hat{E}^{2}$ as well $\hat{E}^{4}$ errors are both approximately twice smaller for $J=40$ than for $J=10$. For $q=2$, this is correctly predicted by our theoretical bound (71). But for $q=4$, our theoretical bound is too pessimistic, since it predicts that $\hat{E}^{4}$ should be about $1.4$ times smaller for $J=40$ than for $J=10$.

Figure 2 depicts the $L^{2}$ and $L^{4}$ errors on the log-log scale for partition sizes $J(\varepsilon)=\varepsilon^{-1}$ and $J(\varepsilon)=\varepsilon^{-2}$. The graphs of estimation errors in these two situations are nearly perfect straight lines with slope 1/2 as soon as $\varepsilon$ is small enough. Figure 2 demonstrates quite convincingly that the two types of estimation errors $\hat{E}^{2}$ and $\hat{E}^{4}$ do scale like $\varepsilon^{1/2}$ for $J(\varepsilon)\sim 1/\varepsilon$ as well as for $J(\varepsilon)\sim 1/\varepsilon^{2}$. So our numerical simulations of the joint Heston SDEs support the following conjecture about the asymptotic behaviour (as $\varepsilon\to 0$) of the $L^{4}$ error

[TABLE]

Our simulations indicate that for $q=2$ and $q=4$ convergence speeds $\|Y_{t}^{\varepsilon}-V_{t}\|_{q}\sim\sqrt{\varepsilon}$ can be achieved for fixed $t$ when the realized volatility $Y_{t}^{\varepsilon}$ is computed with partition sizes $J(\varepsilon)\sim 1/\varepsilon$. This also implies that for the partition size $J(\varepsilon)\sim 1/\varepsilon$, the lagged covariances of realized volatilities should converge to true lagged covariances at $L^{2}$-speeds $\sim\sqrt{\varepsilon}$. We would like to point out that these results are obtained for finite $\varepsilon\geq 0.01$. It is extremely time-consuming to extend these results for smaller values of $\varepsilon$ and it is possible that the asymptotic behavior of $L^{q}$ errors might change for $\varepsilon\ll 0.01$. Surprisingly, this sub-sampling scheme with $J\sim 1/\varepsilon$ gives much fast convergence rate for the $L^{4}$ norm compared to our analytical estimates. As discussed above, it is possible that for extremely small values of $\varepsilon$ numerical simulations would become consistent with our analysis and yield the convergence rate of $\varepsilon^{1/4}$. Improved convergence speed might be related to the ratio of constants in our analytical estimated for the $L^{q}$ speed of convergence for the moment estimates However, for such small values of $\varepsilon$ numerical computations becomes extremely costly. For practical values of $\varepsilon$ considered here we obtain estimated convergence rate of $\varepsilon^{1/2}$.

We also performed numerical simulations with $\rho=0.3$ and $\rho=0.7$ (not displayed here for brevity), where $\rho$ is the correlation between the two Brownian Motions $B_{t}$ and $Z_{t}$ driving the joint Heston SDEs. Our numerical results with $\rho>0$ are almost identical to those for $\rho=0$. This is consistent with our proof of Theorem 1, which explores the autonomous Heston SDE (3) driving the true volatility $V_{t}$, without ever using the Heston SDE (2) for the rate of return process. Another key ingredient of our proof is the study of conditional expectations ${\mathbb{E}}(Y|X)$ when $X$ and $Y$ are polynomial functions of a finite number of $V_{t}$ values. Again this analysis does not use the Heston SDE (2). Constants introduced in Theorem 1 may depend on $\rho$, but our numerical simulations indicate that this dependence is fairly weak.

10 Effective convergence speeds for observable estimators of Heston parameters

In this section we evaluate numerical convergence speeds for our observable estimators $\hat{\theta}_{\varepsilon}$, $\hat{\kappa}_{\varepsilon}$, $\hat{\gamma}_{\varepsilon}$ of the Heston volatility SDE parameters. Recall that these estimators are based on estimated covariances of realized volatilities. This set of simulations is performed as outlined in section 9.2 with the following four values of $\varepsilon=0.1$, $0.05$, $0.02$, $0.01$. Realized volatilities are computed with two different partition sizes

[TABLE]

In order to compute estimators we use the sup-sampling regime

[TABLE]

which is a particular case of our general regime in (70).

Numerical estimates for the $L^{2}$ errors of estimation for the Heston SDE parameters are computed using a Monte-Carlo approach with 1000 long trajectories consistent with the sub-sampling regime outlined above. Each long trajectory yields one set of estimated parameter values computed using (54).

The lag $u^{\varepsilon}$ is chosen to be approximately 0.6. However, since in our discrete formulas the lag is an integer multiple of $\Delta$, i.e. $u^{\varepsilon}=r\times\Delta(\varepsilon)$ the lag changes slightly for different values of $\varepsilon$. The values of the lag for simulations with different values of $\varepsilon$ are chosen to be

[TABLE]

Numerical estimates for the $L^{2}$-errors of parameter estimators are presented in Figure 3. The $L^{2}$ error $\|\hat{\theta}^{\varepsilon}-\theta\|_{2}$ is depicted in the upper-left part of Figure 3. Since $\hat{\theta}^{\varepsilon}$ estimates the empirical mean of the volatility process, expression (61) is directly applicable in this case. Figure 3 demonstrates that, although the sub-sampling regime (79) is identical in both cases, the number of points for computing the realized volatility, $J$, significantly affects the behavior of parameter estimators. First, the numerical error is reduced significantly (approximately 10 times) for $J=\varepsilon^{-2}$ compared to $J=\varepsilon^{-1}$. Second, the asymptotic behavior for the $L^{2}$ error seems also to be affected by the choice of $J$ which is most evident for parameter $\theta$. For the sub-sampling regime (79) and (77) the decay of $L^{2}$ error is much faster than $\varepsilon^{1/2}$ for all three parameters. However, with the choice of $J$ in (78) errors in parameter estimators are almost the same as for the estimators computed under direct observability and the error is proportional to $\varepsilon^{1/2}$. We would like to point out that numerical simulations presented here are for finite values of $\varepsilon\in[0.01,\ldots,0.1]$. We conjecture that for smaller values of $\varepsilon<0.01$ the convergence rate of all parameter estimators computed with $J=\varepsilon^{-1}$ should change to $\varepsilon^{-1/2}$ and asymptote to the black line corresponding to the estimators computed under direct observability.

Our numerical simulations have important practical consequences. In particular, our numerical results suggest that it is important to follow the regime $J=1/\varepsilon^{2}$ for larger values of $\varepsilon$. However, one can switch to a different regime (e.g. $J\sim\varepsilon^{-3/2}$ or even $J\sim\varepsilon^{-1}$) for smaller values of $\varepsilon$ to reduce the computational overhead. This is motivated by rather fast rate of convergence for parameter estimators computed with $J=\varepsilon^{-1}$.

In the regime with $J=\varepsilon^{-1}$ for smaller values of $\varepsilon=0.01,0.02$ errors in all parameter estimators decay significantly compared to $\varepsilon=0.05,0.1$. Behavior of parameter estimators themselves is depicted in Figure 4. It is obvious that the sub-sampling regime $J=\varepsilon^{-1}$ results in very large errors for larger values of $\varepsilon=0.05$, $0.1$. Parameter $\theta$ is estimated more accurately under the computational scheme with $J=\varepsilon^{-1}$ for $\varepsilon=0.01$, but there is still approximately 10% relative error in estimating parameters $\kappa$ and $\gamma$ in this regime. On the other hand, relative errors in estimating all three parameters are much smaller for $J=\varepsilon^{-2}$ and $\varepsilon=0.1$. Therefore, the most beneficial strategy is to use a bigger window, $\varepsilon$, for computing the realized volatility with a large number of points $J=\varepsilon^{-2}$ for the return process.

Asymptotic behaviour of our observable estimators for the Heston parameters $\kappa$ and $\gamma$ strongly depends on the behavior of the lagged covariances $K^{\varepsilon}(u)$ of $Y_{t}^{\varepsilon}$. Thus, we also present our numerical results for the estimation of $\hat{K}^{\varepsilon}(u)$ with two particular time lags $u=0$ and $u\approx 0.6$.

The mean and lagged covariances of $Y_{t}^{\varepsilon}$ are approximated by their empirical estimators, given by

[TABLE]

where for lag $u=0$ the integer $s$ is $s=0$, and is chosen such that $s\Delta\approx 0.6$ when lag $u\approx 0.6$. Recall that the stationary moments of true volatilities are given by (49), (50).

The $L^{2}$ errors for the lagged covariances are computed from Monte-Carlo simulations as

[TABLE]

where the sum involves $MC=1000$ independent evaluations of $\hat{K}^{\varepsilon}(u)$. Results for the covariance estimation are displayed in Figure 5.

Behavior of $L^{2}$ errors for estimated second moments is consistent with the behavior of parameter estimators discussed earlier. In particular, for the range of $\varepsilon\in[0.01,\ldots,0.1]$ convergence rate of $\hat{K}^{\varepsilon}(u)$ computed with $J=\varepsilon^{-1}$ appears to be much faster than $\varepsilon^{-1/2}$, especially for $\hat{K}(0)$. Similar to the behavior of parameter estimators, we conjecture that this is due to the finite range of $\varepsilon$.

The choice of the lag $u^{\varepsilon}$ is motivated by some practical considerations. In particular, one should perform an a-posteriori check after the parameter estimator $\hat{\kappa}$ is computed and ensure that the estimated lagged correlation $K(u^{\varepsilon})$ is not too close to 0 or 1, for instance by checking that $e^{-\hat{\kappa}u^{\varepsilon}}$ lies between $0.3$ and $0.7$. Apart from such practical constraint above, the choice of $u^{\varepsilon}$ is otherwise arbitrary. We performed numerical simulations (not reported here) investigating several other choices of the lag $u^{\varepsilon}$. In particular, we considered $u^{\varepsilon}\approx 0.3$ and the “vanishing lag” case $u^{\varepsilon}=\Delta\equiv\sqrt{\varepsilon}$. Our numerical simulations indicate that for the specific Heston SDE parameters considered here the choice $u^{\varepsilon}\approx 0.6$ yielded near-optimal asymptotic behavior of both, observable moments estimators and parameter estimators.

11 Conclusions

We carried out an extensive analytical and numerical investigation of the Heston joint SDEs driving jointly the squared volatilities $V_{t}$ and the rate of returns $R_{t}$. Since the volatility process $V_{t}$ cannot be observed directly, realized volatilities $Y_{t}^{\varepsilon}$ computed from the return process $R_{s}$ with $s$ in the sliding window $\left[t-\varepsilon,t\right]$ provide classical observable approximations of the unobservable $V_{t}$.

The main goal of this paper is to define and study observable estimators of the Heston SDEs parameters computed from the $Y_{t}^{\varepsilon}$, and exhibiting asymptotic consistency as $\varepsilon\to 0$. This context fits our general framework of indirect observability where parameter estimators for the dynamics of an unobservable process $X_{t}$ can only be computed from observations of a process $Y_{t}^{\varepsilon}$ approximating $X_{t}$ as $\varepsilon\to 0$. Computing realized volatilities $Y_{t}^{\varepsilon}$ from the rates of returns $R_{t}$ requires partitioning the window $\left[t-\varepsilon,t\right]$ into $J(\varepsilon)$ time intervals. For the Heston SDEs we prove precise bounds for $L^{q}$ norms $\|Y_{t}^{\varepsilon}-V_{t}\|_{q}$ in terms of $J(\varepsilon)$ and $\varepsilon$. In particular we show that $\|Y_{t}^{\varepsilon}-V_{t}\|_{4}\leq C\sqrt{\varepsilon}$ provided $J(\varepsilon)\sim 1/\varepsilon^{2}$. However, for small window sizes, $\varepsilon$, partition sizes $J(\varepsilon)\sim 1/\varepsilon^{2}$ are not very practical since they require an overwhelming number of points for small window size $\varepsilon$. Our numerical simulations indicate that it is possible to obtain reasonable numerical estimates in $L^{2}$ sense with more practical partition sizes $J(\varepsilon)\sim 1/\varepsilon$. However, $L^{4}$ errors $\|Y_{t}^{\varepsilon}-V_{t}\|_{4}$ are more sensitive to the choice of the partition size.

Our observable estimators of the Heston SDEs parameters are defined as explicit functions of the empirical mean and two empirical lagged covariances computed from $N(\varepsilon)$ observations $Y_{j\Delta(\varepsilon)}^{\varepsilon}$, $j=1,\ldots,N(\varepsilon)$ of the realized volatility, sub-sampled with time step $\Delta(\varepsilon)$. We prove that for fastest convergence speed of observable parameter estimators to true parameters, the optimal sub-sampling regime is provided by $N(\varepsilon)\sim\varepsilon^{-3/2}$ and $\Delta(\varepsilon)\sim\varepsilon^{1/2}$ with $J(\varepsilon)\sim\varepsilon^{-2}$. Our sub-sampling scheme (69) provides a needed balance between the $L^{2}$ errors of estimation on empirical covariances and the $L^{2}$ difference between true and realized volatilities. This optimal sub-sampling scheme corresponds to a total observational time $T(\varepsilon)=N(\varepsilon)\Delta(\varepsilon)\sim 1/\varepsilon$ and a total number $n(\varepsilon)$ of observed returns rate values $n(\varepsilon)=N(\varepsilon)J(\varepsilon)=1/\varepsilon^{7/2}$.

Surprisingly, our numerical simulations indicate a much faster speed of convergence for the sub-sampling regime $N(\varepsilon)\sim\varepsilon^{-3/2}$ and $\Delta(\varepsilon)\sim\varepsilon^{1/2}$ with $J(\varepsilon)\sim\varepsilon^{-1}$ for a rather wide range of $\varepsilon\in[0.01,0.1]$. Convergence rates for all three parameter estimators are close to $\varepsilon^{1}$. Parameter estimators computed under indirect observability are inferior to estimators computed from directly observed time-series of the volatility process. Therefore, we conjecture that the convergence rate for the sub-sampling regime with $J=1/\varepsilon$ should change for $\varepsilon\ll 0.01$. However, verifying this with numerical simulations is extremely computationally costly. In addition, we also observe that relative errors under the sub-sampling regime $J=1/\varepsilon$ are much larger compared to the sub-sampling regime with $J=1/\varepsilon^{2}$. Therefore, to reduce the computational cost, the optimal estimation strategy is to use a larger window $\varepsilon$ for computing the realized volatility with a large number of points $J=1/\varepsilon^{2}$ for the return process.

When one imposes a bound on the total observational time $T(\varepsilon)$, our theory and numerical simulations indicate that there is a lower $L^{2}$ bound on the estimation errors for the parameters of the Heston volatility SDE. An upper bound on $T(\varepsilon)$ essentially forces a lower-bound on $\varepsilon$. Therefore, in practice, it is then not beneficial to keep over-refining the partition of the sliding time window used to compute realized volatilities. Moreover, when $T(\varepsilon)$ is bounded, decreasing the size of the sliding window, $\varepsilon$, constrains the number of observations of the return process inside this window to decrease, and this generates more inaccurate approximations of true volatilities by realized volatilities.

Our theoretical analysis and numerical simulations of the Heston SDEs presented here provide practical guidelines for fitting joint Heston SDEs to practical observations of stock prices. In particular our results should help define adequate choices for the size $\varepsilon$ of the sliding windows used to compute realized volatilities, as well as for the selection of an efficient sub-sampling time step of returns rate observations.

Acknowledgements. I.T. and R.A. were supported in part by the NSF Grant DMS-1109582. I.T. is also partially supported by the NSF Grant DMS-1620278.

Appendix A Polynomial functions of volatilities and Theorem 2

We evaluate conditional moments for polynomial functions of squared volatilities $V_{t}$. Let ${\cal F}_{s}$ be the filtration generated by the Brownian $B_{t}$ driving the Heston volatility SDE. Note that conditioning by ${\cal F}_{s}$ gives the same results when the volatility process starts at any fixed $V_{0}=y>0$ or when it is the only stationary process driven by the volatility Heston SDE.

Recall the statement of Theorem 2. Fix any polynomial $h$ of total degree $n$ in $k$ variables $(x_{1},\ldots,x_{k})$. Let $0=u(0)<u(1)<\ldots<u(k)$ be any sequence of $k+1$ lag instants. For $T>0$, define random variables $H$ and $H_{T}$ by

[TABLE]

Recall that $\nu_{T}=e^{-T\kappa}$. Define $w_{j}=e^{-\kappa(u(j+1)-u(j))}$ for $j=0,\ldots,k-1$. There is then a polynomial $POL$ in $k+2$ variables such that for all $T>0$ and all $y>0$

[TABLE]

The degree and coefficients of POL are determined by the integers $n,k$, the coefficients of $h$, and the vector $\boldsymbol{\theta}$. The asymptotic polynomial moments are then given by

[TABLE]

For any integer $q\geq 1$ there is a positive constante $C$, and an integer $p\geq 1$, determined only by $q,k,\boldsymbol{\theta}$ and the polynomial $h$ such that, for all positive $T$ and $y$, and all $0=u(0)<u(1)<\ldots<u(k)$

[TABLE]

In particular for $q=1$ one has

[TABLE]

Remarks. Equation (84) also implies that as $T\to\infty$, the random polynomial functionals $H_{T}$ converge in $L^{q}$-norm to the constants ${\mathbb{E}}_{\psi}(H)$, where $L^{q}$-norms are computed under ${\mathbb{E}}_{y}$. Note also, that all the constants introduced in the theorem and in its proof below do not depend on the actual lags $u(0)<u(1)<\ldots<u(k)$.

Proof of Theorem 2:

Proof.

By linearity, we only need to consider the case when $h$ is a monomial in $k$ variables. For $k=1$, the result was proved by (44). Proceeding by recurrence on $k$, assume the result is true for monomials in $k-1$ variables $(x_{2},\ldots,x_{k})$. Any monomial $h$ in $k$ variables can be written as $h=x_{1}^{m}g(x_{2},\ldots,x_{k})$. Define

[TABLE]

The recurrence hypothesis provides a polynomial $R$ in $(k+1)$ variables such that, for all $T$

[TABLE]

where the coefficients of $R$ are determined by $g,\boldsymbol{\theta}$. By the Markov property we thus get

[TABLE]

Since ${\mathbb{E}}_{y}[H_{T}]={\mathbb{E}}_{y}[V_{u(1)+T}^{m}{\mathbb{E}}(G_{T}|{\cal F}_{u(1)+T})]$ we then obtain

[TABLE]

Each monomial $M$ of $R$ is of the form $\nu_{T}^{p}(V_{u(1)+T}\nu_{T})^{j}S(w_{1},w_{2},\ldots,w_{k-1})$ for some $p$, $j$ and some polynomial $S$. Then in the right-hand side of (83), $M$ contributes a term of the form

[TABLE]

Due to (44) with $q=m+j$, this last conditional expectation is a polynomial in the two variables

[TABLE]

with coefficients depending only on $m+j$ and $\boldsymbol{\theta}$. Hence $\Gamma(M)$ is a polynomial in $\nu_{T}$ and $y\nu_{T}$, with coefficients which are polynomials in $(w_{0},w_{1},w_{2},\ldots,w_{k-1})$, fully determined by $m$, $j$, $\boldsymbol{\theta}$. The same property must then hold for the sum ${\mathbb{E}}_{y}(H_{T})$ of all the $\Gamma(M)$ contributed by the monomials $M$ of $R$. This completes the proof of (83) by recurrence on $k$.

Write $POL$ in (83) as a polynomial $POL(z)$ in the $k+2$ variables $z_{i}$. The vector $z(T)=(\nu_{T},y\nu_{T},w_{0},\ldots,w_{k-1})$ tends to $z(\infty)=(0,0,w_{0},\ldots,w_{k-1})$ as $T\to\infty$. The polynomial $Q(z(T))=POL(z(T))-POL(z(\infty))$ can be written for some integer $p$

[TABLE]

where for $s=0,\ldots,p$, each $A_{s}$ is a polynomial in the $(k+1)$ variables $(\nu_{T},w_{0},\ldots,w_{k-1})$. Since all these positive $(k+1)$ variables are inferior to $1$, then each $|A_{s}|$ remains bounded for all $T\geq 0$ and all $u(0)<u(1)<\ldots<u(k)$. Hence there is a constant $C$ such that

[TABLE]

For all $s\geq 1$ we have $\nu_{T}^{s}\leq\nu_{T}=e^{-T\kappa}$, and hence the expansion of $Q(z(T))$ provides a new constant $C_{1}$ such that, for all $u(0)<u(1)<\ldots<u(k)$,

[TABLE]

This proves (85).

Let ${\overline{H}}={\mathbb{E}}_{\psi}(H)$. Expand $\beta(T)=(H_{T}-{\overline{H}})^{q}$ as a linear combination of terms of the form ${\overline{H}}^{q-j}H_{T}^{j}$ for $j=0,\ldots,q$. Recall that $h$ is a polynomial in $x_{1},x_{2},\ldots,x_{k}$. For $j$ fixed, $\sigma_{j}=h^{j}$ is also a polynomial in $x_{1},x_{2},\ldots,x_{k}$. By definition (82), we can express both $\Sigma=H^{j}$ and $\Sigma_{T}=H_{T}^{j}$ as

[TABLE]

For each $j$, equation (85) applied to the polynomial $\sigma=h^{j}$ provides a constant $C_{j}$ and an integer $p(j)$ such that

[TABLE]

and hence there are constants $c_{j}$ such that

[TABLE]

Applying this to $j=0,\ldots,q$ and using the Newton binomial formula yields, for some new constant $C$,

[TABLE]

which completes the proof of (84).