Estimation of Tempered Stable Lévy Models of Infinite Variation

Estimation of Tempered Stable L\'evy Models of Infinite Variation

Introduction

L\'{e}vy processes have experienced a revival in the past $20$ years, propelled by the need for more realistic modeling of irregular behavior in many phenomena of nature and society. These fundamental building blocks of stochastic modeling have been widely applied in many fields, including statistical physics, meteorology, seismology, insurance, finance, and telecommunication. While, in principle, L\'{e}vy models offer ideal conditions for estimation purposes, two main bottlenecks complicate their estimation. Firstly, their marginal distributions often lack tractable or closed-form representations. In those situations, the marginal distributions must be approximated by Fourier, Monte Carlo, or other numerical methods, which makes the estimation slower and noisier. The second issue comes from the need to handle high-frequency sampling data of the process. This type of data has been widely available in finance during the last $15$ years and is increasingly {more common} in other fields. The two just-mentioned issues have rendered traditional statistical methods such as likelihood and Bayesian estimation unfeasible. {\color{Blue} We refer the reader to ContTankov:2004 for more information about L\'evy processes and their application in finance, Masuda for a survey on frequentist parametric estimation of L\'evy process, and FigueroaQiKuffner for more information about Bayesian estimation methods.}

In this paper, we study {a new method} for the estimation of the parameters of a L\'{e}vy model. A semiparametric model is considered in which the jump component is assumed to exhibit small jumps that behave like those of {\color{Blue} a $Y$-stable} L\'{e}vy process. Specifically, the class of tempered stable processes introduced in FigueroaLopezGongHoudre:2016\footnote{The term “tempered stable" is understood here in {\color{Red} a more} general sense than in several classical sources of financial mathematics (e.g., Applebaum:2004, ContTankov:2004, KyprianouSchoutensWilmott:2005) and even more general than in Rosinski:2007. In fact, such class of L\'{e}vy processes is called {\color{Red} the} tempered-stable-like L\'{e}vy processes in FigueroaLopezGongHoudre:2016.} and FigueroaLopezOlafsson:2019(1) is considered. We focus on models of infinite variation (i.e., {$Y\in(1,2)$}), which are arguably the most relevant for financial applications (see AitSahaliaJacod:2009, Belomestny:2010, and FigueroaLopez:2012). The estimation of semiparametric L\'{e}vy models of infinite jump variation under high-frequency data is not well developed. Jacod and Todorov JacodTodorov:2014 were the first to introduce an efficient estimator of the integrated volatility of an It\^{o} semimartingale model in the presence of a L\'{e}vy jump model of infinite variation with Blumenthal-Getoor index {$\beta\in(1,3/2)$} or when the jump component is symmetric. {\color{Blue} Their estimator is based on locally estimating the volatility from the empirical characteristic function of the increments of the process over time blocks of decreasing length.} Recently, Mies Mies:2019 proposed an efficient estimation method for L\'{e}vy models based on a type of approximate semiparametric method of moments with scaling. Specifically, for some suitable moment functions $f_{1},f_{2},\dots,f_{m}$ and a scaling factor $u_{n}\rightarrow\infty$, Mies:2019 {proposed to look} for the parameters $\widehat{\boldsymbol{\theta}}=(\widehat{\theta}_{1},\ldots,\widehat{\theta}_{m})$ such that

where $\widetilde{Z}$ is the superposition of a Brownian motion and independent stable L\'{e}vy processes closely approximating $X$ in a certain sense. The distribution measure $\mathbb{P}_{\boldsymbol{\theta}}$ of $\widetilde{Z}$ depends on some parameters $\boldsymbol{\theta}$, including the volatility $\sigma$ of $X$, and $\mathbb{E}_{\boldsymbol{\theta}}(\cdot)$ denotes the expectation with respect to $\mathbb{P}_{\boldsymbol{\theta}}$. Above, $\Delta_{i}^{n}L:=L_{t_{i}}-L_{t_{i-1}}$ is the $i$-th increment of a generic process $(L_{t})_{t\geq 0}$ given $n$ evenly spaced random samples $L_{t_{0}},\ldots,L_{t_{n}}$ over a fixed time interval $[0,T]$ (i.e., $t_{i}=ih_{n}$ with $h_{n}=T/n$). If $X$ were assumed to {\color{Red} follow} a parametric L\'{e}vy model and we replaced $\mathbb{E}_{\widehat{\boldsymbol{\theta}}}(f_{j}(u_{n}\Delta^{n}_{i}\widetilde{Z}))$ with $\mathbb{E}_{\widehat{\boldsymbol{\theta}}}(f_{j}(u_{n}\Delta^{n}_{i}X))$ in (ref), we will recover a standard Method of Moment Estimator (MME). However, we are assuming that $X$ is semiparametric and that it can be approximated closely enough by {a} parametric L\'{e}vy model $\widetilde{Z}$. The scaling $u_{n}$, which is taken to converge to $\infty$ at the order of $1/\sqrt{\ln(n)/n}$, is also a new feature of this method compare to the standard MME.

The moment functions $f_{1},\ldots,f_{m}$ and the scaling factor $u_{n}$ in (ref) critically affect the performance of the estimators. To determine an appropriate scaling $u_{n}$, we connect it to the threshold parameter $\varepsilon_{n}$ of a Truncated Realized Quadratic Variation (TRQV),

