Volatility Analysis with Realized GARCH-Ito Models

TL;DR

This paper develops a unified realized GARCH-Ito model for high-frequency financial data, capturing both continuous and jump components, and proposes estimation methods validated through simulations and empirical analysis.

Contribution

It introduces the realized GARCH-Ito model embedding discrete realized GARCH in continuous volatility, with new estimation techniques and empirical validation.

Findings

Model effectively captures volatility dynamics with jumps.

Proposed estimation methods show good finite sample performance.

Empirical application demonstrates practical usefulness.

Abstract

This paper introduces a unified approach for modeling high-frequency financial data that can accommodate both the continuous-time jump-diffusion and discrete-time realized GARCH model by embedding the discrete realized GARCH structure in the continuous instantaneous volatility process. The key feature of the proposed model is that the corresponding conditional daily integrated volatility adopts an autoregressive structure where both integrated volatility and jump variation serve as innovations. We name it as the realized GARCH-Ito model. Given the autoregressive structure in the conditional daily integrated volatility, we propose a quasi-likelihood function for parameter estimation and establish its asymptotic properties. To improve the parameter estimation, we propose a joint quasi-likelihood function that is built on the marriage of daily integrated volatility estimated by high-frequency data and nonparametric volatility estimator obtained from option data. We conduct a simulation study to check the finite sample performance of the proposed methodologies and an empirical study with the S&P500 stock index and option data.

Tables

Table 1. Table 1 : The mean squared errors (MSEs) for the jump process parameters ω L subscript 𝜔 𝐿 \omega_{L} and λ 𝜆 \lambda given n = 125 , 250 , 500 , 1000 𝑛 125 250 500 1000 n=125,250,500,1000 and m = 390 , 780 , 2340 , 23400 𝑚 390 780 2340 23400 m=390,780,2340,23400 .

		MSE
		${\omega }_L$					$\lambda$
$n$ \ $m$		$390$	$780$	$2340$	$23400$		$390$	$780$	$2340$	$23400$
$125$		$5.461\times {10}^{-4}$	$1.288\times {10}^{-4}$	$1.317\times {10}^{-5}$	$4.149\times {10}^{-8}$		457.691	329.315	180.923	1.868
$250$		$5.335\times {10}^{-4}$	$1.244\times {10}^{-4}$	$1.231\times {10}^{-5}$	$3.960\times {10}^{-8}$		456.705	327.808	177.606	1.609
$500$		$5.232\times {10}^{-4}$	$1.213\times {10}^{-4}$	$1.190\times {10}^{-5}$	$3.921\times {10}^{-8}$		453.558	327.112	176.499	1.480
$1000$		$5.182\times {10}^{-4}$	$1.193\times {10}^{-4}$	$1.159\times {10}^{-5}$	$3.859\times {10}^{-8}$		450.895	325.991	175.006	1.227

Table 2. Table 2 : The mean squared errors (MSEs) for the QMLE-HL and QMLE-HLO methods on estimating realized GARCH volatility parameters for n = 125 , 250 , 500 , 1000 𝑛 125 250 500 1000 n=125,250,500,1000 and m = 390 , 780 , 2340 , 23400 𝑚 390 780 2340 23400 m=390,780,2340,23400 .

			MSE $\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$
			QMLE-HL					QMLE-HLO
$n$	$m$		${\omega }^g$	${\alpha }^g$	${\beta }^g$	$\gamma$		${\omega }^g$	${\alpha }^g$	${\beta }^g$	$\gamma$	$a$	$b$	${\sigma }_e$
125	390		22.514	83.956	401.751	80.986		6.349	70.967	220.209	73.972	13.829	7.301	2.707
	780		13.874	59.539	263.366	64.944		1.946	52.000	64.073	55.349	5.973	5.622	1.456
	2340		12.759	27.847	154.776	32.416		1.549	21.814	55.391	27.515	3.915	5.341	0.896
	23400		11.414	9.172	100.052	12.784		1.430	2.500	36.110	2.417	1.801	2.574	0.057
250	390		9.197	76.620	266.625	75.162		3.862	69.965	169.024	65.612	11.865	3.784	2.018
	780		4.645	50.045	146.061	58.639		1.106	46.422	34.844	50.483	4.224	3.189	1.426
	2340		3.604	20.791	73.631	25.116		0.850	19.384	27.946	20.633	2.154	2.947	0.723
	23400		3.089	4.571	47.478	5.838		0.762	1.356	16.774	1.209	1.135	1.557	0.029
500	390		4.552	71.620	187.886	69.883		2.633	65.360	140.817	60.363	10.524	2.300	1.895
	780		1.767	46.471	71.798	53.530		0.561	42.864	18.012	45.275	2.939	1.983	1.357
	2340		1.232	17.835	42.019	18.183		0.421	13.502	16.107	15.762	1.288	1.873	0.597
	23400		1.108	2.127	24.276	2.645		0.390	0.718	8.779	0.609	0.611	0.841	0.014
1000	390		2.544	69.202	139.960	60.467		1.808	60.128	126.474	52.889	9.530	1.694	1.646
	780		0.706	44.901	34.083	44.476		0.293	38.988	8.461	36.868	1.942	1.569	1.174
	2340		0.522	16.317	23.354	13.971		0.271	10.613	7.610	8.862	0.855	1.222	0.500
	23400		0.454	1.087	13.779	1.301		0.247	0.366	4.518	0.306	0.325	0.436	0.007

Table 3. Table 3 : The mean squared prediction errors (MSPEs) of the realized GARCH volatility predictors h i ​ ( θ ) subscript ℎ 𝑖 𝜃 h_{i}(\theta) proposed in realized GARCH-Itô model with the QMLE-HL and the QMLE-HLO methods, the GARCH volatility predictor h i ​ 0 ​ ( θ 0 ) subscript ℎ 𝑖 0 subscript 𝜃 0 h_{i0}(\theta_{0}) proposed in unified GARCH-Itô model (Kim and Wang,, 2016 ) , and the benchmark jump-adjusted MSRV method for n = 125 , 250 , 500 , 1000 𝑛 125 250 500 1000 n=125,250,500,1000 and m = 390 , 780 , 2340 , 23400 𝑚 390 780 2340 23400 m=390,780,2340,23400 .

		MSPE $\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟐}}$
		Realized GARCH-Itô		Unified GARCH-Itô	Jump-adjusted
$n$	$m$	QMLE-HL	QMLE-HLO	QMLE-HL	MSRV
125	390	4.017	3.303	7.560	7.869
	780	2.119	1.839	7.570	5.287
	2340	1.296	1.141	8.284	3.229
	23400	0.578	0.459	8.806	1.205
250	390	3.819	3.240	7.957	7.959
	780	1.990	1.715	8.088	5.346
	2340	1.206	1.035	9.182	3.284
	23400	0.500	0.438	9.593	1.231
500	390	3.657	3.101	8.127	8.004
	780	1.860	1.657	8.138	5.478
	2340	1.007	0.911	8.483	3.286
	23400	0.438	0.396	9.664	1.202
1000	390	3.501	2.998	8.052	7.963
	780	1.775	1.601	8.378	5.403
	2340	0.903	0.852	8.474	3.165
	23400	0.401	0.389	9.141	1.235

Table 4. Table 4 : The mean squared prediction errors (MSPEs) of the realized GARCH-Itô estimates with the QMLE-HL and the QMLE-HLO, the unified GARCH-Itô estimates with the QMLE-HL, and the jump-adjusted MSRV estimates.

	MSPE $\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟗}}$
Forecast Origin	Realized GARCH-Itô		Unified GARCH-Itô	Jump-adjusted
	QMLE-HL	QMLE-HLO	QMLE-HL	MSRV
$h=376$	2.527	2.323	3.141	2.655
$h=397$	3.024	2.770	3.744	3.177
$h=420$	3.851	3.510	4.766	4.040
$h=439$	5.005	4.536	6.189	5.251
$h=462$	4.052	3.913	6.813	4.134
$h=483$	6.628	5.073	12.559	6.578

Equations

$$d{X}_t={\mu }_{t}dt+{\sigma }_t(\theta )d{B}_t+L_{t}d{\Lambda }_t,$$

$${\sigma }_t^2(\theta )$$

$$L_t^2={\omega }_L+M_t,$$

$${\sigma }_n^2(\theta )=\omega +\gamma {\sigma }_{n-1}^2(\theta )+\alpha {\int }_{n-1}^n{\sigma }_s^2(\theta )ds+\beta {\int }_{n-1}^n{L}_s^{2}d{\Lambda }_s,$$

