State Space Model of Realized Volatility under the Existence of Dependent Market Microstructure Noise

State Space Model of Realized Volatility under the Existence of Dependent Market Microstructure Noise

\titlerunning{The model of RV under the Existence of Dependent MN}

Introduction

Volatility means the degree of variation of a stock price. Integrated volatility (IV) is an integral of continuously changing instantaneous volatility over a specified period, for example, a day. IV is important in finance, for example, option pricing, risk management, and optimal portfolio construction. Volatility is an unobservable variable and needs to be estimated from observable data such as stock price. There are classical volatility models, Generalized Autoregressive Conditional Heteroskedastic (GARCH) model and Stochastic Volatility model. These are models of autocovariance structure of volatility and can estimate the volatility based on the model from daily return. For these, see Bollerslev et al[2010] for example.

In recent years, Realized volatility (RV) has attracted attention as a method of estimating IV. RV is calculated using high frequency data such as seconds and minutes return, and is a model-free estimator that does not assume volatility model like GARCH or SV model. IV can be observed by RV, and modeling that can describe the statistical properties of RV and developments of improvement statistics of RV closer to IV are studied. Research on RV has dramatically increased in recent years due to expansion of availability of high frequency data. In regard to the survey of RV modeling, for example, see MacAleer and Medeiros[2008] and its references.

Since RV is based on the quadratic variation, in order to make estimation that is closer to IV, we have to increase the frequency data used for RV calculation. However, it is actually confirmed that when the RV is calculated using the ultra high frequency data, the RV value skyrocket. The reason why this phenomenon occurs is explained as follows. The theoretical correctness of RV is assumed that there is no measurement error in the data, but the actual data includes observation error called market microstructure noise (MN), whereby RV includes bias. (Zhou[1996]) An improved estimator of RV that converges to IV under the existence of MN assumed to be independent noise. However, empirical studies of Hansen and Lunde[2006] have reported that MN is not independent: auto-correlate and correlate with return. Under this “dependent” MN, bias remains in the estimator of Zhou[1996]. This was solved by the Realized Kernel (RK) of Barndorff-Nielsen et al[2011], and thereafter, study to improve this convergence rate and empirical research using RK are continuing. For the survey of the improved estimators of RV, see Mukherjee et al[2019] for example.

Barndorff-Nielsen and Shephard[2002] proposed a model approach of IV estimation with RV as observation. In other words, by formulating the RV as a state space time series model representing the sum of IV and observation error, it is possible to estimate and predict IV based on the model from RV as data by Kalman filter. Meddahi[2003] generalizes this model and assumes a general class continuous time SV model (Square root stochastic autoregressive volatility model) for instantaneous volatility and derives model of IV and clarifies the relationship between SV model and IV parameters.

Nagakura and Watanabe[2009, 2015] expanded this model to the case where the observation price includes MN which is an independent noise and derives the covariance structure of the bias of RV due to MN. As a result, the bias follows the MA$(1)$ process. They further derive the relationships between this bias and the MN parameters, uniquely estimate the parameters of the SV model and MN from the RV data and state space model and propose a method to estimate and predict IV and bias due to MN.

In this paper, we will relax the assumption of Nagakura and Watanabe[2009, 2015], and try to extend their model assuming Hansen and Lunde [2006] type dependent MN.

Integrated Volatility and Realized Volatility

Let $p (t)$ be the logarithmic value of the asset price at time $t$ and let this follow the next Ito process.