which is known to be a consistent estimator for the integrated volatility of a general semimartingale model. {\color{Red} Again, above $\Delta_i^nX=X_{t_i}-X_{t_{i-1}}$ and we are assuming regular sampling observations $X_{t_1},\dots,X_{t_n}$ with $t_i=ih_n$ and $h_n=T/n$. Next, note that} by taking $f_{1}(x)=x^{2}{\bf 1}_{\{|x|\leq 1\}}$ in (ref), we {\color{Red} recover the TRQV (ref), which suggests the relationship} $u_{n}=1/\varepsilon_{n}$. That is, $1/u_{n}$ plays the same role as the threshold in TRQV.

Recently, FigueroaLopezMancini:2019 studied the problem of optimal thresholding of TRQV {\color{Red} (ref)} under the mean-square error. Specifically, in the case of a L\'{e}vy process with volatility $\sigma$, it is shown that the threshold $\varepsilon=\varepsilon^{\star}_{n}$ that minimizes the mean-square error, $\mathbb{E}((\text{TRQV}_{n}(\varepsilon)-\sigma^{2}T)^{2})$, solves the equation:

where {\color{Red} $b_{1,h_{n}}(\varepsilon):=X_{h_n}^{2}{\bf 1}_{\{|X_{h_n}|\leq\varepsilon\}}$}. By analyzing the small-time asymptotic behavior of $\mathbb{E}(b_{1,h_{n}}(\varepsilon))$ (i.e., when $n\rightarrow\infty$ so that $h_{n}\rightarrow 0$), FigueroaLopezMancini:2019 {\color{Red} proved} that the optimal threshold $\varepsilon^{\star}_{n}$ for a L\'{e}vy process with a $Y$-stable jump component behaves like

where {\color{Red} $h_{n}=T/n$} is the time span between observations and, as usual, $a_{n}\sim b_{n}$ means $a_{n}/b_{n}\rightarrow 1$ as $n\rightarrow\infty$. The proportionality constant $\sqrt{2-Y}$ roughly tells us that the higher the jump activity is, the lower the optimal threshold has to be if we want to discard the higher noise represented by the small jumps. This fact opens the door to an iterative method to estimate $\sigma^{2}$. We can first estimate $Y$ and $\sigma^{2}$ using, for instance, the method of moments (ref). We can then use the TRQV with the threshold $\widehat{\varepsilon}^{\,\star}_{n}=\sqrt{(2-\widehat{Y})\widehat{\sigma}^{2}h_{n}\ln(1/h_{n})}$.

In this paper, we first extend the result of FigueroaLopezMancini:2019 to allow for a general tempered stable L\'{e}vy process. Furthermore, we {propose a new approximation for} $\varepsilon^{\star}_{n}$ of the form:

where $C$ controls the overall intensity of jumps. The approximation (ref) says that if $C$ is small (relative to $\sigma$) then the threshold can be loosened up (in fact, $\widetilde{\varepsilon}_{n}^{\,\star}\nearrow\infty$ as $C\searrow 0$ as it should be). In practice $C$ is small compare to $\sigma$ and (ref) provides a significant correction compare to (ref). We then proceed to {devise} a new method to estimate the volatility, the index of jump activity $Y$, and $C$ by combining a variation of the {approximate semiparametric} method of moments in Mies:2019, TRQVs, and the approximate optimal threshold (ref). Compared to Mies:2019 we introduce simpler moment functions $f_{1},\ldots,f_{m}$, and a systematic and objective method to tune the scaling factor $u_{n}$ in (ref). The performance of the proposed procedure is superior to the efficient methods of JacodTodorov:2014 and Mies:2019. Finally, as in JacodTodorov:2014, we use a localization technique to estimate the integrated volatility of an It\^{o} semimartingale. Specifically, the idea is to split the time horizon into small blocks where the process is approximately L\'{e}vy and, hence, its volatility level can be estimated using our method. For values of $Y\geq 1.5$, our method outperforms the method proposed by JacodTodorov:2014.

The rest of this paper is organized as follows. Section (ref) provides the framework and assumptions as well as some known preliminary results from the literature. Section (ref) obtains the asymptotic behavior of $\mathbb{E}(b_{1,h_{n}}(\varepsilon))$ and derives (ref). The second-order approximation (ref) is derived in Section (ref) as well as a numerical assessment of the approximations in the case of a CGMY jump component. The new method to estimate the parameters of a tempered stable L\'{e}vy model is presented in Section (ref) together with an analysis of its performance via Monte Carlo simulations. The proofs are deferred to an appendix section.

The Model and Some Preliminary Results

Throughout, $\mathbb{R}_{+}:=[0,\infty)$ and $\mathbb{R}_{0}:=\mathbb{R}\backslash\{0\}$, and we let $(\Omega,\mathscr{F},\mathbb{F},\mathbb{P})$ be a complete filtered probability space on which all stochastic processes are defined, where $\mathbb{F}:=(\mathscr{F}_{t})_{t\in\mathbb{R}_{+}}$ satisfies the usual conditions. We consider a L\'{e}vy process $X:=(X_{t})_{t\in\mathbb{R}_{+}}$ of the form