$${\int }_{n-1}^n{\sigma }_t^2(\theta )dt=h_n(\theta )+D_{n}a.s.,$$

$$h_n(\theta )={\omega }^g+\gamma h_{n-1}(\theta )+{\alpha }^g{\int }_{n-2}^{n-1}{\sigma }_s^2(\theta )ds+{\beta }^g{\int }_{n-2}^{n-1}{L}_t^{2}d{\Lambda }_t,$$

$${\omega }^g=\gamma ({\rho }_1-{\varrho }_2+2{\varrho }_3){\omega }_1-({\varrho }_1-\gamma {\varrho }_2+2\gamma {\varrho }_3){\omega }_2+(1-\gamma )\{({\varrho }_2-2{\varrho }_3)\nu +{\varrho }_2\beta \lambda {\omega }_L\},$$

$${\alpha }^g=\left({\rho }_1-{\rho }_2+2\gamma {\varrho }_3\right)\alpha ,{\beta }^g=\left({\rho }_1-{\rho }_2+2\gamma {\varrho }_3\right)\beta ,\theta =\left({\omega }^g,{\alpha }^g,{\beta }^g,\gamma \right),$$

$${\rho }_1={\alpha }^{-1}(e^{\alpha }-1),{\rho }_2={\alpha }^{-2}(e^{\alpha }-1-\alpha ),{\rho }_3={\alpha }^{-3}(e^{\alpha }-1-\alpha -\frac{{\alpha }^2}{2}),$$

$$D_n=D_n^c+D_n^J,$$

$$D_n^c=2\nu {\alpha }^{-2}{\int }_{n-1}^n\left\{\alpha (n-t-{\alpha }^{-1})e^{\alpha (n-t)}+1\right\}{Z}_{t}d{Z}_t,$$

$$D_n^J=\beta {\alpha }^{-1}\left\{{\int }_{n-1}^n\left(e^{\alpha (n-t)}-1\right)M_{t}d{\Lambda }_t+{\omega }_L{\int }_{n-1}^n\left(e^{\alpha (n-t)}-1\right)(d{\Lambda }_t-\lambda dt)\right\}$$

$$E\left[\int_{n-1}^{n}\sigma^{2}_{t}(\theta)dt\Bigg{|}\mathcal{F}_{n-1}\right]=h_{n}(\theta)\quad a.s.,$$

$$E[h_n(\theta )]=\frac{{\omega }^g+{\beta }^g\lambda {\omega }_L}{1-{\alpha }^g-\gamma },E[{\sigma }_n^2]=\frac{(\omega +\beta \lambda {\omega }_L)(1-{\alpha }^g-\gamma )+\alpha ({\omega }^g+{\beta }^g\lambda {\omega }_L)}{(1-{\alpha }^g-\gamma )(1-\gamma )},$$

$$Y_{t_{i,j}}=X_{t_{i,j}}+{ϵ}_{t_{i,j}},$$

$$L_{n,m}^{GH}(\theta )=-\sum _{i=1}^n\left[\log (h_i(\theta ))+\frac{R{V}_i}{h_i(\theta )}\right].$$

$$\begin{array}{rl}{h}_i(\theta )=& {\omega }^g+\gamma h_{i-1}(\theta )+{\alpha }^g{\int }_{i-2}^{i-1}{\sigma }_t^2(\theta )dt+{\beta }^g{\int }_{i-2}^{i-1}{L}_t^{2}d{\Lambda }_t\\ =& \sum _{l=1}^{i-1}{\gamma }^{l-1}\left\{{\omega }^g+{\alpha }^g{\int }_{i-l-1}^{i-l}{\sigma }_t^2(\theta )dt+{\beta }^g{\int }_{i-l-1}^{i-l}{L}_t^{2}d{\Lambda }_t\right\}+{\gamma }^{i-1}{h}_1(\theta ),i=2,\dots ,n.\end{array}$$

$$h_1(\theta )=\frac{{\omega }^g+{\beta }^g\lambda {\omega }_L}{1-{\alpha }^g-\gamma }.$$

$${\hat{h}}_i(\theta )=\sum _{l=1}^{i-1}{\gamma }^{l-1}\left\{{\omega }^g+{\alpha }^{g}R{V}_{i-l}+{\beta }^{g}J{V}_{i-l}\right\}+{\gamma }^{i-1}{h}_1(\theta ),i=2,\dots ,n.$$

$${\hat{L}}_{n,m}^{GH}(\theta )=-\sum _{i=1}^n\left[\log ({\hat{h}}_i(\theta ))+\frac{R{V}_i}{{\hat{h}}_i(\theta )}\right].$$

$${\hat{\theta }}^{GH}=\underset{\theta \in \Theta }{\mathrm{argmax}}\text{mbox}{\hat{L}}_{n,m}^{GH}(\theta ),$$

$$\Theta =\{({\omega }^g,{\alpha }^g,{\beta }^g,\gamma ):{\omega }_l^g<{\omega }^g<{\omega }_u^g,{\alpha }_l^g<{\alpha }^g<{\alpha }_u^g,{\beta }_l^g<{\beta }^g<{\beta }_u^g,{\gamma }_l<\gamma <{\gamma }_u,{\alpha }^g+\gamma <1\},$$

$${\Vert {\hat{\theta }}^{GH}-{\theta }_0\Vert }_{max}=O_p\left(m^{-1/4}+n^{-1/2}\right).$$

$$\sqrt{n}\left({\hat{\theta }}^{GH}-{\theta }_0\right)\stackrel{d}{\to }N\left(0,B^{-1}{A}^{GH}{B}^{-1}\right),$$

$$A^{GH}$$

$$B=\frac{1}{2}E\left[\frac{\partial h_{1}(\theta)}{\partial\theta}\frac{\partial h_{1}(\theta)}{\partial\theta^{T}}\Bigg{|}_{\theta=\theta_{0}}h_{1}^{-2}(\theta_{0})\right].$$

$${\hat{\kappa }}_T(k_{\ell })={\kappa }_T(k_{\ell })+{\epsilon }_{\ell },$$

$$N{V}_i=\frac{-2}{Tu}\mathcal{R}\left(\log \left({\hat{f}}_i(u)\wedge T\right)\right),$$

$${\hat{f}}_i(u)=1-(u^2+\sqrt{-1}u)\sum _{\ell =2}^N{e}^{(\sqrt{-1}u-1)k_{\ell -1}-\sqrt{-1}u{X}_i}{\hat{\kappa }}_T(k_{\ell -1}){\Delta }_{\ell },$$

$$N{V}_{i-1}=b+a{h}_i(\theta )+e_i,i=1,\dots ,n,$$

$${\hat{L}}_{n,m}^{GHO}(\varphi )=-\sum _{i=1}^n\left[\log ({\hat{h}}_i(\theta ))+\frac{R{V}_i}{{\hat{h}}_i(\theta )}\right]-\sum _{i=1}^n\left[\log ({\sigma }_e^2)+\frac{(N{V}_{i-1}-b-a{\hat{h}}_i(\theta ){)}^2}{{\sigma }_e^2}\right].$$

$${\hat{\varphi }}^{GHO}=\underset{\varphi \in \Phi }{\text{mbox}argmax}{\hat{L}}_{n,m}^{GHO}(\varphi ),{\hat{\theta }}^{GHO}=\text{mbox}thefirstfourcoordinatesof{\hat{\varphi }}^{GHO},$$

Full text

JEL classification: C10, C22, C58

Keywords: High-frequency financial data, option data, quasi-maximum likelihood estimation, stochastic differential equation, volatility estimation and prediction.

1 Introduction