where $W(t)$ is standard Brownian motion and $\mu(t), \sigma(t)$ are c\`adl\`ag and $\mathcal{F}(t)$-adapted process.

Here, the Integrated Volatility (IV) at time $t$ is defined as

where the unit of $t$ is determined depending on research objective. In this paper, $t$ is interpreted as a day. $IV_t$ is the theoretical true volatility, it represents the price variability at time $t$, it is an important risk indicator in the finance. However, this is unobservable. Our aim is to estimate/predict this IV from observable data.

As data, $m$ intraday returns $r_t$ discretely (equal intervals) observed at time $t$. where

By intraday returns, the Realized Volatility (RV) at time $t$ is calculated as

Under the above, from the classical stochastic calculus, we can see that the following holds. (see, e.g. Barndorff-Nielsen and Shephard[2002].)

Market Microstructure Noise

On the above, RV is the consistent estimator of IV, but it is confirmed that if RV is actually calculated in a situation where $ m $ is quite large, the value of RV will skyrocket.

The following figure shows the time interval set as 30 seconds$(m=2880)$ , 1 minute$(m=1440)$, 2 minute$(m=720)$, 5 minute$(m=288)$, 10 minute$(m=144)$, 20 minute$(m=72)$, 25 minute$(m=57.6)$ and 30 minute$(m=48)$ respectively, using the dollar/yen rate data from May 1, 2009 to April 29, 2016, calculate the RV and plot the average.

Since the average of RV increases in proportion to $m$, the upward bias against RV can be confirmed, suggesting a gap between the theory of the previous section and the reality. Zhou[1996] proposed market microstructure noise (MN) as a model to explain this phenomenon. In this model, the actual observation price $p^*_t$ deviates from true log price $p(t)$ of the equation ((ref)) by MN as follows.

where $\varepsilon_t$ is MN and is explained as observation error due to market microstructure. The observation intraday return $r^{*(m)}_t$ calculated from $p^*_t $ is as follows:

where $\varepsilon_t-\varepsilon_{t-\frac{1}{m}}=e^{(m)}_t$. $RV^{*(m)}_t$ calculated by observation intraday return $r^{*(m)}_t$ is expressed as follows:

where $\displaystyle u_t^{(m)}=2\sum^m_{i=1}r^{(m)}_{t-1+i/m}e^{(m)}_{t-1+i/m}+\sum^m_{i=1}e^{(m)2}_{t-1+i/m}.$

This $RV^{*(m)}_t$ is called noise-contaminated RV (NCRV). As mention above, when $m\to\infty$ is $RV^{(m)}_t \to IV_t \ {\rm in \ P}$, but $u_t^{(m)} \nearrow \infty$ can also be proved. So if $m$ is large, the influence of bias $u_t^{(m)}$ by MN will be large. On the other hand, if $m$ is small, then the variance of the estimation error becomes large due to the decrease in the number of samples.

Assumption of asset price

The asset price model of Assumption (ref) is a general model and includes representative multi factor stochastic volatility models e.g. Heston model and GARCH diffusion model. (see, Meddahi[2003])

Meddahi[2003] showed that under the Assumption (ref), $IV_t$ follows the ARMA$(p, p)$ process:

where ${\rm WN}(\sigma_\eta^2)$ is a white noise of variance $\sigma_\eta^2$\footnote{$\eta_t$ is white noise is that: ${\rm E}[\eta_t]=0, {\rm Var}[\eta_t]=const<\infty , {\rm Cov}[\eta_t,\ \eta_{t-s}]=0,\ (s\neq 0).$}.

RV has an estimated error $d_t^{(m)}$ for IV. Barndorff-Nielsen and Shephard[2002] suggested that $d_t^{(m)}$ is white noise and proposed a model to estimate and predict IV from RV as follows. From the equation ((ref)),

where $d_t^{(m)}, \eta_t$ are the white noise of variance $\sigma_d^{(m)2}, \sigma_\eta^2$. Equation((ref)), ((ref)) are called BSM($p$) model. The parameters $(c_{IV}, \{\phi_i\}_1^p, \{\theta_i\}_1^p, \sigma^{(m)2}_d, \sigma^2_\eta)$ can be estimated by the Kalman filter and the quasi-maximum likelihood method, and we can estimate and predict $IV_ {1: T}$ from $RV^{(m)}_{1: T}$ by Kalman smoother and Kalman filter.

Assumption of MN

Previous study (Nagakura and Watanabe[2015])

Nagakura and Watanabe[2015] assumed the following for MN $\varepsilon_t$.

Under this assumption, Nagakura and Watanabe[2009, 2015] proved the following.

From the above, the observed $RV_t^{*(m)}$ is decomposed into a true process $IV_t$, its estimation error $d_t^{(m)}$ and the bias due to MN $u_t^{(m)}$. From the equations((ref)), ((ref)), it is represented by the following linear state space model.

This is called the NW$(p)$ model. The parameters\\ $(c_{IV}, \{\phi_i\}_1^p, \{\theta_i\}_1^p, c_u^{(m)}, \theta^{(m)}_u, \sigma^{2(m)}_d, \sigma^{2(m)}_\eta, \sigma^{2(m)}_\xi)$ can be estimated by the Kalman filter and quasi-maximum likelihood method. We can estimate and predict of $IV_{1:T}, u_{1:T}^{(m)}$ from $RV^{*(m)}_{1:T}$ by Kalman smoother and Kalman filter.

Nagakura and Watanabe[2009] did an empirical study using high frequency data of dollar/yen rate and found that, in the case of $p=1,2$, they compared the IV estimation accuracy and the 1-ahead prediction accuracy of the NW$(p)$ model and BSM$(p)$ model, and NW(2) was the best in their results.

Dependent MN

According to Hansen and Lunde[2006], MN has the following properties:

So, MN is NOT independent noise. We assume MN satisfying HL.1 and HL.2 and derive the properties of $u_t^{(m)}$.

Model identification

Estimation of ${\rm E}\left[u_t^{(m)}\right], \sigma^2$

From Theorem (ref), $r_t^*\sim{\rm MA}(q+1)$. Let its auto-covariance be $G_n$.

From these,

Also, the following holds.

From these, we can obtain $\sigma^2, {\rm E}[u_t]$ from $r_t^*$ and $RV_t^*$.

Estimation of ${\rm Cov}[\varepsilon^{(2)}_t, \varepsilon^{(2)}_{t\pm n}]$

From equation ((ref)), $\Omega_i, (i=1,\ldots, q)$ can be calculated in order, and the following can be known inductively.

Model identification and parameter estimation of $u_t^{(m)}$ in case of $f(m)=0$

When $f(m)=0, \sigma_\delta=0$, that is, when MN follows normal MA$(q)$ process, from Theorem (ref) ${\rm Cov}\left[u_s^{(m)}, d_t^{(m)}\right]={\rm Cov}\left[u_s^{(m)}, IV_t\right]=0, \forall s,t$ and ${\rm Cov}\left[u_t^{(m)}, u_{t\pm n}^{(m)}\right]=mC_{mn}^{(3)}+\sum^{m-1}_{k=1}(m-k)(C_{mn+k}^{(3)}+C_{mn-k}^{(3)})$ . Next, we'll calculate $n$ that ${\rm Cov}\left[u_t^{(m)}, u_{t\pm n}^{(m)}\right]=0$. In the above equation, when $k=m-1$, $(m-m+1)C_{mn-m+1}^{(3)}=C_{m(n-1)+1}^{(3)}$ is the last term which becomes $0$ with respect to $n$ and since $C^{(3)}_n=0, n\geq q+2$ holds, $m(n-1)+1\geq q+2,$ thus, $n\geq \cfrac{q+1}{m}+1.$ So, ${\rm Cov}\left[u_t^{(m)}, u_{t\pm n}^{(m)}\right]$ becomes $0$ for the first time $n=\left\lceil \cfrac{q+1}{m}\right\rceil+1$. ($\lceil\cdot\rceil$ is round up.) that is,

Therefore, $u_t^{(m)}\sim$MA$\left(\left\lceil \cfrac{q+1}{m}\right\rceil\right)$.

Nagakura and Watanabe[2015] proved $u_t^{(m)}\sim$MA$(1)$, assuming MN is IID noise. This is the case when $q=0$ in Assumption (ref),which is $\left\lceil \cfrac{q+1}{m}\right\rceil=1.$ This is true when $q\leq m-1$. On the contrary, $u_t^{(m)}\sim$MA$(Q), Q\geq 2$ is when $q> m-1$, which means that the autocorrelation of MN continues for more than one day. In the empirical study of Ubukata and Oya[2009], autocorrelation of MN continues for about 60 seconds long, so in reality it can be considered to be $q\leq m-1$.

Now if $q\leq m-1$ then $u_t^{(m)}\sim{\rm MA}(1)$ so, there are parameters $c_u^{(m)}, \theta^{(m)}_u, \sigma^{2(m)}_\xi$ that can be expressed as

Under this expression, the following holds.

On the other hand, from Theorem (ref),

where $C_n^{(3)}$ is in the Theorem (ref) with $\sigma_\delta=0, f(m)=0$.

Since $\sigma^2$ and $\Omega_i, i=0,\ldots, q$ can be calculated from $r_t^*$, the left side of Equation ((ref)) is obtained from $r_t^*$. From the second and third equations of Equation ((ref)), letting ${\rm Var}\left[u_t^{(m)}\right]=:v, {\rm Cov}\left[u_t^{(m)}, u_{t\pm 1}^{(m)}\right]=:c$, $ v=(1+\theta^{(m)2}_u)\sigma^{2(m)}_\xi, c=\theta^{(m)}_u\sigma^{2(m)}_\xi. $ Since $\sigma^{2(m)}_\xi>0$ and $v^2-4c^2=8\sigma^2(\Omega_0-\Omega_1)+8m(\Omega_0-\Omega_1)^2\geq 0$, $ \theta^{(m)}_u=\cfrac{v\pm\sqrt{v^2-4c^2}}{2c}. $ In order to satisfy the stationary condition of $u^{(m)}_t$:$|\theta^{(m)}_u|<1$, let $ \theta^{(m)}_u=

$ be set. And $ \sigma^{2(m)}_\xi=\cfrac{c}{\theta^{(m)}_u}. $ To summarize the above, it is as follows:

The parameters of linear state space model ((ref)), ((ref)), ((ref)) and ((ref)) can be estimated as follows.

Model identification and parameter estimation of $u_t^{(m)}$ in case of $f(m)=o\left(m^\frac{1}{2}\right)$

Here, we assume $f(m)=o\left(m^{\frac{1}{2}}\right)$, $\sigma_\delta=0$ and $q<m$. Furthermore, to make explicit $\displaystyle {\rm Cov}\left[\int^t_{t-\frac{1}{m}}\sigma(s)ds, \int^{t\pm\frac{n}{m}}_{t\pm\frac{n\pm1}{m}}\sigma(s)ds\right]$ and $\displaystyle {\rm Cov}\left[\int^{t-1+\frac{i}{m}}_{t-1+\frac{i-1}{m}}\sigma^2(u)du, \int^{t-1\pm n+\frac{j}{m}}_{t-1\pm n+\frac{i\pm1}{m}}\sigma(u)du\right]$ that appear in Theorems (ref) and (ref). These can not be obtained under the assumption (ref). Because SR-SARV model defines $\sigma^2(t)$, but it requires $\sigma(t)$ for these calculations.

So, we use an approximation technique called the delta method. Function $g(x,y)$ that depends on two variables, $x,y$, the Taylor series to second order about the point (a, b) is $$g(x,y)\approx g(a,b)+\left[(x-a)\frac{\partial}{\partial x}+(y-b)\frac{\partial}{\partial y}\right]g(a,b)+\frac{1}{2}\left[(x-a)\frac{\partial}{\partial x}+(y-b)\frac{\partial}{\partial y}\right]^2g(a,b).$$ Let $x, y$ are random variables$X, Y$, $a, b$ are ${\rm E}[X], {\rm E}[Y]$ and take expected values on both sides, $${\rm E}[g(X,Y)]\approx g(a, b)+\frac{1}{2}g_{xx}(a,b){\rm Var}(X)+\frac{1}{2}g_{yy}(a,b){\rm Var}(Y)+g_{xy}(a,b){\rm Cov}[X,Y]$$

Let $g(x,y)=\sqrt{xy}$, we have

Here, using the approximation $\int^{t\pm \frac{n}{m}}_{t\pm \frac{n\pm 1}{m}} \sigma(s) ds\approx \frac{1}{m}\sigma(t\pm \frac{n}{m})$, we have

Then, from Theorem (ref),

Also, the aforementioned approximation $${\rm Cov}\left[\int^t_{t-\frac{1}{m}}\sigma(s)ds, \int^{t\pm\frac{n}{m}}_{t\pm\frac{n\pm1}{m}}\sigma(s)ds\right]\approx \frac{1}{m^2}{\rm Cov}[\sigma(t\pm \tfrac{n}{m}), \sigma(t)]=O(\{f(m)\}^2m^{-2}\exp(\lambda_{\rm min} \tfrac{n}{m}))$$ where $\lambda_{\rm min}=\min_{i=1,\ldots, p}\lambda_i$.

Using this, we approximate part $C^{(1)}$ of ${\rm Cov}\left[u_t^{(m)}, u_{t\pm n}^{(m)}\right]$ in Theorem (ref):$4mC_{mn}^{(1)}+4\sum^{m-1}_{k=1}(m-k)(C_{mn+k}^{(1)}+C_{mn-k}^{(1)})$ as follows.

Similarly, $C^{(4)}_n$ and $C^{(5)}_n$ in Theorem (ref) are as follows:

Finally, applying the above to Theorems (ref),(ref), we have

where

where $\gamma_n=2(4\Omega_{n}^2+\Omega_{n-1}^2+\Omega_{n+1}^2-4\Omega_{n}\Omega_{n+1}-4\Omega_{n-1}\Omega_{n}+2\Omega_{n-1}\Omega_{n+1}),$ and also,

Ignoring $o\left(1\right)$, when $f(m)$ and ${\rm E}\left[\sigma(t)\right]$ are given, from previous subsection, using the data of $r_t^*$, $\Omega_i, i=0,\ldots, q$ can be calculated. And from the above equation ${\rm Var} \left[u_t^{(m)}\right], {\rm Cov}\left[u_t^{(m)}, u_{t\pm 1}^{(m)}\right]$ can be calculated.

Here, ${\rm E}[\sigma(t)]$ appears. So, we use next delta method for ${\rm E}[\sigma(t)]$.

Random variable $X$, for second order differential function $g$, consider the following Taylor expansion around the average ${\rm E}[X]$ of $X$:

Taking the expected value on both sides, we have

We use this to compute ${\rm E}[\sigma(t)]$.

By this approximation, ${\rm E}[\sigma(t)]$ can be represented by $\sigma^2, \{\omega^2_i\}_1^p$. When parameter estimation is performed below, the above expression is used instead of ${\rm E}[\sigma(t)]$.

From Equation ((ref)), we have

For $\Omega_0>0$, we restrict the range of $f(m), {\rm E}[\sigma(t)]$ as follows:

where $A=-2mG_0-4m\Omega_1+2\sigma^2.$ Then, the parameters of linear state space model ((ref)), ((ref)), ((ref)) and ((ref)) can be estimated as follows:

Model comparison

We compare the IV estimation and 1-ahead prediction accuracy of the NW model and the proposal model, according to Nagakura and Watanabe[2015], by the following method.

Comparison by simulation

We use this $p^*_t$ as data and calculate NCRV from data. The following figure (ref) is a plot of calculated NCRV (red dots) and true IV (gray line).

NCRV is greatly influenced by bias due to MN, and is estimated to be overall larger than IV.

Estimate the parameters of the NW model and proposal model with the factor number $p=1$, using $p^*_t$ every minute for $N$ days and $RV_t^{*(1440)}$ calculated from data. In the case of NW model, all parameters are estimated by QMLE according to the method of Nagakura and Watanabe[2015] direct approach (p55) with $RV_t^{*(1440)}$ as data. In the case of proposal model, zero-f and weak-f estimate the parameters by methods of previous section with $RV_t^{*(1440)}$ and $r_t^*$ as data. The following table (ref) is the estimation result of each model and comparison with a true value.

Estimate and predict IV according to each estimated model and compare the IV estimation error with the 1-ahead prediction error. Figures (ref), (ref) and (ref) are plots of the estimation result of $IV_t$ and $u_t^{(m)}$ in the NW, zero-f and weak-f respectively. The gray line is the true value, and the red line is the estimated value in the model.

In NW model, the averages of $IV_t$ and $u_t^{(m)}$ ($\sigma^2$ and $c_u^{(m)}$) are estimated by QMLE. Since this accuracy is bad, $IV_t, u_t^{(m)}$ are shifted from the true value.\footnote{When maximizing the quasi-likelihood function, we tried changing the initial value, but $\sigma^2$ and $c_u^{(m)}$ could not make it converge to “good” value.} In the proposal model (zero-f model and weak-f model), these are calculated by Equation ((ref)) and ((ref)) from data $r_t^*$ which is also estimated to be close to the true value. However, the variance of $u_t^{(m)}$ seems to be underestimated, the true variance of $u_t^{(m)}$ is 3.6641$\times 10^{-4}$, but the value estimated by zero-f model is 1.258$\times 10^{-5}$, weak-f model is 1.819$\times 10^{-5}$. NW model is 1.225$\times 10^{-4}$, which is also underestimated, but it is closest to the true value. This part is a term $4mC_0^{(1)}+8\sum^{m-1}_{k=1}(m-k)C_k^{(1)}+4mC_0^{(2)}+4\sum^{m-1}_{k=1}(m-k)C_k^{(2)}$ in the Theorem (ref), where HL.1 contributes to the variance of $u_t^{(m)}$, and zero-f model ignores this term with $f(m)$, so this part is underestimated. In the weak-f model, somehow $f(m)$ is estimated to be smaller than the true value, so it is close to the result of the zero-f model.

The tables (ref) and (ref) are comparison results of estimation accuracy and prediction accuracy.

In Table (ref), RK is added as the IV estimation method. Barndorff-Nielsen et al[2011] constructed an estimator, Realized Kernel $RK_t^{(m)}$ which has asymptotic normality under the existence of MN satisfying HL.1, HL.2 and HL.3:

and Kernel function $K(x) : \mathbb{R}^1\to\mathbb{R}^1$ is $C^2$ and satisfies $K(0)=1, K'(0)=0$ and appropriate definiteness and integrability. The parameter $H$ is called bandwidth, when $H\propto m^{\frac{3}{5}}$, $RK_t^{(m)}$ has asymptotic normality with convergence rate $m^{\frac{1}{5}}$. Barndorff-Nielsen et al[2011] recommended $K(x)$ be Parzen kernel:

and they proved bandwidth is optimized as follows: $H=c\xi^{\frac{4}{5}}m^{\frac{3}{5}}, c=3.5134,\xi^2=\cfrac{{\rm Var}[\varepsilon_t]}{\sqrt{T\int^T_0\sigma^4(t)dt}}$. However, since an unobservable variable is used for $\xi^2$, the optimal value is not obtained. Here, quote Ubukata and Watanabe[2011], estimated every day with $\widehat{{\rm Var}[\varepsilon_t]}=\cfrac{RV_t^{(m)}}{2m}$ and $\int^T_0\sigma^4(t)dt$ is estimated every day by Realized Quarticity(Barndorff-Nielsen and Shephard[2002]):$\displaystyle RQ_t^{(m)}=\cfrac{m}{3}\sum_{i=1}^mr^{(m)4}_{t-1+\frac{i}{m}}$ with $m=96$. We call $H$ estimated by this method “rough". In Table (ref), “RK” using “rough” bandwidth $H$ and “RK (true $H$)” using the optimal bandwidth $H$ calculated by using true $\sigma^2(t), \varepsilon_t$, which is actually impossible to obtain.

In Table (ref), for MSE and $R^2$, zero-f model and weak-f model are equal and the best. For QLIKE, RK(true $H$) is the minimum and zero-f model is the second. From these results, it can be seen that the proposal models zero-f model and weak-f model dominates the existing method NW model and RK in the IV estimation accuracy.

In Table (ref), for MSE, zero-f model and weak-f model are equal and minimum. For QLIKE, weak-f model is minimum. For $R^2$, zero-f model is the best. Thus, proposal models zero-f model and weak-f model dominates the NW model in the IV 1-ahead prediction accuracy.

Model comparision by actual data

Data

Data are taken from the middle price (1749 days) of the dollar-yen rate every 30 seconds (m=2880) from May 1, 2009 to April 29, 2016. \footnote{The data are extracted every 30 seconds from the tick data downloaded from {\url{https://pepperstone.com/en/}}.} The sample size is 4905910. Following Anderson et al[2001], Nagakura, Watanabe[2009, 2015], we remove from the data the trading days that satisfy the following conditions as inactive trading days

The following 184 days were removed from the data. All dates in black are Sundays, and dates in red are other days of the week.\\[0.3cm] {\scriptsize 2009-05-03, 2009-05-10, 2009-05-17, \textcolor{red}{2009-05-21}, 2009-05-24, 2009-05-31, 2009-06-07, 2009-06-14, 2009-06-21, 2009-06-28, 2009-07-05, 2009-07-12, 2009-07-19, 2009-07-26, 2009-08-02, 2009-08-09, 2009-08-23, 2009-08-30, 2009-09-06, 2009-09-13, 2009-09-20, 2009-09-27, 2009-10-04, 2009-10-11, 2009-10-18, 2009-10-25, 2009-11-01, 2009-11-08, 2009-11-15, 2009-11-22, 2009-12-06, 2009-12-13, 2009-12-20, \textcolor{red}{2009-12-25}, 2009-12-27, 2010-01-03, 2010-01-10, 2010-01-17, 2010-01-24, 2010-01-31, 2010-02-07, 2010-02-21, 2010-02-28, \textcolor{red}{2010-03-03}, 2010-03-07, 2010-03-14, 2010-03-21, 2010-03-28, 2010-04-04, 2010-04-11, 2010-04-18, 2010-04-25, 2010-05-02, 2010-05-09, 2010-05-16, 2010-05-23, 2010-05-30, 2010-06-06, 2010-06-13, 2010-06-20, 2010-06-25, 2010-07-04, 2010-07-11, 2010-07-18, 2010-07-25, 2010-08-01, 2010-08-08, 2010-08-15, 2010-08-22, 2010-08-29, 2010-09-05, 2010-09-12, 2010-09-19, 2010-09-26, 2010-10-03, 2010-10-10, 2010-10-17, 2010-10-24, 2010-10-31, 2010-11-07, 2010-11-14, 2010-11-21, 2010-11-28, 2010-12-05, 2010-12-12, 2010-12-19, 2010-12-26, 2011-01-02, 2011-01-09, 2011-01-16, \textcolor{red}{2011-01-28}, 2011-01-30, 2011-02-06, 2011-02-13, 2011-02-20, 2011-02-27, 2011-03-06, 2011-03-13, 2011-03-20, 2011-03-27, 2011-04-03, 2011-04-10, 2011-04-17, 2011-04-24, 2011-05-01, \textcolor{red}{2011-05-05}, \textcolor{red}{2011-05-06}, 2011-05-08, 2011-05-15, 2011-05-22, 2011-05-29, 2011-06-05, 2011-06-12, 2011-06-19, 2011-06-26, 2011-07-03, 2011-07-10, 2011-07-17, 2011-07-24, 2011-08-07, 2011-08-14, 2011-08-21, 2011-08-28, 2011-09-04, 2011-09-11, 2011-09-18, 2011-09-25, 2011-10-02, 2011-10-09, 2011-10-16, 2011-10-23, 2011-10-30, 2011-11-06, 2011-11-13, 2011-11-20, 2011-11-27, 2011-12-04, 2011-12-11, 2011-12-18, 2011-12-25, 2012-01-08, 2012-01-15, 2012-01-17, 2012-01-22, 2012-01-29, 2012-01-31, 2012-02-05, 2012-02-12, 2012-02-19, 2012-02-26, \textcolor{red}{2012-02-29}, 2012-03-04, 2012-03-11, 2012-03-18, 2012-03-25, 2012-04-01, 2012-04-08, 2012-04-15, 2012-04-22, 2012-04-29, 2012-05-06, \textcolor{red}{2012-05-11}, 2012-05-13, 2012-05-20, 2012-05-27, 2012-07-01, 2012-07-08, 2012-07-15, 2012-07-22, 2012-07-29, 2012-11-18, 2012-12-25, \textcolor{red}{2013-01-01}, \textcolor{red}{2013-12-25}, \textcolor{red}{2014-01-01}, \textcolor{red}{2014-09-30}, 2014-11-09, \textcolor{red}{2014-11-26}, \textcolor{red}{2014-12-02}, \textcolor{red}{2014-12-03}, \textcolor{red}{2014-12-25}, \textcolor{red}{2015-01-01}, 2015-05-31, \textcolor{red}{2015-12-25}}\\[0.3cm]

The next Flg. (ref) is a plot of the daily returns for the data period used.

From the data, we computed the observed returns of equation ((ref)) every 30 seconds ($m=2880$), 1 minute ($m=1440$), 5 minutes ($m=288$), 10 minutes ($m=144$) and 30 minutes ($m=48$) multiplied by 100. Table (ref) is the basic statistics of the computed observation returns.

In all cases, the mean is almost equal to zero. The larger $m$ is, the smaller the mean and variance are. As for the autocorrelation, SAC(1) is relatively large, while those after SAC(2) are quite small. In addition, the larger $m$ is, the larger the first-order autocorrelation is. Therefore, for each autocorrelation at $m$, we performed the $k$-order Ljung-Box test with $k=1,2,\ldots,1000$, and all the autocorrelations are significant even at the significance level of 1$\%$. In addition, $k$-order autocorrelations were tested for each $k$. The SAC \textcolor{blue}{*} are the values that are significant at the 1$\%$ level of significance in this test. The next table(ref) is the autocorrelation $\Omega_i/\Omega_0$ calculated from the autocovariance $\Omega_i,\ i=0, 1, \ldots, q$ of MN calculated from the Theorem(ref).

For $\hat{q}+1$, the starting point of $k$ where SAC($k$) took three consecutive non-significant values in this test is considered as the truncation of the autocorrelation, and this $k-1$ is set as $\hat{q}+1$. With this $\hat{q}$, we determined the MA order $q$ of MN. The larger $m$ is, the larger $q$ is. The variance appears to decrease with $m$. With the exception of $m=2880$, the autocorrelation for each order is larger for larger $m$.

From these observed returns, we calculated each RV using the Formula ((ref)). The following figure shows a plot of the calculated RVs.

Table (ref) is the basic statistics for each calculated RV and RK.

As seen in Fig.(ref), we can confirm that the larger $m$ is, the larger the average of RV is.

Comparison method

Estimation errors are obtained by estimating the parameters of the model and IV using all the data. The 1-ahead forecasting error uses 1500 days of data, estimates the parameters of the model, and forecasts IV 1-ahead. This is done 500 times, shifting the data interval by one day. Since IV is unobservable, we use the Realized Kernel of Barndorff-Nielsen et al [2011] as a proxy variable:

Barndorff-Nielsen et al[2011] constructed an estimator Realized Kernel $RK_t^{(m)}$ with asymptotic normality in the presence of MNs satisfying HL.1, HL.2 and HL.3:

$K(x)$:$\mathbb{R}^1\to\mathbb{R}^1$ is a Kernel function, which satisfies appropriate non-negative definiteness and integrability in the $C^2$ class, $K(0)=1, K'(0)=0$. The parameter $H$ is the bandwidth, and when $H\propto m^{\frac{3}{5}}$, $RK_t^{(m)}$ is asymptotically normal with the convergence rate $m^{\frac{1}{5}}$. Barndorff-Nielsen et al [2011], $K(x)$ is a Parzen kernel:

The optimal parameters recommended are $H=c\xi^{\frac{4}{5}}m^{\frac{3}{5}}, c=3.5134,\xi^2=\cfrac{{\rm Var}[\varepsilon_t]}{\sqrt{T\int^T_0\sigma^4(t)dt}}$. Here, following Ubukata, Watanabe [2011], $\widehat{{\rm Var}[\varepsilon_t]}=\cfrac{RV_t^{(m)}}{2m}$ is estimated daily from $\frac{RV_t^{(m)}}{2m}$, $\int^T_0\sigma^4(t)dt$ is estimated daily from the Reaized Quarticity(Barndorff-Nielsen, Shephard [2002]): $\displaystyle RQ_t^{(m)}=\cfrac{m}{3}\sum_{i=1}^mr^{(m)4}_{t-1+\frac{i}{m}}$, estimated daily with $m=96$ and We use $RK_t^{(2880)}$, calculated from the observation returns every 30 seconds, as a proxy variable for $IV_t$.

The following tables(ref) and (ref) are the estimated SR-SARV model and MN parameters for each model for $p=1$ and $2$, respectively.

For each model, the parameters $\sigma^2, \omega_1^2, \lambda_1$ of the SR-SARV model have similar values regardless of $m$, which is consistent with the theory that the parameters of the SR-SARV model are independent of $m$.

Although there seems to be no regularity in the zero-f model with respect to $\Omega_0$, the NW model takes similar values regardless of $p, m$. This is consistent with the results of Nagakura and Watanabe [2009, 2015]. On the other hand, the weak-f model increases as $m$ decreases, and is estimated to be larger than the NW model. Furthermore, $\Omega_0$ is always estimated larger for the weak-f model than for the zero-f model. This suggests that the method of theorem ref{thm3} underestimates $\Omega_0$. Also, $f(m)$ are all negative, this implies a negative correlation between returns $r_t$ and MN $\varepsilon_t$.

The following Tables (ref) and (ref) are the estimation results of the State Space Form for each model with $p=1, 2$, respectively.

As in the SR-SARV model, the parameters $c_{IV}, \phi_1, \phi_2, \theta_1, \theta_2, \sigma_\eta^2$ of IV have similar values independent of $m$.

The bias and estimation error parameters $c_u, \theta_u, \sigma_\xi^2, \sigma_d^2$ seem to depend on $m$, but they take similar values for $p$. The $c_u$ of the NW model increases as $m$ increases. This is consistent with the results of the volatility signature plot (Fig (ref)). The $c_u$ of the zero-f and weak-f models also increases with the value of $m$ except for $m=2880$, while $\sigma_\xi^2$ conversely decreases. The $\theta_u$ of the NW and weak-f models are larger for larger $m$.

The $\sigma_d^2$ of the NW model, zero-f and weak-f models for $p=1$ increases as $m$ decreases. This is consistent with the asymptotics of RV. On the other hand, the same value is obtained for all $m$ for zero-f and weak-f models $p=2$. For the NW model and the zero-f and weak-f models for $p=1$, it is a natural interpretation that the variance of estimation error $\sigma_d^2$ increases as the sample size $m$ gets smaller. Meanwhile, it is interesting to note that the $sigma_d^2$ of the zero-f and weak-f models for $p=2$ take the same value of 1.527. This value is smaller than in any other case.

The following tables (ref), (ref) are the basic statistics of the estimated IV for each model for $p=1, 2$.

In the zero-f model, ${\rm Var}\left[u_t\right]$ is very small, and the result is like a parallel shift of NCRV by $c_u$, so that the auto-covariance is almost the same as that of the original NCRV. On the other hand, the values of the NW model and the weak-f model are similar, and their autocovariances after the first order are larger than those of the original NCRV. This is because the variance of NCRV is larger than that of IV due to the variation of $u_t$. ${\rm Var}[IV_t]/{\rm Var}[RV_t^*], {\rm Var}[u_t]/{\rm Var}[RV_t^*]$ is the ratio of biases due to IV and MN in NCRV. For $p=1$, the ratio of IV between the NW model and the weak-f model is less than half, and decreases with $m$ except for $m=1440$. For $p=2$, the ratio of IVs of the NW model and the weak-f model behaves similarly to that for $p=1$, but is slightly larger. For the zero-f model, the ratio of IV is almost unity. In other words, the variation of bias due to MN autocorrelation is quite small. The ${\rm Corr}[r_t, \varepsilon_t]$ is obtained from the parameter estimates of the weak-f model by the Equation ((ref)) ${\rm Corr}[r_t, \varepsilon_t]=\cfrac{f(m){\rm E}[\sigma(t)]}{\sqrt{m\ sigma^2\left(\cfrac{(f(m))^2}{m}+\Omega_0\right)}}$ is the correlation between returns and MN. The absolute value of $p=2$ is larger than $p=1$. Furthermore, the absolute value decreases as $m$ increases except for $m=2880$.

The following figure(ref) shows the $\widehat{IV}_t$ estimation results for each model at each $m$ (top, solid line: $RV_t^{*(m)}$, dashed line: $\widehat{IV}_t$), $\hat{u}_t^{(m)}$ (middle), $\widetilde{IV}_t$ one period ahead (bottom).

The following table shows the IV estimation and one period ahead forecasting errors for each model with $IV_t=RK_t ^{(2880)}$ for each model. The best value is marked with \textcolor{red}{*}, the second and third best values are marked with \textcolor{red}{**} and \textcolor{red}{***}, respectively.

From the table (ref), (ref), (ref), (ref), we know the following.

For IV estimation, zero-f(2) seems to be the most dominant.

From the tables (ref), (ref), (ref), and (ref), we know the following.

For IV one period ahead forecasts, weak-f(2) seems to be superior.

The zero-f model has by far the best score in IV estimation, but the one-period ahead forecasting is quite poor. The weak-f model is superior in one period ahead prediction, and the score of the weak-f model is not bad in IV estimation. The NW model has a score similar to that of the weak-f model in IV estimation, but is clearly inferior in one period ahead prediction.

Conclusion

In this paper, we relaxed the assumptions of Nagakura and Watanabe[2009, 2015], and extended their model assuming Hansen and Lunde[2006] type dependent MN that auto-correlate and correlate with return. Under these assumptions, we derived the properties of the bias due to MN $u_t^{(m)}$. Based on these, we proposed the parameters of $u_t^{m}$ and dependent MN estimation method from NCRV and observation return data.

In model comparison by simulation, proposal model was superior to the existing method in the accuracy of parameter estimation, IV estimation and 1-ahead prediction.

In model comparison by actual data on the dollar-yen rate, we newly found the following