where $W:=(W_{t})_{t\in\mathbb{R}_{+}}$ is a Wiener process and $J:=(J_{t})_{t\in\mathbb{R}_{+}}$ is an independent pure-jump {tempered stable} L\'{e}vy process with L\'{e}vy triplet $(b,0,\nu)$. The L\'{e}vy measure $\nu$ is assumed to be absolutely continuous with a density $s:\mathbb{R}_{0}\rightarrow\mathbb{R}_{+}$ of the form

Here, $C_{\pm}>0$, $Y\in(1,2)$, and $q:\mathbb{R}_{0}\rightarrow\mathbb{R}_{+}$ is a bounded Borel-measurable function. Concretely, we make the following assumptions on $q$.

Using a density transformation technique in Sato:1999, we can change the probability measure from $\mathbb{P}$ to another locally absolutely continuous measure $\widetilde{\mathbb{P}}$, under which $J$ is a $Y$-stable L\'{e}vy process and $W$ is a standard Brownian motion independent of $J$. Concretely, let

Note that $\widetilde{\nu}$ is the L\'{e}vy measure of a $Y$-stable L\'{e}vy process and, also,

Next, define $\widetilde{\mathbb{P}}$ such that, for any $t\in\mathbb{R}_{+}$,

By virtue of Sato:1999, a necessary and sufficient condition for the measure transformation from $\mathbb{P}$ to $\widetilde{\mathbb{P}}$ to be well defined is given by

which can be shown to follow from Assumption (ref)$-$(i) & (ii) (cf. FigueroaLopezOlafsson:2019(2)). Under $\widetilde{\mathbb{P}}$, $J$ is a L\'{e}vy process with L\'{e}vy triplet $(\widetilde{b},0,\widetilde{\nu})$, and $W$ is a standard Brownian motion which is independent of $J$. In particular, under $\widetilde{\mathbb{P}}$, the centered process $Z:=(Z_{t})_{t\in\mathbb{R}_{+}}$, given by

is a strictly $Y$-stable process with its skewness, scale, and location parameters given by $(C_{+}-C_{-})/(C_{+}+C_{-})$, $\left\{(C_{+}+C_{-})\Gamma(-Y)|\cos(\pi Y/2)|\right\}^{1/Y}$, and $0$, respectively. Let $p_{Z}$ denote the marginal density of $Z_{1}$ under $\widetilde{\mathbb{P}}$. It is well known (cf. Sato:1999 and references therein) that

so that

The processes $U:=(U_{t})_{t\in\mathbb{R}_{+}}$ and $Z$ can be expressed in terms of the jump-measure $N(dt,dx)$ of the process $J$ and its compensator $\widetilde{N}(dt,dx):=N(dt,dx)-\widetilde{\nu}(dx)dt$ (under $\widetilde{\mathbb{P}}$), as follows:

where

The existence of the integral defining $\eta$ follows from Assumption (ref)$-$(i) & (ii). Clearly, $Z^{+}:=(Z^{+}_{t})_{t\in\mathbb{R}_{+}}$ and $-Z^{-}:=(-Z^{-}_{t})_{t\in\mathbb{R}_{+}}$ are independent one-sided $Y$-stable processes with scale, skewness, and location parameters given by $\left(C_{\pm}|\Gamma(-Y)\cos(\pi Y/2)|\right)^{1/Y}$, $1$, and $0$, respectively, so that

Moreover, it can be shown that (cf. FigueroaLopezGongHoudre:2017) there {\color{Blue} exists a universal constant} $K\in(0,\infty)$, such that for any $z>0$,

{\color{Red} Combining (ref) and (ref), we deduce that there exists a constant $\widetilde{K}\in(0,\infty)$ such that, for any $z>0$,

Furthermore, using (14.34) in Sato:1999 and an argument similar to that in the proof of FigueroaLopezGongHoudre:2017, we can show that:

where above, without loss of generality, we use the same constant $\widetilde{K}$ as in (ref).}

Main Result

The TRQV, defined as

is one of the most {commonly used estimators} for the integrated volatility of an It\^{o} semimartingale. Above, $\Delta_{i}^{n}X:=X_{t_{i}}-X_{t_{i-1}}$ for $i=1,\ldots,n$, where $X_{t_{0}},X_{t_{1}},\ldots,X_{t_{n}}$ are evenly spaced {observations} of $X$ over a fixed time horizon $[0,T]$, so that $t_{i}=t_{i,n}=ih_{n}$ for $i=0,1,\ldots,n$, with $h_{n}:=T/n$. One of its drawbacks is the necessity of tuning the threshold $\varepsilon$ up, which strongly affects the performance of the estimator. It is shown in FigueroaLopezMancini:2019 that, for a L\'{e}vy process $X$ with volatility $\sigma>0$, there exists a unique threshold $\varepsilon=\varepsilon^{\star}_{n}$, which minimizes the mean-square error, $\mathbb{E}((\widehat{\sigma}^{2}_{n}(\varepsilon)-\sigma^{2})^{2})$. {Furthermore, the minimizer $\varepsilon^{\star}_{n}$ is such that}

and {solves} the equation