In modern financial markets, volatility measures the degree of dispersion for assets and plays a crucial role in portfolio allocation, performance evaluation, and risk management. Low-frequency and high-frequency stock data are widely adopted to model the dynamic evolution of daily volatilities. Option data provide one more natural source for the more precise forecast of volatilities and have been investigated thoroughly since the seminal work of Black and Scholes, (1973). In traditional volatility analysis, researchers employ discrete parametric econometric models and low-frequency data. Examples include the generalized autoregressive conditional heteroskedasticity (GARCH) models (Bollerslev,, 1986; Engle,, 1982) which adopt squared daily log returns as innovations in the conditional volatilities. However, when the volatility changes rapidly to a new level, it is often difficult to catch up with the new level immediately using only the daily log returns as the innovations (Andersen et al.,, 2003). On the other hand, high-frequency financial data that refer to intra-daily observations such as tick-by-tick stock prices became available thanks to advances in information technology. Major challenges in estimating volatilities with high-frequency data are the market microstructure noises and price jumps. Without the presence of price jumps, Zhang et al., (2005) proposed two-time scale realized volatility (TSRV) which is a consistent estimator for daily variation while Zhang, (2006) further improved the TSRV to multi-scale realized volatility (MSRV) so that it can achieve the optimal convergence rate. Other forms of estimators that can achieve the optimal convergence rate only in the presence of market microstructure noises are kernel realized volatility (KRV) (Barndorff-Nielsen et al.,, 2008), quasi-maximum likelihood estimator (QMLE) (Aït-Sahalia et al.,, 2010; Xiu,, 2010), pre-averaging realized volatility (PRV) (Jacod et al.,, 2009), and robust pre-averaging realized volatility (Fan and Kim,, 2018). Empirical studies support the existence of price jumps, and decomposition of daily variation into its continuous and jump components can improve volatility forecasts (Aït-Sahalia et al.,, 2012; Andersen et al.,, 2007; Barndorff-Nielsen and Shephard,, 2006; Corsi et al.,, 2010). For example, Mancini, (2004) studied a threshold method for jump-detection and presented the order of an optimal threshold, and Davies and Tauchen, (2018) further examined a data-driven type threshold method. Also Fan and Wang, (2007) and Zhang et al., (2016) employed wavelet method to identify the jumps given noisy high-frequency data. We refer to the estimators of daily variation based on high-frequency data as the realized volatility estimators. Such estimators are more informative compared to simple squared daily log returns as the innovations, which may help to catch up with rapid changes in the volatility process better.

Efforts made for volatility estimation usually employ low- and high-frequency data independently. However, the inter-correlation between low- and high-frequency data gathered at the two different time scales cannot be ignored as low-frequency data present high-frequency data in an aggregated form. There are several attempts to bridge the gap between the two types of data. For example, multiple studies proposed new GARCH type models, which include realized volatilities as innovations in the conditional volatilities (Engle and Gallo,, 2006; Shephard and Sheppard,, 2010; Hansen et al.,, 2012). On the other hand, Wang, (2002) showed that the standard GARCH model and its diffusion limit are nonequivalent asymptotically, which discredits the direct application of statistical inferences derived for the GARCH model to its diffusion limit. Thus, Kim and Wang, (2016) introduced the unified GARCH-Itô model by embedding the standard GARCH volatility structure in the instantaneous volatilities of an Itô diffusion process. The unified GARCH-Itô model is a continuous-time process at the high-frequency timescale and when restricted to the low-frequency timescale, retains the standard GARCH structure.

In this paper, we expand the unified GARCH-Itô model (Kim and Wang,, 2016) so that features of financial data at both frequencies can be better captured as follows. First, price jumps that are well-documented in empirical studies are allowed, and we incorporate squared price jumps into the volatility dynamics by a structure similar to the ones introduced in the COGARCH model (Klüppelberg et al.,, 2004) and the jump-driven volatility model (Todorov,, 2011). Second, we embed the realized GARCH volatility structure (Hansen et al.,, 2012) in the instantaneous volatilities of a jump-diffusion process, which employs the more informative high-frequency data-based innovations. Third, the well-known intra-day U-shape volatility pattern is accounted for (Admati and Pfleiderer,, 1988; Andersen et al.,, 1997, 2019; Hong and Wang,, 2000). We name the proposed model as the realized GARCH-Itô model. The key feature of the proposed model is that its conditional volatility has integrated volatility and jump variation as innovations. Based on the structure of the conditional volatility process, we propose a quasi-likelihood function for estimating model parameters. Specifically, the quasi-likelihood function that is usually adopted in the standard GARCH type models is employed, and the realized volatility estimators are used as the proxy for conditional volatilities. We call the proposed estimator the quasi-maximum likelihood estimator based on high-frequency data and low-frequency structure (QMLE-HL). The proposed model and this estimating approach are constructed purely based on stock data. We as well harness option data to improve the model parameter estimation. In specific, Todorov, (2019) developed nonparametric volatility estimator based on a portfolio of short-dated option contracts given a general setting where jumps are present. As stated in Todorov, (2019), the estimator can be viewed as the option counterpart of high-frequency data-based volatility estimators. To incorporate the option-based nonparametric volatility estimator, we construct a joint quasi-likelihood function. We call the proposed estimator the quasi-maximum likelihood estimator based on high-frequency data, low-frequency structure and additional option data (QMLE-HLO). Both the QMLE-HL and the QMLE-HLO present good consistency and asymptotic properties. In numerical analysis, we further demonstrate that the joint estimation method QMLE-HLO performs better in estimation and prediction than the QMLE-HL.

This paper is organized as follows. Section 2 introduces the realized GARCH-Itô model. We demonstrate its connection with the realized GARCH model and discuss its advantages comparing to the unified GARCH-Itô model. Section 3 introduces quasi-likelihood estimation methods and investigates their asymptotic behaviors. Section 4 conducts a simulation study to check the finite sample performance for the proposed estimators. Section 5 carries out an empirical analysis with S&P500 stock and option data to demonstrate the advantage of the proposed model in volatility analysis. We collect all the proofs in the Appendix.

2 Realized GARCH-Itô model

The realized GARCH-Itô model is an innovated jump-diffusion process that can incorporate high-frequency based volatility model (Shephard and Sheppard,, 2010) and realized GARCH model (Hansen et al.,, 2012) structures. Let $\mathbb{R}_{+}=[0,\infty)$ and $\mathbb{N}$ be the set of all non-negative integers. Our proposed model is formulated as follows.

Definition 1**.**

Log stock price $X_{t}$, $t\in\mathbb{R}_{+}$, obeys a realized GARCH-Itô model if it satisfies

[TABLE]

[TABLE]

where $\lceil t-1\rceil$ denotes the ceiling of $t-1$, $Z_{t}=\int_{\lceil t-1\rceil}^{t}dW_{t}$, $B_{t}$ and $W_{t}$ are standard Brownian motions with respect to filtration $\mathcal{F}_{t}$ with $dW_{t}dB_{t}=\rho dt$ a.s., $\mu_{t}$ is a predictable process that is known as the drift, and $\sigma_{t}(\theta)$ is the volatility process that is adapted to $\mathcal{F}_{t}$. For the jump part, $\Lambda_{t}$ is the standard Poisson process with constant intensity $\lambda$ and $L_{t}$ denotes the i.i.d. jump sizes which are independent of the Poisson and continuous diffusion processes.

Remark 1**.**

The i.i.d. assumption on jump sizes can be rewritten as

[TABLE]

where $M_{t}$’s are i.i.d. random variables with mean zero and variance $\zeta^{2}$, $\omega_{L}+M_{t}$ is restricted to be positive. For instance, if the jump sizes $L_{t}$’s obey the Normal distribution with mean $\delta$ and variance $\eta$, then the corresponding $\omega_{L}$ takes value $\delta^{2}+\eta$ while $M_{t}$ has mean zero and variance $4\delta^{2}\eta+2\eta^{2}$.

The instantaneous volatility $\sigma^{2}_{t}(\theta)$ in (2.3) is defined at all times for $t\in\mathbb{R}_{+}$ and also retains some U-shape pattern within the intra-day. Specifically, when considering the deterministic process part of the instantaneous volatility, it is convex with respect to time $t$ and for an appropriate parameter, it has the smallest value in the middle section of the day. This U-shape instantaneous volatility pattern is often observed in empirical data and supported by financial market (Admati and Pfleiderer,, 1988; Andersen et al.,, 1997, 2019; Hong and Wang,, 2000). Moreover, random fluctuations are accounted for in the instantaneous volatility process. We note that when the process is restricted to integer times, it employs the realized GARCH model type structure (Hansen et al.,, 2012) with an additional jump innovation term as follows:

[TABLE]

where $\omega=\gamma\omega_{1}-\omega_{2}$ and $n\in\mathbb{N}$. Therefore, the instantaneous volatility process is affected by both the integrated volatilities and the jump variations of the stock price process. In comparison to the unified GARCH-Itô model (Kim and Wang,, 2016), the realized GARCH-Itô model considers price jumps, accounts for intra-day U-shape volatility pattern, and adopts a richer volatility dynamics with random fluctuations.

For statistical inferences, we study the integrated volatilities obtained from the realized GARCH-Itô model over consecutive integers, that is, $\int_{n-1}^{n}\sigma^{2}_{t}(\theta)dt$.

Proposition 1**.**

Iterative relationship exists in integrated volatilities for the realized GARCH-Itô model defined in Definition 1 and when condition (2.4) is met.

1. (a)

For $0<\alpha<1$ and $n\in\mathbb{N}$, the realized GARCH-Itô model implies that

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

are all martingale differences. 2. (b)

For $0<\alpha<1$ and $n\in\mathbb{N}$,

[TABLE]

where $h_{n}(\theta)$ is defined in (2.7). 3. (c)

For $0<\alpha^{g}+\gamma<1$ and $n\in\mathbb{N}$,

[TABLE]

where $\omega^{g}$, $\alpha^{g}$ and $\beta^{g}$ are defined in (2.8).

Proposition 1 (a) indicates that the daily integrated volatility can be decomposed into the realized GARCH volatility $h_{n}(\theta)$ and the martingale difference $D_{n}$, where the GARCH volatility $h_{n}(\theta)$ can be further explained by historical integrated volatilities and jump variations. We utilize this model feature to build up parameter estimation methods. Moreover, this paper uses the integrated volatilities as proxy to develop an estimation procedure for the GARCH parameter $\theta=\left(\omega^{g},\alpha^{g},\beta^{g},\gamma\right)$ in Section 3. This is because without the spot volatility estimation, we cannot distinguish the interceptor parameters $\omega_{1}$, $\omega_{2}$, and $\nu$.

3 Parameter estimation

In this section, we first discuss the model set-up and review nonparametric estimation methods for the integrated volatility in the presence of market microstructure noises given the jump-diffusion process. With the well-performing realized volatility and jump variation estimators, we construct quasi-maximum likelihood estimation procedures and investigate their asymptotic behaviors.

3.1 The model set-up and realized volatility estimators

Let $n$ be the total number of low-frequency observations and $m_{i}$ be the total number of high-frequency observations during the $i$th low-frequency period, for example, the $i$th day. We further denote $m=\sum_{i=1}^{n}m_{i}/n$. The underlying log price process is assumed to obey the realized GARCH-Itô model as described in Definition 1. The low-frequency data are the true log prices at integer times, $X_{i},i=0,1,\ldots,n$. The high-frequency data are observations between integer times and are contaminated by market microstructure noises. Major sources for the market microstructure noises are bid-ask bounce, discreteness of price change, and infrequent trading that only play a role in high-frequency trading (Ait-Sahalia and Yu,, 2009). We let $t_{i,j}$ be the high-frequency observed time points during the $i$th low-frequency period such that $i-1=t_{i,0}<t_{i,1}<\cdots<t_{i,m_{i}}=t_{i+1,0}=i$. In this regard, we take the well-agreed assumption in high-frequency literature such that

[TABLE]

where $\epsilon_{t_{i,j}}$’s are market microstructure noises that are some stationary random variables with $E(\epsilon_{t_{i,j}})=0$. Moreover, we note that the effect of the drift term $\mu_{t}$ on high-frequency data based volatility estimators is negligible asymptotically, so we take $\mu_{t}=0$ to highlight on modeling the volatility and jump processes.

Without the presence of price jumps, researchers have constructed nonparametric realized volatility estimators that take advantage of sub-sampling and local-averaging techniques to remove the effect of market microstructure noises so that the integrated volatility can be estimated consistently and efficiently. Such estimators include the multi-scale realized volatility estimator (Zhang,, 2006, 2011), the pre-averaging realized volatility estimator (Christensen et al.,, 2010; Jacod et al.,, 2009), and the kernel realized volatility estimator (Barndorff-Nielsen et al.,, 2008). To identify the jump locations given noisy high-frequency data, Fan and Wang, (2007) and Zhang et al., (2016) proposed wavelet methods to detect jumps and applied the MSRV method to jump-adjusted data. They demonstrated that the estimator of jump variation has the convergence rate of $m^{-1/4}$, which further helps the estimator of integrated volatility to achieve the optimal convergence rate of $m^{-1/4}$. In this paper, we let $JV_{i}$ to be the estimator of jump variation for the $i$th day and $RV_{i}$ to be the corresponding estimator of daily integrated volatility that is robust to microstructure noises and price jumps, where both estimators can achieve the convergence rate $m^{-1/4}$.

3.2 Quasi-maximum likelihood estimation based on high-frequency data and low-frequency structure

3.2.1 Estimation procedure

Recall that the integrated volatility over the $i$th period can be decomposed into the realized GARCH volatility $h_{i}(\theta)$ and martingale difference $D_{i}$ as described in Proposition 1 (a). We harness this information for making inferences on the true parameter $\theta_{0}=(\omega^{g}_{0},\alpha^{g}_{0},\beta_{0}^{g},\gamma_{0})$. Specifically, using the likelihood of the standard GARCH model and the low-frequency structure of the realized GARCH-Itô model, we define the following quasi-likelihood function

[TABLE]

Under some technical conditions, the impact of the martingale difference term $D_{i}$ is negligible in the asymptotic sense. Therefore, the realized volatility estimators $RV_{i}$’s based on data from (3.1) can be considered as the observed value for $h_{i}(\theta)$’s and are employed as the proxy. To harness the proposed quasi-likelihood function (3.2), we first need to evaluate the realized GARCH term $h_{i}(\theta)$. Recall the iterative relationship in the realized GARCH term $h_{i}(\theta)$ as described in Proposition 1 (a):

[TABLE]

The initial $h_{1}(\theta)$ is selected to be $E[h_{1}(\theta)]$ that is given in Proposition 1 (c). Specifically, we take

[TABLE]

The true integrated volatilities and jump variations are not observed so that we adopt their estimators $RV_{i}$ and $JV_{i}$, respectively. Specifically, let

[TABLE]

With the realized GARCH volatility estimator $\widehat{h}_{i}(\theta)$ in (3.3), the quasi-likelihood function (3.2) is updated to the following:

[TABLE]

We estimate the true parameter $\theta_{0}$ by maximizing the quasi-likelihood function $\widehat{L}_{n,m}^{GH}(\theta)$ in (3.4),

[TABLE]

and call the maximizer $\widehat{\theta}^{GH}$ in (3.5) the quasi-maximum likelihood estimator based on high-frequency data and low-frequency structure combined (QMLE-HL).

3.2.2 Asymptotic theory

This section establishes the consistency and asymptotic distribution for the proposed estimator $\widehat{\theta}^{GH}$. We first define some notations. For any given random variable $X$ and $p\geq 1$, define $\|X\|_{L_{p}}=\left\{E[|X|^{p}]\right\}^{1/p}$. For a matrix $A=\left(A_{i,j}\right)_{1\leq i\leq k^{\prime},1\leq j\leq k}$, let $\|A\|_{max}=\mbox{max}_{i,j}|A_{i,j}|$. Let $C$’s be positive generic constants whose values are free of $\theta$, $n$, and $m_{i}$, and may change from occurrence to occurrence. To investigate the asymptotic behaviors of proposed estimation method, we require the following technical assumptions.

Assumption 1**.**

1. (a)

Let

[TABLE]

where $\omega_{l}^{g},\omega_{u}^{g},\alpha_{l}^{g},\alpha_{u}^{g},\beta_{l}^{g},\beta_{u}^{g},\gamma_{l},\gamma_{u}$ are known positive constants. 2. (b)

We have $\underset{t\in\mathbb{R_{+}}}{\max}\mbox{ }E\left\{\sigma^{4}_{t}(\theta_{0})\right\}<\infty$ and $E(\epsilon_{t_{i,j}}^{4})<\infty$. 3. (c)

There exist some fixed constants $C_{1}$ and $C_{2}$ such that $C_{1}m\leq m_{i}\leq C_{2}m$, and $\sup_{1\leq j\leq m_{i}}|t_{i,j}-t_{i,j-1}|=O(m^{-1})$ and $n^{2}m^{-1}\rightarrow 0$ as $m,n\rightarrow\infty$ . 4. (d)

One of the following conditions is satisfied.

• (d1)

There exists a positive constant $\delta$ such that $E\left[\left(\frac{R_{i}^{2}}{h_{i}(\theta_{0})}\right)^{2+\delta}\right]\leq C$ for any $i\in\mathbb{N}$, where $R_{i}=\int_{i-1}^{i}\sigma_{t}(\theta_{0})dB_{t}$.

• (d2)