where $b_{1,h_{n}}(\varepsilon):=X_{h_{n}}^{2}{\bf 1}_{\{|X_{h_{n}}|\leq\varepsilon\}}$. Therefore, in order to determine the asymptotic behavior of the optimal threshold $\varepsilon^{\star}_{n}$, we need to study the asymptotic behavior of $\mathbb{E}(b_{1,h}(\varepsilon))$ as both $h\rightarrow 0+$ and $\varepsilon=\varepsilon(h)\rightarrow 0+$ in such a way that $\varepsilon(h)/\sqrt{h}\to\infty$, as $h\rightarrow 0$. {Our main theoretical result} accomplishes this for the tempered stable L\'{e}vy processes of Section (ref), {and its} proof is deferred to Appendix (ref).

The following result gives the asymptotic behavior of the optimal threshold $\varepsilon_{n}^{\star}$. Its proof is similar to that of FigueroaLopezMancini:2019 and is outline below for completeness and also to motivate some approximation methods proposed below.

The proportionality constant $\sqrt{2-Y}$ of the previous result is intuitive and roughly tells us that the higher the jump activity is, the lower the optimal threshold has to be if we want to discard the higher noise represented by the jumps and to catch information about the Brownian component.

Other Approximations and Illustration for a CGMY Model

In this section, we introduce other approximations to the optimal threshold derived from the formulas in Theorem (ref) and the proof of Corollary (ref). We then illustrate their performance in the case of a L\'{e}vy process with a CGMY jump component $J$ (cf. CarrGemanMadanYor:2002). The CGMY model is considered a prototypical jump process of infinite activity in finance. In the notation of the L\'{e}vy density (ref), a CGMY model is given by

Thus, the conditions of Assumption (ref) are satisfied with $\alpha_{+}=-M$ and $\alpha_{-}=G$. We adopt the parameter setting

These values are similar to those used in FigueroaLopezOlafsson:2019(2)\footnote{FigueroaLopezOlafsson:2019(2) considers the asymmetric case $\nu(dx)=C(x/|x|)\bar{q}(x)|x|^{-1-Y}\,dx$ with $C(1)=0.015$ and $C(-1)=0.041$. Here, we take $C=(C(1)+C(-1))/2$ in order to simplify the simulation of the model. Our values of $G$, $M$, and $Y$ are the same as in FigueroaLopezOlafsson:2019(2).}, who themselves took them from an empirical study in Kawai:2010. We take $T=1$ year and $n=252(6.5)(60)$, which corresponds to a frequency of $1$ minute (assuming $252$ trading days and $6.5$ trading hours per day).

To compute $\mathbb{E}(b_{1,h}(\varepsilon))$, we use Monte Carlo and the change of probability measure (ref). Concretely, under $\widetilde{\mathbb{P}}$, we have the following representation:

where $Z_{h}^{+}$ and $-Z_{h}^{-}$ are independent one-sided $Y$-stable random variables with common scale, skewness, and location parameters given by $C|\Gamma(-Y)\cos(\pi Y/2)|h^{1/Y}$, $1$, and $0$, respectively. Such a distribution can be simulated efficiently\footnote{In our code, we use the R package {\tt stabledist} to generate them.}.

We consider two different approximations of the equation (ref) defining the optimal threshold $\varepsilon^{\star}_{n}$. For the first approximation, we replace {$\mathbb{E}(b_{1}(\varepsilon,h_{n}))$ (where $h_{n}=1/n$)} in (ref) with its leading order terms as given by Theorem (ref), namely,

For the second approximation, we take a simplified version of (ref), only keeping those terms that are found to be {significant:}

{Interestingly, as $C\rightarrow 0$, we} have $\widetilde{\varepsilon}_{n}^{\,\star}\rightarrow\infty$, which makes sense. The approximation (ref) says that if $C$ is small (relative to $\sigma$) then the threshold can be loosened up.

Figure (ref) shows the graphs of the left-hand expressions of (ref) (solid blue) and the approximation (ref) (dashed red) against $\varepsilon$ for three different values of $\sigma$: $0.1$, $0.2$, and $0.4$. The solid blue vertical line is the “true" optimum threshold $\varepsilon=\varepsilon^{\star}_{n}$, the dotted brown vertical line shows $\varepsilon=\widetilde{\varepsilon}_{n}^{\,\star}$ with $\widetilde{\varepsilon}_{n}^{\,\star}$ given as in approximation ((ref)), and the dotted/dashed vertical green line is the approximation $\varepsilon=\varepsilon_{n}:=\sqrt{(2-Y)\sigma^{2}h_{n}\ln(1/h_{n})}$ derived in {(ref) of} Corollary (ref). We also show the vertical line passing at the root of (ref) (vertical dashed red). It is evident that for the considered values of $Y$ and $\sigma$, the root of (ref) and $\widetilde{\varepsilon}_{n}^{\,\star}$ are reasonably good approximations of $\varepsilon^{\star}_{n}$. However, we cannot say the same about $\varepsilon_{n}=\sqrt{(2-Y)\sigma^{2}h_{n}\ln(1/h_{n})}$, which is a good approximation of $\varepsilon^{\star}_{n}$ only for small values of $\sigma$ and, otherwise, it underestimates $\varepsilon^{\star}_{n}$.

Next, we consider the value of $Y=1.5$, while all the other CGMY parameter values remain unchanged. Figure (ref) below shows the graphs of the left-hand expressions of (ref) (solid blue) and (ref) (dashed red), against $\varepsilon$ for three different values of $\sigma$: $0.1$, $0.2$, and $0.4$. The Equation (ref) derived from Theorem (ref) is a relatively accurate approximation of (ref), especially for larger values of $\sigma$. As before, the approximation (ref) established in Corollary (ref) is accurate for small and medium values of $\sigma$ but not for larger values. The approximation (ref) is reasonably accurate for all considered values of $\sigma$.

Finally, we consider the value of $Y=1.7$. All the other CGMY parameter values remain the same. The approximations are shown in Figure (ref). We deduce that for such a large value of $Y$, the approximation (ref) derived from {Theorem} (ref) is not accurate anymore, though it improves as $\sigma$ gets larger. On the other hand, {the other suggested approximation (ref) is still relatively} accurate to approximate the optimal threshold {$\varepsilon_{n}^{\star}$} (the root of (ref)). We again have that for small {and medium} values of $\sigma$, the approximation (ref) is good, which is not the case for large values of $\sigma$.

To summarize, while for values of $Y\leq 1.5$, the approximation ((ref)) may be the most accurate, this is not the case anymore for larger values of $Y$. On the other hand, the approximation (ref) is reasonably good for a large range of values of $Y$. Due to this reason, in our simulations of Section (ref), we use {(ref) to assess the finite sample performance of the proposed estimation method below.

A New Method To Estimate The Volatility

In this section, we propose a new method for estimating the volatility $\sigma^{2}$ and other parameters of a tempered-stable L\'{e}vy process using the TRQV (ref) and the approximations of the optimal threshold derived in Section (ref). Then, we illustrate the method in the case of a CGMY L\'{e}vy process. Finally, using a localization technique, we adapt our method to estimate the integrated variance under a Heston stochastic volatility model with CGMY jumps and compare it to the method proposed by JacodTodorov:2014, which is known to be efficient when $Y\leq 1.5$.

Estimation of Stable-Like L\'{e}vy Measures

As shown by (ref) and the asymptotic expansion of Theorem (ref), the optimal threshold $\varepsilon^{\star}_{n}$ depends on the volatility, and vice versa. It is then natural to consider an iterative method to estimate $\varepsilon^{\star}_{n}$. But before this, we need to estimate $C_{\pm}$ and $Y$. Several methods have been proposed in the literature for this purpose (see, e.g., AitSahaliaJacod:2009, Bull:2016, and Reiss:2013). Mies Mies:2019 recently proposed an efficient method using the method of moments. In this part, we adapt and modify this method and apply it in combination with the approximations of Theorem (ref) to estimate the optimal threshold $\varepsilon^{\star}_{n}$ of the TRQV and subsequently the other parameters $\sigma$, $Y$, and $C_{\pm}$.

Consider a L\'{e}vy process $X:=(X_{t})_{t\in\mathbb{R}_{+}}$ with characteristic triplet $(\mu,\sigma^{2},\nu)$. The approach of Mies:2019 builds on the assumption that $\nu$ can be well approximated by the superposition of stable L\'{e}vy measures in the sense that

for some $L,\rho\in(0,\infty)$, where $\widetilde{\nu}$ is given by

for some $N\in\mathbb{N}$, $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{N})\in(0,2)^{N}$, and $\boldsymbol{r}=(r_{1}^{+},r_{1}^{-},\ldots,r_{N}^{+},r_{N}^{-})\in\mathbb{R}_{+}^{2N}$ such that

for some $\alpha_{0}\in(0,\infty)$. We want to estimate $\boldsymbol{\theta}:=(\sigma^{2},\boldsymbol{r},\boldsymbol{\alpha})$ given $n$ observations, $X_{t_{1}},X_{t_{2}},\ldots,X_{t_{n}}$, of the process $X$ at known times $0=t_{0}<t_{1}<\cdots<t_{n}=T$. As before, we assume the sampling times are evenly spaced and we done the time step between observations as $h_{n}:=T/n$. Conditions (ref), (ref), and (ref) essentially say that we can approximate $X$ by a fully specified L\'{e}vy process $\widetilde{Z}:=(\widetilde{Z}_{t})_{t\in\mathbb{R}_{+}}$ with characteristic triplet $(0,\sigma^{2},\widetilde{\nu})$ and, hence, with the decomposition

where $W:=(W_{t})_{t\in\mathbb{R}_{+}}$ is a standard Brownian motion and $S^{m}:=(S_{t}^{m})_{t\in\mathbb{R}_{+}}$, $m=1,\ldots,N$, are independent $\alpha_{m}$-stable processes, independent of $W$, each with L\'{e}vy density $\alpha_{m}|z|^{-1-\alpha_{m}}(r_{m}^{+}{\bf 1}_{\{x>0\}}+r_{m}^{-}{\bf 1}_{\{x<0\}})$, respectively.