$\frac{E[R_{i}^{4}|\mathbf{\mathcal{F}}_{i-1}]}{h^{2}_{i}(\theta_{0})}\leq C$* a.s. for any $i\in\mathbb{N}$.* 5. (e)

$\sup\limits_{i\in\mathbb{N}}\left\|RV_{i}-\int_{i-1}^{i}\sigma_{s}^{2}(\theta_{0})ds\right\|_{L_{2}}\leq Cm^{-1/4}$* and $\sup\limits_{i\in\mathbb{N}}\left\|JV_{i}-\int_{i-1}^{i}L_{s}^{2}d\Lambda_{s}\right\|_{L_{2}}\leq Cm^{-1/4}$.* 6. (f)

For any $i\in\mathbb{N}$, $E\left[RV_{i}|\mathcal{F}_{i-1}\right]\leq C\,E\left[\int_{i-1}^{i}\sigma_{s}^{2}ds|\mathcal{F}_{i-1}\right]+C$ a.s. 7. (g)

$\left(D_{i},\int_{i-1}^{i}\sigma^{2}_{t}(\theta_{0})dt,R^{2}_{i}\right)$* is a stationary ergodic process.*

Remark 2**.**

The parameters of interests are related to volatilities (the 2nd moment), thus, to study their asymptotic behaviors, we require some finite 4th moment conditions such as Assumption 1 (b) and (d). Therefore, these conditions are not restrictive at all. Assumption 1 (c) is a well-known key condition in high-frequency data based volatility analysis. Under the finite 4th moment condition, Kim et al., (2016) showed that the realized volatility estimators satisfy Assumption 1 (e). Finally, the stationary ergodic condition Assumption 1 (g) is used to obtain asymptotic normality for the QMLE-HL.

The following theorems establish the convergence rate and asymptotic normality for the QMLE-HL $\widehat{\theta}^{GH}$ defined in (3.5).

Theorem 1**.**

Under Assumption 1 (a)-(f) (except for $n^{2}m^{-1}\rightarrow 0$ in Assumption 1 (c)), we have

[TABLE]

Theorem 2**.**

Under Assumption 1, we have as $m,n\rightarrow\infty$,

[TABLE]

where

[TABLE]

and

[TABLE]

Remark 3**.**

Theorem 1 shows that the convergence rate of $\widehat{\theta}^{GH}$ is $m^{-1/4}+n^{-1/2}$. The rate $n^{-1/2}$ is coming from the usual parametric convergence rate based on the low-frequency structure while the rate $m^{-1/4}$ is due to the high-frequency volatility and jump variation estimations and is known as the optimal convergence rate for estimating integrated volatilities with the presence of market microstructure noises and price jumps. Theorem 2 provides the asymptotic normal distribution for $\widehat{\theta}^{GH}$. When deriving the asymptotic normality, the condition $n^{2}m^{-1}\rightarrow 0$ in Assumption 1 (c) is imposed so that the high-frequency estimation errors of order $m^{-1/4}$ are negligible in comparison with the low-frequency estimation errors of order $n^{-1/2}$. When the condition $n^{2}m^{-1}\rightarrow 0$ is not satisfied, the asymptotic normality may depend on $m^{1/4}(RV_{i}-\int_{i-1}^{i}\sigma_{s}^{2}(\theta_{0})ds)$, which is the quantity related to high-frequency estimation. For example, if $m^{1/4}(RV_{i}-\int_{i-1}^{i}\sigma_{s}^{2}(\theta_{0})ds)$ is some martingale difference sequence, we can relax the condition $n^{2}m^{-1}\rightarrow 0$ to $nm^{-1}\rightarrow 0$. We also note that if the true stock prices are observed (i.e., without the microstructure noises), we only need the typical condition $nm^{-1}\rightarrow 0$ instead of $n^{2}m^{-1}\rightarrow 0$ to obtain the asymptotic normality (see Todorov, (2009)).

Remark 4**.**

We note that when replacing $m^{-1/4}$ in Assumption 1 (e) by $m^{-\xi}$ for some positive constant $\xi\in(0,1/4]$, the convergence rate in Theorem 1 will change to $m^{-\xi}+n^{-1/2}$. On the other hand, the condition $n^{2}m^{-1}\rightarrow 0$ in Assumption 1 (c) will be relaxed to $n^{2}m^{-4\xi}\rightarrow 0$ for deriving the asymptotic normality in Theorem 2.

3.3 Quasi-maximum likelihood estimation based on based on high-frequency data, low-frequency structure, and additional option data

3.3.1 Estimation procedure

In this section, we discuss how to incorporate additional option data information in parameter estimation. The famous Black-Scholes model indicates that option prices are determined by several factors such as time to expiration, strike price, underline asset price, and its volatility, and so one can deduce the volatility from option data. For example, the VIX presents the stock market’s general expectation of volatility. However, we usually find that the VIX is different from the historical nonparametric realized volatility. This may be because of the jumps in stock prices and the wedge between the risk-neutral and statistical probabilities. Recently, Todorov, (2019) proposed a nonparametric volatility estimator based on a portfolio of noisy short-dated option contracts with different strike prices. This estimator is robust to price jumps and does not require any assumption on the wedge between risk-neutral and statistical probabilities. Specifically, let $T$ be the time to expiration for an option contract, $k_{\ell}$ be the $\ell$th log strike price, where $k_{1}<k_{2}<\cdots<k_{N}$ and $\Delta_{\ell}=k_{\ell}-k_{\ell-1}$ for $\ell=2,\ldots,N$. Let $\kappa_{T}(k_{\ell})$ be the true option price given expiration $T$ and log-strike $k_{\ell}$. Due to observation errors in empirical derivatives pricing, the observed option price $\widehat{\kappa}_{T}(k_{\ell})$ obeys

[TABLE]

where the noises $\varepsilon_{\ell}$’s are random variables with mean zero and satisfy the technical conditions in Todorov, (2019). Given this set-up, Todorov, (2019) proposed the following nonparametric volatility estimator

[TABLE]

where

[TABLE]

$\mathcal{R}(A)$ is the real part of a complex number $A$, and $u$ is a tuning parameter.

Under some technical conditions, as $T$ goes to zero, this nonparametric volatility estimator $NV_{i}$ converges to the true spot volatility $\sigma^{2}_{i}(\theta_{0})$ (Todorov,, 2019). However, option contracts from traditional data sources such as the OptionMetrics are often quoted at the market open or close on each trading day so that the minimum choice of $T$ is $1$ business day. In this sense, $NV_{i}$ may contain integrated volatility for the remaining period from time $i$. Also Todorov, (2019) showed that the estimates $NV_{i}$’s hold a close relationship with the jump-robust realized type volatility estimates $RV_{i}$’s in his empirical study. Based on his results, we assume that the nonparametric volatility estimator $NV_{i-1}$ and the conditional daily integrated volatility $h_{i}(\theta)$ have the following linear relationship:

[TABLE]

where $b$ and $a$ are the intercept and slope coefficients, respectively. Moreover, $e_{i}$’s are martingale differences with mean zero and variance $\sigma_{e}^{2}$, and they are independent of the price process and the microstructure component.

Let $\varphi=(\omega^{g},\alpha^{g},\beta^{g},\gamma,a,b)$ and $\phi=(\omega^{g},\alpha^{g},\beta^{g},\gamma,a,b,\sigma^{2}_{e})$. Note that $\theta$ corresponds to the first four coordinates of $\varphi$ and $\phi$. We generalize (3.4) to propose the following joint quasi-likelihood function based on high-frequency and option data for estimating the true parameter $\phi_{0}=(\omega^{g}_{0},\alpha^{g}_{0},\beta^{g}_{0},\gamma_{0},a_{0},b_{0},\sigma^{2}_{e0})$

[TABLE]

We maximize $\widehat{L}^{GHO}_{n,m}(\phi)$ in (3.7) to obtain parameter estimators, that is,

[TABLE]

where $\Phi$ is the parameter space of $\phi$. We call the proposed estimator $\widehat{\phi}^{GHO}$ (or $\widehat{\theta}^{GHO}$) in (3.8) the quasi-maximum likelihood estimator based on high-frequency data, low-frequency structure, and additional option data combined (QMLE-HLO).

3.3.2 Asymptotic theory

To establish the asymptotic behaviors of the proposed estimation method, we require the following additional assumptions.

Assumption 2**.**

1. (a)

Let

[TABLE]

where $a_{l},a_{u},b_{l},b_{u},\sigma_{e_{l}}^{2},\sigma_{e_{u}}^{2}$ are known positive constants. 2. (b)