Mies Mies:2019 proposed to estimate the parameters, $\boldsymbol{\theta}=(\sigma^{2},\boldsymbol{r},\boldsymbol{\alpha})$, of the approximating process $\widetilde{Z}$ using the method of moments. We now proceed to briefly review her method. The first step is to choose $3N+1$ moment functions $\boldsymbol{f}=(f_{1},\ldots,f_{3N+1})^{\text{T}}$, one for each parameters of $\widetilde{Z}$, and a suitable scaling factor $u_{n}\propto 1/\sqrt{h_{n}\ln(1/h_{n})}$, where “$\propto$" hereafter means “proportional to". Next, define the MME $\widehat{\boldsymbol{\theta}}_{n}$ to be a solution of the following equation

where $\boldsymbol{0}=(0,\ldots,0)^{\text{T}}\in\mathbb{R}^{3N+1}$ and $\mathbb{E}_{\boldsymbol{\theta}}(\boldsymbol{f}(u_{n}\widetilde{Z}_{h_{n}}))$ denotes the expectation such that $\widetilde{Z}_{h_{n}}$ is determined by the parameter vector $\boldsymbol{\theta}$. Since $\widetilde{Z}_{h_{n}}$ is fully specified and, thus, its characteristic function is available, $\mathbb{E}_{\boldsymbol{\theta}}(\boldsymbol{f}(u_{n}\widetilde{Z}_{h_{n}}))$ can be computed by, e.g., Fourier methods.

Our idea is to combine a version of Mies' method with our results in Section (ref) to improve our estimation of $\sigma^{2}$ and $Y$. Concretely, we propose to first find the roots of $\boldsymbol{F}_{n}(\boldsymbol{\theta})$ and plug them into a suitable approximation of Equation (ref) to obtain an estimate of the optimal threshold $\varepsilon^{\star}_{n}$. This can in turn be used to estimate the volatility via thresholding. To solve the $3N+1$ equations (ref), we propose to solve an optimization problem with objective function of the form

where

This particular weights are motivated by the scaling of the Central Limit Theorem for $\boldsymbol{F}_{n}(\boldsymbol{\theta})$ established in Mies:2019.

For simplicity, suppose we only want to estimate $\alpha=\alpha_{1}$, $r_{1}^{\pm}$, and $\sigma^{2}$ (the method can easily be adapted to estimate more parameters of $\widetilde{\nu}$). We then {\color{Blue} propose} the following {procedure}:

Estimation of a L\'{e}vy Process with CGMY Jump Component

In this subsection, we apply the method introduced in the previous subsection to the case of a CGMY jump component and compare it to the estimators of Mies Mies:2019 and Jacod and Todorov JacodTodorov:2014. Specifically, we work with simulated data from the model (ref) where $J$, is a pure-jump CGMY L\'{e}vy process, independent of the Brownian motion $W$, with L\'{e}vy measure

We use the same values of $C$, $G$, and $M$ as in (ref), but with different values of $\sigma$ and $Y$. We consider observations of a $5$ minutes frequency over a one-year ($252$ days) time horizon with a trading time of $6.5$ hours per day (so that {\color{Blue} $n=252\times 6.5\times 12=19656$}). It should be clarified that we are indeed in the same setting as that of Subsection (ref) since, as $x\rightarrow 0$, $\nu_{CGMY}(x)=C|x|^{-1-Y}+O(|x|^{-Y})$. This suggests us to take $N=1$ in (ref) and to use a $Y$-stable process to approximate the CGMY process because only the parameters $\sigma^{2}$, $C$, and $Y$ are of primary interest. Then Assumptions (ref)$-$(ref) are satisfied with $\rho=Y$ and $\widetilde{Z}_{t}=\sigma W_{t}+S_{t}$, $t\in\mathbb{R}_{+}$, where $(S_{t})_{t\in\mathbb{R}_{+}}$ is a $Y$-stable process with L\'{e}vy measure $\widetilde{\nu}(dx):=C|x|^{-1-Y}dx$. The parameters of the approximating model are $\boldsymbol{\theta}=(\sigma^{2},C,Y)$.

Next, we choose the $3$ moment functions $\boldsymbol{f}=(f_{1},f_{2},f_{3})^{\text{T}}$ as

and a suitable scaling factor $u_{n}$ to be specified below. These functions are simpler than the ones proposed in Mies:2019 and were chosen because of their superior performance. Even though the moment functions (ref) do not meet the strict constraints imposed in Mies:2019 (see Assumptions (F1)$-$(F2) therein), we believe that most of the assumptions therein are not needed for the validity of the asymptotic theory in Mies:2019. This will be investigated in a future work together with an objective and systematic method to calibrate the moment functions.

To determine a suitable scaling factor $u_{n}$, we will connect it to the threshold parameter $\varepsilon$ of the TRQV estimator (ref). The key observation is to analyze the moment equation corresponding to the function $f_{3}$, namely,

which, after some trivial simplifications, can be written as