$\sup_{i\in\mathbb{N}}E\left[e_{i}^{4}\right]<\infty$. 3. (c)

$\left(D_{i},\int_{i-1}^{i}\sigma^{2}_{t}(\phi_{0})dt,R^{2}_{i},e_{i}\right)$* is a stationary ergodic process.*

The following theorems establish the convergence rate and asymptotic normality for the QMLE-HLO $\widehat{\phi}^{GHO}$ defined in (3.8).

Theorem 3**.**

Under Assumption 1 (a)–(f) (except for $n^{2}m^{-1}\rightarrow 0$ in Assumption 1 (c)) and Assumption 2 (a)–(b), we have

[TABLE]

Theorem 4**.**

Under Assumption 1 and Assumption 2, we have as $m,n\rightarrow\infty$,

[TABLE]

where

[TABLE]

[TABLE]

and $f_{i}(\varphi)=b+ah_{i}(\theta)$ for $i=1,\ldots,n$. Here $\mathbf{0}_{i\times j}$ denotes an $i$-by-$j$ matrix of zeros.

Remark 5**.**

Theorem 3 shows that the convergence rate for the QMLE-HLO is the same as the QMLE-HL. Theorem 4 provides the asymptotic normal distribution for the QMLE-HLO.

4 Simulation study

In this section, we conducted a simulation study to check the finite sample performance of the estimators $\widehat{\theta}^{GH}$ and $\widehat{\phi}^{GHO}$ given by (3.5) and (3.8) respectively, as well as to investigate the prediction performance of the realized GARCH volatilities $\widehat{h}_{i}(\widehat{\theta}^{GH})$ and $\widehat{h}_{i}(\widehat{\theta}^{GHO})$, which was also compared with the performance of the GARCH volatilities used in Kim and Wang, (2016). Here $\widehat{h}_{i}(\cdot)$ is defined in (3.3). The true log prices $X_{t_{i,j}}$, $t_{i,j}=i-1+j/m$, $i=1,\ldots,n$, $j=1,\ldots,m$, were generated based on the proposed realized GARCH-Itô model defined in (2.1) and (2.3) with the following set of parameters $\omega_{1}=5.816$, $\omega_{2}=1.228$, $\alpha=0.765$, $\beta=0.482$, $\nu=0.6$, $\gamma=0.225$, and $\rho=-0.6$. For the jump process, we took the intensity $\lambda$ to be 26 and generated $L_{t}^{2}$ such that $L_{t}^{2}=\omega_{L}+M_{t}$, where $\omega_{L}=0.005$ and $M_{t}$ follows the normal distribution with mean zero and standard deviation 0.001. Each jump $L_{t}$ was further assigned to be either positive or negative randomly. The chosen parameters resulted in the following target parameter $\theta=(\omega^{g},\alpha^{g},\beta^{g},\gamma)=(0.0122,0.717,0.452,0.225)$ for modeling the dynamics in conditional integrated volatilities. We note that the parameter $\omega^{g}$ was scaled by 10000 times compared to its empirical counterpart while the rest parameters remained the same. Scaling in this simulation study was done in order to avoid the generation of any negative value for the instantaneous volatilities due to the U-shape intra-day pattern. Initial values for the simulation were chosen to be $X_{0}=10$ and $\sigma_{0}^{2}=E(\sigma_{1}^{2})=1.4$. For the high-frequency data $Y_{t_{i,j}}$’s from (3.1), market microstructure noises were added to simulated log prices $X_{t_{i,j}}$’s between integer times, and the noises were modeled by i.i.d normal random variables with mean [math] and standard deviation $0.005$. For the option model described in (3.6), we took $a=0.812$, $b=0.072$, $\sigma_{e}=0.04$, where the intercept $b$ and standard deviation $\sigma_{e}$ were scaled by roughly 10000 times comparing to their empirical estimates. We took $n=125,250,500,1000$ and $m=390,780,2340,23400$. For each combination of $n$ and $m$, we repeated the simulation procedure for 2000 times. We followed the procedure as described in Fan and Wang, (2007) to detect the jump locations, estimate the jump variations, and compute the jump-adjusted MSRV estimators. Model parameter estimators were obtained by maximizing the proposed quasi-likelihood functions $\widehat{L}^{GH}_{n,m}(\theta)$ and $\widehat{L}^{GHO}_{n,m}(\phi)$ defined in (3.4) and (3.7), respectively.

Table 1 reports the mean squared errors (MSEs) for the jump parameters $\omega_{L}$ and $\lambda$. We find that the MSEs decrease as the number of high-frequency observations increases for each $n$, and larger $n$ often helps to locate the jumps and to estimate the parameters $\omega_{L}$ and $\lambda$ better. Table 2 presents the MSEs for the QMLE-HL and QMLE-HLO. The proposed estimating procedures present good finite sample performances and support the theoretical results derived in Section 3. For each estimation method, as the number of low-frequency or high-frequency observations increases, the MSEs decrease. When comparing the two methods, the QMLE-HLO has smaller MSE than the QMLE-HL. Thus, it is reasonable to conclude that additional option data help to enhance the estimation of model parameters.

[FIGURE:]

[FIGURE:]

The major motivation of our model proposal is to predict future volatilities by taking advantage of the imposed autoregressive type of model structure at the low-frequency. So we examined the finite sample performance of the proposed predictors $\widehat{h}_{i}(\widehat{\theta}^{GH})$ and $\widehat{h}_{i}(\widehat{\theta}^{GHO})$, where $\widehat{\theta}^{GH}$ and $\widehat{\theta}^{GHO}$ are defined in (3.5) and (3.8), respectively, and $\widehat{h}_{i}(\cdot)$ is given by (3.3). For comparison purpose, we as well investigated the prediction performance of the unified GARCH-Itô model proposed by Kim and Wang, (2016), and denote the predictor by $\widehat{h}_{i0}(\widehat{\theta}_{0}^{GH})$. Specifically, we evaluated the mean squared prediction errors (MSPEs) by

[TABLE]

where $\widehat{H}_{i}$ is one of the followings: $\widehat{h}_{i}(\widehat{\theta}^{GH})$, $\widehat{h}_{i}(\widehat{\theta}^{GHO})$, or $\widehat{h}_{i0}(\widehat{\theta}_{0}^{GH})$. As a benchmark, we as well considered the prediction of $h_{i}(\theta)$ using $RV_{i-1}$. We let the initial forecast origin to be $h=n-20$ and expanded the observation window by one low-frequency period at a time. Each time, the model parameters were estimated and the predictors were obtained.

Table 3 summarizes the MSPEs and Figure 1 presents the log MSPEs against the number of high-frequency observations. Overall, the MSPE for the realized GARCH-Itô approach decreases as the number of low-frequency or high-frequency observations increases. Moreover, the QMLE-HLO method presents the best performance regarding the MSPE. That is, the numerical results indicate that utilizing information contained in an additional data source can improve both the estimation and prediction performance of the proposed methodology. On the other hand, the unified GARCH-Itô model is not capable of explaining the rich dynamics in order to predict the conditional integrated volatilities. This may be because it takes into account neither the realized volatility nor the jump variation as an innovation. The benchmark method does not perform well because the realized GARCH-Itô model has rich dynamics that cannot be fully captured by the jump-adjusted MSRV method.

[FIGURE:]

5 Empirical analysis

In this section, we illustrate the proposed estimation methods with trading data in second for S&P500 stock index and option data quoted at the market opening on each trading day, where S&P500 stock index is the underline asset. The data sets were obtained from the TAQ and the CBOE database, respectively. We examined the period from January 3rd, 2017 to December 31th, 2018 so that the number of low-frequency periods is $n=502$. The high-frequency data are available between open and close of the market so that the number of high-frequency observations for a full trading day is $m=23400$. We followed the procedure given in Fan and Wang, (2007) to detect jumps, as well as to compute the jump variation estimates $JV_{i}$’s and jump-adjusted MSRV estimates $RV_{i}$’s. We estimated the intensity $\lambda$ by the daily averaged number of price jumps, and the parameter $\omega_{L}$ by the sample median of all squared price jumps because the sample median better described the center of the distribution formed by squared jumps. The estimated values are $\widehat{\lambda}=25.938$ and $\widehat{\omega}_{L}=3.675\times 10^{-8}$. For the option data, we followed the procedure presented in Todorov, (2019) as their empirical study covered a similar period and considered the S&P500 index as well. Specifically, we took the option contracts where the time to expiration ranges from 1 to 2 business days and skipped the contracts that were settled on a holiday. The average number of strikes per date was $62.843$ and the values of the tuning parameters were set to be the same as in Todorov, (2019). Denote the option-based nonparametric volatility estimates by $NV_{i}$’s. Figure 2 displays the auto- and cross-correlation functions (Brockwell and Davis,, 2016) for the $RV_{i}$’s, $JV_{i}$’s, and $NV_{i}$’s, which provides promising evidence for explaining the rich dynamics with these innovations. The QMLE-HL estimates are $\widehat{\omega}^{g}=1.224\times 10^{-6},\widehat{\alpha}^{g}=0.717,\widehat{\beta}^{g}=0.452$, and $\widehat{\gamma}=0.225$, and the QMLE-HLO estimates are $\widehat{\omega}^{g}=3.450\times 10^{-7},\widehat{\alpha}^{g}=0.512,\widehat{\beta}^{g}=2.375,\widehat{\gamma}=0.305,\widehat{a}=0.812,\widehat{b}=7.198\times 10^{-6},\widehat{\sigma}_{e}=4.298\times 10^{-6}$. The parameter $\omega^{g}$ denotes the intercept term in the realized GARCH volatility dynamics while the parameter $b$ denotes the intercept term in model (3.6). Their small estimated values reflect the overall level of daily volatilities that can be seen in Figure 3.

Figure 3 displays the jump-adjusted MSRV estimates, the option-based nonparametric volatility estimates, the realized GARCH volatility estimates from the QMLE-HL and the QMLE-HLO. For comparison purpose, we as well present the GARCH volatilities adopted in the unified GARCH-Itô model (Kim and Wang,, 2016). Figure 3 shows that the nonparametric jump-adjusted MSRV and the option-based nonparametric volatility estimates are both volatile, and the realized GARCH volatility estimates from the QMLE-HL and QMLE-HLO methods can account for these dynamics well. Moreover, when comparing with the unified GARCH-Itô estimates, the proposed realized GARCH-Itô estimates are closer to the jump-adjusted MSRV estimates. This may be because the realized GARCH-Itô model includes realized volatilities and jump variations as innovations while the unified GARCH-Itô model comprises squared daily log returns as innovations. That is, the proposed structure in the realized GARCH-Itô model helps to capture the market dynamics promptly.

To investigate the prediction performance of the proposed methodologies, we employed the MSPE criteria again. Denote the forecast origin by $h$. To further examine the dependency of split points, we took $h=376,397,420,439,462,483$, where each value corresponds to the last trading day of June, July, August, September, October, and November in the year of 2018. Since the exact conditional daily integrated volatilities are unknown for empirical data, we used the jump-adjusted MSRV estimates instead and evaluated the following MSPE:

[TABLE]

where $\widehat{H}_{i}$ is one of the followings: $\widehat{h}_{i}(\widehat{\theta}^{GH})$, $\widehat{h}_{i}(\widehat{\theta}^{GHO})$, $\widehat{h}_{i0}(\widehat{\theta}_{0}^{GH})$, or $RV_{i-1}$, and $\widehat{h}_{i}(\cdot)$ is defined in (3.3).

Table 4 summarizes the MSPEs from the realized GARCH-Itô, the unified GARCH-Itô, and the jump-adjusted MSRV estimates. Overall, the proposed realized GARCH-Itô estimates outperform the other methods in terms of the MSPE across various split points. When comparing the realized GARCH-Itô estimates, the QMLE-HLO presents smaller MSPE than the QMLE-HL. The empirical results indicate that the realized GARCH-Itô model holds advantages in predicting future volatilities as it utilizes the autoregressive structure in daily integrated volatilities and emphasizes high-frequency based information by using both realized volatilities and jump variations as innovations. Moreover, incorporating option-based nonparametric volatility estimates could help to predict future volatilities.

[FIGURE:]

6 Conclusion

In this paper, we introduce a novel realized GARCH-Itô model based on a jump-diffusion process which embeds the discrete realized GARCH model structure (Hansen et al.,, 2012) in its instantaneous volatility process. When the model is restricted to the low-frequency period, it employs an autoregressive type structure to explain the co-dynamics in the integrated volatilities and jump variations. Model parameters in the realized GARCH-Itô model are estimated by maximizing a quasi-likelihood function. To improve the statistical performance of the proposed estimating approach and to incorporate additional information from option data, we as well connect the nonparametric volatility estimator proposed by Todorov, (2019) with the conditional integrated volatility from the proposed model. A joint quasi-likelihood function is then adopted and we show that this method helps to improve accounting for the market dynamics in the numerical analysis.

We also leave some open issues for future study. For example, we may observe some heterogeneous variance in model (3.6). One possible approach is to generalize the homogeneous variance in (3.6) to heterogeneous variance such as replacing $\sigma_{e}^{2}$ by $\sigma_{e}^{2}h_{i}^{\zeta}(\theta)$, where parameter $\zeta>0$ is used to adjust the level of heteroscedasticity with $\zeta=0$ corresponding to the homogeneous case. We replace $\sigma_{e}^{2}$ by $\sigma_{e}^{2}\widehat{h}_{i}^{\zeta}(\theta)$ in the quasi-likelihood $\widehat{L}^{GHO}_{n,m}(\phi)$ given by (3.7) and then estimate $\zeta$ jointly with the other parameters by maximizing $\widehat{L}^{GHO}_{n,m}(\phi)$. Moreover, it is important to explore further about the optimal approach to combine and model the return and option data for volatility estimation.

Appendix A Appendix

Let $C>0$ and $0<\rho<1$ be generic constants whose values are free of $\theta$, $\phi$, $n$, and $m$ and may change from occurrence to occurrence.

A.1 Proof of Proposition 1

Proof of Proposition 1. For $k,n\in\mathbb{N}$, let

[TABLE]

By the Itô’s Lemma, we have

[TABLE]

Then simple algebraic manipulations show

[TABLE]

Since

[TABLE]

we have

[TABLE]

where $\omega^{g}$, $\alpha^{g}$ and $\beta^{g}$ are defined in (2.8). Thus, we have

[TABLE]

where $D_{n}=D_{n}^{c}+D_{n}^{J}$. Since the integrand of $D_{n}^{c}$ is predictable, $D_{n}$ is a martingale difference. Proposition 1 (b) and (c) can be showed immediately following the results of Proposition 1 (a). $\blacksquare$

A.2 Proof of Theorem 1

Maximizing $\widehat{L}^{GH}_{n,m}$ proposed in Section 3.2 is equivalent to maximizing

[TABLE]

We focus on $\widehat{L}^{GH}_{n,m}$ defined above in this proof. Define

[TABLE]

To ease notations, we denote derivatives of any given function $g$ at $x_{0}$ by

[TABLE]

Lemma 1 in Kim and Wang, (2016) shows that the dependence of $h_{i}(\theta)$ on the initial value decays exponentially. Thus, we may use the true initial value $\sigma^{2}_{0}(\theta_{0})$ during the rest of the proofs.

Lemma 1**.**

Under Assumption 1 (a)-(f), we have

1. (a)

$E\left(R_{i}^{2}\right)=E\left(\int_{i-1}^{i}\sigma_{t}^{2}(\theta_{0})dt\right)=E\left\{h_{i}(\theta_{0})\right\},$* $\sup_{i\in\mathbb{N}}E(R^{2}_{i})\leq\frac{\omega^{g}_{0}+\beta_{0}^{g}\lambda\omega_{L}}{1-\alpha^{g}_{0}-\gamma_{0}}+E(h_{1}(\theta_{0}))<\infty,$ and $\sup_{i\in\mathbb{N}}E(\sup_{\theta\in\Theta}h_{i}(\theta))<\infty;$* 2. (b)

for any $p\geq 1$,

[TABLE]

for any $j,k,l\in\{1,2,3,4\}$, where $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})=(\omega^{g},\alpha^{g},\beta^{g},\gamma)$.

Proof of Lemma 1. The statements can be showed similar to the proofs of Lemma 2 (Kim and Wang,, 2016). $\blacksquare$

Lemma 2**.**

Under Assumption 1 (a)-(d), we have

1. (a)

there exists a neighborhood $B(\theta_{0})$ of $\theta_{0}$ such that

[TABLE]

for any $j,k,l\in\{1,2,3,4\}$ where $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})=(\omega^{g},\alpha^{g},\beta^{g},\gamma)$; 2. (b)