This suggests that $1/u_{n}$ has a similar role to that of the threshold $\varepsilon$ in the TRQV estimator; namely, the choice of $u_{n}$ should ensure that $\sum_{i=1}^{n}f_{3}(u_{n}\Delta_{i}^{n}X)$ is dominated by the Brownian component or, equivalently, to eliminate the increments in which the jump component $J$ of $X$ dominates the Brownian component. Hence, in what follows, we will fix $u_{n}$ as $1/\varepsilon_{n}$, where $\varepsilon_{n}=\sqrt{2{\sigma}_{0}^{2}h_{n}\ln(1/h_{n})}$ and $\sigma_{0}^{2}$ is a suitable initial estimate of $\sigma^{2}$. We consider the following initial values for $\sigma^{2}$:

where

Broadly, we recommend to use the loose estimator $\widehat{\sigma}_{n,01}^{2}$ as our initial value $\sigma_{0}^{2}$ if the volatility is “large" (say, 0.4 or larger), and, otherwise, use the tighter estimator $\widehat{\sigma}_{n,02}^{2}$.

For the moment functions $\boldsymbol{g}:=(g_{1},g_{2})^{\text{T}}$ in Step (ref) of the algorithm in Subsection (ref) (the ones used to correct estimates of $\boldsymbol{r}$ and $\boldsymbol{\alpha}$ while fixing that of $\sigma^{2}$), we choose

Finally, we use the approximation (ref) in Steps (ref) and (ref) of the algorithm outlined in Subsection (ref). For clarity and easy reference, we outline below the precise estimation procedure for the case of the CGMY model.

We compare the simulated performance of our estimator $(\widehat{\sigma}_{n}^{\star})^{2}$ to the estimator $\widehat{\sigma}_{n,1}^{2}$ (which could be considered the plain estimator proposed by Mies:2019) and the estimator $\widehat{\sigma}_{n,\text{JT}}^{2}$ in JacodTodorov:2014. In the latter one, we use the equation (5.3) therein with $\zeta=1.5$ and {\color{Blue} $k_{n}=252\times 6.5\times 12=19656$}, which is reasonable since the volatility is constant and there is no need to localize the estimator (so we only need one block). We take the scaling factor $u_{n}=(\ln(1/h_{n}))^{-1/30}$, as proposed in the simulation portion of Mies:2019, and $\bar{u}_{n}=(8/3)u_{n}$ for the term $S_{T}^{n}$ of equation (5.3) in JacodTodorov:2014. JacodTodorov:2014 suggests to use $u_{n}=(\ln(1/h_{n}))^{-1/30}/\sqrt{BV}$ and $\bar{u}_{n}=0.3\,u_{n}$, where $BV$ is the bipower variation, which we also tried in our simulation, but obtained worst results. In fact, we tried different parameters settings for $\zeta$, $k_{n}$, $u_{n}$, and $\bar{u}_{n}$, and select the values with the best performance. In each simulation, we divided the one-year data into $12$ months and compute the estimate of $\sigma^{2}$ for each month, and then take the average of these $12$ monthly estimators as $\widehat{\sigma}_{n,\text{JT}}^{2}$.

The results are summarized in Tables (ref)$-$(ref) for different parameter settings. The tables report the sample means, standard deviations (SDs), sample mean and SD of relative errors, and MSEs for different parameter settings based on $2000$ simulations. We also report the {TRQV} estimator, denoted by $(\widetilde{\sigma}_{n}^{\star})^{2}$, using {the threshold $\widetilde{\varepsilon}_{n}^{\,\star}$ given in (ref)} with the true values of $\sigma^{2}$, $C$, and $Y$. Finally, we also {report the TRQV} estimator, denoted by $(\sigma^{\star}_{n})^{2}$, corresponding to the true {optimal threshold} $\varepsilon^{\star}_{n}$ obtained by solving (ref) after finding $\mathbb{E}(b_{1}(\varepsilon))$ via a large scale Monte Carlo experiment.

Tables (ref)$-$(ref) show that, when $\sigma=0.2$, the MSEs of $\sigma_{0}^{2}=\widehat{\sigma}_{n,02}^{2}$, $\widehat{\sigma}_{n,1}^{2}$, $\widehat{\sigma}_{n,2}^{2}$, and $(\widehat{\sigma}_{n}^{\star})^{2}$ are getting smaller in each step. The MSE of $(\widehat{\sigma}_{n}^{\star})^{2}$ is about $81.8\%$, $71.8\%$, and $56.8\%$ lower than the MSE of $\widehat{\sigma}_{n,1}^{2}$, for $Y=1.7,\,1.5,\,1.35$, respectively, while this is $98.5\%$, $23.4\%$, and $21.7\%$ lower than the MSE of $\widehat{\sigma}_{n,\text{JT}}^{2}$, for $Y=1.7,\,1.5,\,1.35$, respectively. Similarly, as shown in Tables (ref)$-$(ref), when $\sigma=0.4$, the MSEs of $\sigma_{0}^{2}=\widehat{\sigma}_{n,02}^{2}$, $\widehat{\sigma}_{n,1}^{2}$, $\widehat{\sigma}_{n,2}^{2}$, and $(\widehat{\sigma}_{n}^{\star})^{2}$ are also getting smaller in each step. The MSE of $(\widehat{\sigma}_{n}^{\star})^{2}$ is about $35.3\%$, $47.3\%$, and $56.8\%$ lower than the MSE of $\widehat{\sigma}_{n,1}^{2}$, for $Y=1.7,\,1.5,\,1.35$, respectively. The MSE of $(\widehat{\sigma}_{n}^{\star})^{2}$ are $96.5\%$, $3.5\%$, and $58.5\%$ lower than the MSE of $\widehat{\sigma}_{n,\text{JT}}^{2}$, for $Y=1.7,\,1.5,\,1.35$, respectively. So the iterative method has a good performance and significantly improves the MSEs of the estimators $\widehat{\sigma}_{n,1}^{2}$ and $\widehat{\sigma}_{n,\text{JT}}^{2}$. Regarding the estimates of $Y$ and $C$, we notice that the second step estimates $\widehat{Y}_{n,2}$ and $\widehat{C}_{n,2}$ (obtained from fixing $\widehat{\sigma}^{2}_{n,2}$ and then applying (ref)) are significantly better than the first step estimates $\widehat{Y}_{n,1}$ and $\widehat{C}_{n,1}$ when $Y\leq 1.5$. When $Y=1.7$, there is no significant improvement.

Integrated Variance Estimation for a Stochastic Volatility Model with Jumps

In this subsection, we apply the method in the previous subsection to estimate the integrated variance under a stochastic volatility model with a CGMY jump component. We also examine the finite sample performance of the resulting estimator and compare it with the estimator of Jacod and Todorov JacodTodorov:2014. The basic idea is to split the time-period $[0,T]$ into smaller subintervals so that $\sigma$ would be approximately constant in each subinterval and, hence, $X$ is approximately L\'{e}vy within that interval. We then apply the method developed in Subsection (ref) to each subinterval to estimate the volatility level in each subinterval and finally aggregate the resulting estimates to estimate the integrated volatility.

Specifically, we consider the following Heston model:

where $(W_{t})_{t\in\mathbb{R}_{+}}$ and $(B_{t})_{t\in\mathbb{R}_{+}}$ are two independent standard Brownian motions and $(J_{t})_{t\in\mathbb{R}_{+}}$ is a CGMY L\'{e}vy process independent of $(W_{t})_{t\in\mathbb{R}_{+}}$ and $(B_{t})_{t\in\mathbb{R}_{+}}$. The parameters of the volatility specification are set as

The values of $\kappa$ and $\xi$ above are borrowed from ZhangMyklandAitSahalia:2005. In the simulation, we experiment with values of $Y=1.5$ and $Y=1.7$, and compute the estimated integrated volatility for one day under two different estimators.

We consider $5$-second observations over a one-year ($252$ days) time horizon with $6.5$ trading hours per day. We set $k_{n}=160$, which corresponds to $30$ blocks per day. As mentioned above, we treat the stochastic volatility as a constant volatility in each block, so that we can estimate the integrated volatility for each block by computing our estimator $(\widehat{\sigma}_{n}^{\star})^{2}$ with all the estimation parameters specified as in Subsection (ref). We then add the integrated volatilities for the $30$ blocks to compute our daily estimator of the integrated volatility $\int_{t}^{t+1/252}V_{s}ds$ for that day. For the estimator of {JacodTodorov:2014}, we use both equations (4.2) and (5.3) therein with $k_{n}=160$ (number of observation in each block), $\xi=1.5$, and $u_{n}=0.05(-\ln h_{n})^{-1/30}/\sqrt{BV}$, where $BV$ is the bipower variation of the previous day. To assess the accuracy of the different methods, we compute the Median Absolute Deviation (MAD) around the true value, $\int_{t}^{t+1/252}V_{s}ds$, over $200$ simulation paths.

Figure (ref) shows the estimated integrated volatility for each day computed by our MME (solid black line) and the JT estimator in JacodTodorov:2014 (dotted blue line) (see Eq. (5.3) therein), and the true daily integrated volatility (dashed red) for $1$ simulation path. When $Y=1.7$, the JT estimator tends to jitter around the true value, while our MME exhibits better performance. This behavior is further corroborated by Table (ref), which shows the MADs of both our MME and JT estimator for $6$ different days based on $200$ simulated paths. However, when $Y=1.5$, the JT estimator outperforms our MME, as shown in the right panel of Figure (ref) and Table (ref). Now, it is important to point out that the daily estimate (5.3) of JacodTodorov:2014 is based on more data than that used in our estimates. Indeed, the estimator (5.3) in JacodTodorov:2014 employs a debiasing procedure, whose debiasing term consists of two components: one that can be computed using the data in each day and another term, depending only on $Y$, that is computed using the data during the whole time horizon (in this case, one year worth of data). To explore the performance of the estimator using only contemporary data, we also analyze the performance of the estimator (4.2) in JacodTodorov:2014, which can be computed using only the data collected in each day. When $Y=1.7$ (left panel of Figure (ref)), the daily estimated integrated volatility (4.2) in JacodTodorov:2014 overestimates the true integrated volatility, and produces extremely large or small estimates at some points. We also observe this behavior when $Y=1.5$ (right panel of Figure (ref)), but the estimate is much more stable and outperforms our MME most of the time, as shown in Table (ref). To summarize, when $Y>1.5$ is large, our approach performs fairly well for the stochastic volatility model.

Acknowledgements

{\color{Red} The authors are grateful to the Editor and two anonymous referees for their multiple suggestions that help to improve the original manuscript.}