$-\triangledown\psi_{n}^{GH}(\theta_{0})$* is a positive definite matrix for $n\geq 5$.*

Proof of Lemma 2. The proof is in the online Appendix.

Lemma 3**.**

Under Assumption 1 (a)-(f), we have

[TABLE]

Proof of Lemma 3. The proof is in the online Appendix.

Proposition 2**.**

Under Assumption 1 (a)-(d), there is a unique maximizer of $L_{n}^{GH}(\theta)$ and as $m,n\rightarrow\infty$, $\widehat{\theta}^{GH}\rightarrow\theta_{0}$ in probability.

Proof of Proposition 2. The statement can be showed similar to the proofs of Theorem 1 (Kim and Wang,, 2016) together with the result of Lemma 3. $\blacksquare$

Proof of Theorem 1. By the mean value theorem and Taylor expansion, there exists $\theta^{*}$ between $\theta_{0}$ and $\widehat{\theta}^{GH}$ such that

[TABLE]

If $-\triangledown\widehat{\psi}_{n,m}^{GH}(\theta^{*})\overset{p}{\rightarrow}-\triangledown\psi_{n}^{GH}(\theta_{0})$ which is a positive definite matrix by Lemma 2 (b), the convergence rate of $\|\widehat{\theta}^{GH}-\theta_{0}\|_{max}$ is the same as that of $\widehat{\psi}_{n,m}^{GH}(\theta_{0})$. Thus, it is enough to show

[TABLE]

and

[TABLE]

First consider $\widehat{\psi}_{n,m}^{GH}(\theta_{0})=O_{p}(m^{-1/4})+O_{p}(n^{-1/2})$. Similar to the proofs of Theorem 2 (Kim and Wang,, 2016), we can show that

[TABLE]

By the application of the Itô’s lemma and Itô’s isometry, we can show for any $j\in\{1,2,3,4\}$,

[TABLE]

where the last inequality is due to Lemma 1 (b). Similar to the proofs of Theorem 2 (Kim and Wang,, 2016) together with the results of Lemma 2 and Proposition 2, we can show

[TABLE]

$\blacksquare$

A.3 Proof of Theorem 2

Proof of Theorem 2. By the mean value theorem and Taylor expansion, we have for some $\theta^{*}$ between $\theta_{0}$ and $\widehat{\theta}^{GH}$,

[TABLE]

where the second equality is due to (A.4). By the ergodic theorem and the result in the proof of Theorem 1, we have

[TABLE]

and $B$ is a positive definite matrix. For any $f\in\mathbb{R}^{4}$, let

[TABLE]

Then $d_{i}$ is a martingale difference with $E(d_{i}^{2})<\infty$.

Since $\left(D_{i},\int_{i-1}^{i}\sigma_{t}^{2}(\theta_{0})dt,R_{i}^{2}\right)$’s are stationary and ergodic processes, $d_{i}$ is also stationary and ergodic. By the martingale central limit theorem and Cram$\acute{\text{e}}$r-Wold device, we have

[TABLE]

Therefore, by Slutsky’s theorem, we conclude that

[TABLE]

$\blacksquare$

A.4 Proof of Theorem 3

Maximizing $\widehat{L}^{GHO}_{n,m}$ is equivalent to maximizing

[TABLE]

where $\widehat{f}_{i}(\varphi)=b+a\widehat{h}_{i}(\theta)$. We focus on $\widehat{L}^{GHO}_{n,m}$ defined above in this proof. Define

[TABLE]

and

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

and

[TABLE]

Lemma 4**.**

Under Assumption 1 (a)–(f) and Assumption 2 (a)–(b),

1. (a)

there exists a neighborhood $B(\phi_{0})$ around $\phi_{0}$ such that

[TABLE]

for any $j,k,l\in\left\{1,2,\ldots,7\right\}$, where $\phi=(\phi_{1},\phi_{2},\phi_{3},\phi_{4},\phi_{5},\phi_{6},\phi_{7})=(\omega^{g},\alpha^{g},\beta^{g},\gamma,a,b,\sigma^{2}_{e})$; 2. (b)

$-\triangledown\psi_{n}^{GHO}(\theta_{0})$* is a positive definite matrix for $n\geq 7$.*

Proof of Lemma 4. The proof is in the online Appendix. $\blacksquare$

Lemma 5**.**

Under Assumption 1 (a)-(f) and Assumption 2 (a)–(b), we have

[TABLE]

Proof of Lemma 5. The proof is in the online Appendix. $\blacksquare$

Proposition 3**.**

Under Assumption 1 (a)-(f) and Assumption 2 (a)–(b), there exists a unique maximizer for $L_{n}^{GHO}(\phi)$. As $m,n\rightarrow\infty$, $\widehat{\phi}^{GHO}\rightarrow\phi_{0}$ in probability, where $\phi_{0}$ is a vector of true parameters.

Proof of Proposition 3. According to the definition of $L_{n}^{GHO}(\phi)$, we have

[TABLE]

Then, similar to the proofs in Theorem 1 of Kim and Wang, (2016), we can show the uniqueness of the solution of $L_{n}^{GHO}(\phi)$, which together with Lemma 5 implies Proposition 3. $\blacksquare$

Proof of Theorem 3. By the mean value theorem and Taylor expansion, we have

[TABLE]

where $\phi^{*}$ is between $\phi_{0}$ and $\widetilde{\phi}^{GHO}$. According to Lemma 4 (b), $-\nabla\psi_{n}^{GHO}(\phi_{0})$ is a positive definite matrix. If $-\nabla\widehat{\psi}_{n,m}^{GHO}(\phi^{*})\xrightarrow{p}-\nabla\psi_{n}^{GHO}(\phi_{0})$, then the convergence rate of $\widehat{\phi}^{GHO}-\phi_{0}$ is the same as the convergence rate of $\widehat{\psi}_{n,m}^{GHO}(\phi_{0})$.

By the similar arguments in the proof of Theorem 1, we can show

[TABLE]

We have

[TABLE]

The arguments in the proof of Theorem 1 shows that the first term of the right side of (A.6) is $O_{p}(n^{-1/2})$. Since $e_{i}$ is independent of $\frac{\partial f_{i}(\varphi_{0})}{\partial\varphi}$, the second term of the right side of (A.6) is also $O_{p}(n^{-1/2})$. Thus, the convergence rate of $\widehat{\psi}_{n,m}^{GHO}(\phi_{0})$ is $n^{-1/2}+m^{-1/4}$.

Similar to the proof of Theorem 1, we can show

[TABLE]

Therefore, the statement is proved. $\blacksquare$

A.5 Proof of Theorem 4

Proof of Theorem 4. Since the mean value theorem and Taylor expansion provides

[TABLE]

where $\phi^{*}$ is between $\phi_{0}$ and $\widehat{\phi}^{GHO}$, we have

[TABLE]

where the equality can be showed similar to the proof of Theorem 1. Since $e_{i}$ is independent of $D_{i}$ and $\left(D_{i},e_{i},Z_{i}^{2}\right)$ is stationary and ergodic, by the Cramér-Wold device and the martingale central limit theorem, we have

[TABLE]

On the other hand, we have

[TABLE]

Therefore, by the Slutsky’s theorem, we have

[TABLE]

$\blacksquare$

Acknowledgements

The research of Xinyu Song was supported by the Fundamental Research Funds for the Central Universities (2018110128), China Scholarship Council (201806485017), and National Natural Science Foundation of China (Grant No. 11871323). The research of Donggyu Kim was supported in part by KAIST Settlement/Research Subsidies for Newly-hired Faculty grant G04170049 and KAIST Basic Research Funds by Faculty (A0601003029). The research of Huiling Yuan was supported by the State Scholarship Fund. The research of Xiangyu Cui was supported by National Natural Science Foundation of China (71671106). The research of Zhiping Lu was supported by Natural Science Foundation of Shanghai (17ZR1409000) and the 111 Project (B14019). The research of Yong Zhou was supported by the National Natural Science Foundation of China (71931004 and 91546202). The research of Yazhen Wang was supported in part by NSF Grants DMS-15-28735, DMS-17-07605, and DMS-19-13149.

We thank the Associate Editor, Viktor Todorov, and two anonymous referees for many constructive suggestions that have significantly improved the paper.

This research was performed using the compute resources and assistance of the UW-Madison Center For High Throughput Computing (CHTC) in the Department of Computer Sciences. The CHTC is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the National Science Foundation, and is an active member of the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy’s Office of Science.