# Modeling Multiple Irregularly Spaced Financial Time Series

Chiranjit Dutta Corresponding Author<EMAIL_ADDRESS>Nalini
Ravishanker  Sumanta Basu  Department of Statistics and Data Science, Cornell
University, Ithaca, NY, USA

###### Abstract

In this paper we propose univariate volatility models for irregularly spaced
financial time series by modifying the regularly spaced stochastic volatility
models. We also extend this approach to propose multivariate stochastic
volatility (MSV) models for multiple irregularly spaced time series by
modifying the MSV model that was used with daily data. We use these proposed
models for modeling intraday logarithmic returns from health sector stocks
data obtained from Trade and Quotes (TAQ) database at Wharton Research Data
Services (WRDS).

##  1 Introduction

In finance, we often encounter time-series data that exhibits irregular
spacing in time meaning that the time intervals between successive data points
are not the same. This type of data is known as irregularly spaced time series
data. We propose univariate and multivariate stochastic volatility models for
irregularly spaced time series. Specifically, we build on the gap time
modeling idea of Nagaraja et al. (2011) to construct useful time series models
that can help better understand volatility patterns in irregularly spaced
financial time series. To achieve this, we modify existing stochastic
volatility models that were originally designed for regularly spaced data.
Additionally, we extend this approach to model multiple irregularly spaced
time series by modifying multivariate stochastic volatility (MSV) model of
(Chib et al., 2009) that was used with daily data. High-frequency financial
data are inherently irregularly spaced since trades can occur at any point in
time. Moreover, the microstructure of financial markets, such as the methods
of placing and executing orders, can also cause unevenly spaced data. For
instance, in electronic markets, trades can be executed rapidly, and traders
may use different strategies and algorithms that result in varying transaction
frequencies. The literature on dynamic statistical modeling is sparse.
Practitioners interested in understanding dynamic evolution of stock
properties require specialized models since the ones designed for regularly
spaced data are inappropriate. One of the early approaches to model the return
series sampled at irregularly spaced time intervals set by the trade arrivals
is ACD-GARCH model by Ghysels and Jasiak (1998). This is a random coefficient
GARCH, or doubly stochastic GARCH, where the stochastic durations between
transactions determine the parameter dynamics for the GARCH equation. They
proposed a Generalized Method of Moments (GMM) based two-step procedure for
estimating the parameters. ACD-GARCH is quite cumbersome to estimate and also
difficult to generalize to multiple time series. Meddahi et al. (2006)
proposed a GARCH type model for irregularly spaced data which is an exact
discretization of continuous time stochastic volatility processes observed at
irregularly spaced times. Their model combines the advantages of ACD-GARCH
Ghysels and Jasiak (1998) and ACD (Engle and Russell, 1998). A continuous
version of GARCH (COGARCH) is another way to model irregularly spaced time
series data (Maller et al., 2008). Recently, Buccheri et al. (2021a) propose
to model intraday log-prices through a multivariate local-level model with
score-driven covariance matrices and to treat asynchronicity as a missing
value problem.

More details are needed

Write a short description about each section?

##  2 Volatility Models for Univariate Irregularly Spaced Time Series

Let $\\{r_{t_{j}}\\}$ be a sequence of log-returns of a financial asset, where
$t_{j}$ denotes the time of the $j^{th}$transaction. Let and
$g_{j}=t_{j}-t_{j-1},j>1$ be the known gap times between consecutive returns.
Unlike regularly spaced time series (such as daily returns), for irregularly
spaced time series (such as transaction level intra-day returns), the gap
times $g_{j},~{}j>1$ are not the same.

###  2.1 Irregularly spaced stochastic volatility (IR-SV) Model

We define an irregular stochastic volatility (IR-SV) model for $r_{t_{j}}$ as

$\displaystyle r_{t_{j}}$
$\displaystyle=\exp\Bigg{(}\frac{h_{t_{j}}}{2}\Bigg{)}\epsilon_{t_{j}},\quad\epsilon_{t_{j}}\sim
N(0,1),$ (1) $\displaystyle h_{t_{1}}$
$\displaystyle=\mu+\eta_{t_{1}},\quad\eta_{t_{1}}\sim
N\Big{(}0,\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)},$ (2) $\displaystyle
h_{t_{j}}$
$\displaystyle=\mu+\phi^{g_{j}}(h_{t_{j-1}}-\mu)+\eta_{t_{j}},\quad\eta_{t_{j}}\sim
N\Bigg{(}0,\dfrac{\sigma_{\eta}^{2}(1-\phi^{2g_{j}})}{1-\phi^{2}}\Bigg{)},\text{
for $j>1$},$ (3)

where $\phi$ is the persistence parameter and $\mu$ is the location parameter.
In (2)-(3), we assume that the log-volatility process $\\{h_{t_{j}}\\}$ is a
stationary Gaussian autoregressive (AR) process with $|\phi|<1$, and $\phi\neq
0$. We also assume that the gap times $g_{j},~{}j>1$ are bounded away from
$0$.

The formulation of the IR-SV model in (1)-(3) can be viewed as an extension to
volatility modeling of the stationary gap time autoregressive (gap-AR) model
for an irregularly spaced time series described in Nagaraja et al. (2011).

######  Proposition 2.1.

$h_{t_{j}}$ is a weakly stationary process and the unconditional distribution
of $h_{t_{j}}$ is

$h_{t_{j}}\sim\text{N}\Big{(}\mu,\frac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)},$

###### Proof.

Let $x_{t_{j}}=h_{t_{j}}-\mu,\forall j\in\mathbb{Z}^{+}$. From (1)-(3),
$x_{t_{j}}$ can be recursively written as

$x_{t_{j}}=\sum_{k=1}^{j-1}\phi^{\sum_{i=k+1}^{j}g_{i}}\eta_{t_{k}}+\eta_{t_{j}}$
(4)

Hence it follows that $E(x_{t_{j}})=0,\forall j\in\mathbb{Z}^{+}$. For
$l=0,1,\cdots$, the gap time $l=t_{j+l}-t_{j}=\sum_{i=j+1}^{j+l}g_{i}$ and the
covariance of $x_{t_{j}}$ and $x_{t_{j+l}}$ is

$\displaystyle Cov(x_{t_{j}},x_{t_{j+l}})$
$\displaystyle=\text{E}\left[\Bigg{(}\sum_{k=1}^{j-1}\phi^{\sum_{i=k+1}^{j}g_{i}}\eta_{t_{k}}+\eta_{t_{j}}\Bigg{)}\Bigg{(}\sum_{k=1}^{j+l-1}\phi^{\sum_{i=k+1}^{j+l}g_{i}}\eta_{t_{k}}+\eta_{t_{j+l}}\Bigg{)}\right]$
$\displaystyle=\sum_{k=1}^{j-1}\phi^{\sum_{i=k+1}^{j}g_{i}+\sum_{i=k+1}^{j+l}g_{i}}\text{E}\left[\eta^{2}_{t_{k}}\right]+\text{E}\left[\eta^{2}_{t_{j}}\right]$
$\displaystyle=\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}\left[\phi^{2\sum_{i=2}^{j}g_{i}+l}+\sum_{k=2}^{j-1}\phi^{2\sum_{i=k+1}^{j}g_{i}+l}(1-\phi^{2g_{k}})+\phi^{l}(1-\phi^{2g_{j}})\right]$
$\displaystyle=\dfrac{\sigma_{\eta}^{2}\phi^{l}}{1-\phi^{2}}\left[\phi^{2\sum_{i=2}^{j}g_{i}}-\phi^{2\sum_{i=2}^{j}g_{i}}+\phi^{2\sum_{i=3}^{j}g_{i}}+\cdots+\phi^{2g_{j}}-\phi^{2g_{j}}+1\right]$
$\displaystyle=\dfrac{\sigma_{\eta}^{2}\phi^{l}}{1-\phi^{2}}$

Thus $Cov(x_{t_{j}},x_{t_{j+l}})$ is a function of gap time $l$ only and the
series $x_{t_{j}}$ and hence $h_{t_{j}}$ is a weakly stationary process. By
stationarity of $h_{t_{j}}$ and (3),

$\displaystyle\text{Var}(h_{t_{j}})=\phi^{2g_{j}}\text{Var}(h_{t_{j}})+\text{Var}(\eta_{t_{j}})$
$\displaystyle\implies\text{Var}(h_{t_{j}})=\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}$

Also, $\text{E}(h_{tj})=\mu$ and hence the unconditional distribution of
$h_{t_{j}}$ is
$\text{N}\Big{(}\mu,\frac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}$. ∎

######  Proposition 2.2.

We discuss properties of the distribution of $\\{r_{t_{j}}\\}$.

  1. 1.

The expectation and variance of the squared returns is

$\displaystyle\text{E}\big{(}r_{t_{j}}^{2}\big{)}=\exp\Big{(}\mu+\frac{\sigma_{\eta}^{2}}{2(1-\phi^{2})}\Big{)}$
(5)
$\displaystyle\text{Var}\big{(}r_{t_{j}}^{2}\big{)}=\exp\Big{(}2\mu+\dfrac{\sigma^{2}_{\eta}}{1-\phi^{2}}\Big{)}\Bigg{(}3\exp\Big{(}\dfrac{\sigma^{2}_{\eta}}{1-\phi^{2}}\Big{)}-1\Bigg{)}$
(6)

  2. 2.

Let the gap time $l=t_{j+l}-t_{j}=\sum_{i=j+1}^{j+l}g_{i}$, $l=0,1,\cdots$,
the autocovariance function of $r_{t_{j}}^{2}$ and $r_{t_{j+l}}^{2}$ is

$\displaystyle\text{Cov}\big{(}r_{t_{j}}^{2},r_{t_{j+l}}^{2}\big{)}=\exp\Big{(}2\mu+\frac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}\left[\exp\Big{(}\dfrac{\sigma_{\eta}^{2}\phi^{l}}{1-\phi^{2}}\Big{)}-1\right]$
(7)

  3. 3.

The kurtosis of the returns is

$\text{K}\big{(}r_{t_{j}}\big{)}=3\exp\Big{(}\dfrac{\sigma^{2}_{\eta}}{1-\phi^{2}}\Big{)}$
(9)

  4. 4.

The sequence of squared returns $r^{2}_{t_{j}}$ is a stationary process.

###### Proof.

The first two even moments of $r_{t_{j}}$ are

$\displaystyle\text{E}\big{(}r_{t_{j}}^{2}\big{)}=\text{E}\big{(}\exp(h_{t_{j}})\epsilon_{t_{j}}^{2}\big{)}=\text{E}\big{(}\exp(h_{t_{j}})\big{)}\text{E}\big{(}\epsilon_{t_{j}}^{2}\big{)}=\exp\Big{(}\mu+\frac{\sigma_{\eta}^{2}}{2(1-\phi^{2})}\Big{)},$
$\displaystyle\text{E}\big{(}r_{t_{j}}^{4}\big{)}=\text{E}\big{(}\exp(2h_{t_{j}})\epsilon_{t_{j}}^{4}\big{)}=\text{E}\big{(}\exp(2h_{t_{j}}\big{)}\text{E}\big{(}\epsilon_{t_{j}}^{4}\big{)}=3\exp\Big{(}2\mu+\frac{2\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}.$

Then,

$\displaystyle\text{Var}\big{(}r_{t_{j}}^{2}\big{)}=$ $\displaystyle
3\exp\Big{(}\mu+\frac{3\sigma_{\eta}^{2}}{2(1-\phi^{2})}\Big{)}-\exp\Big{(}2\mu+\frac{\sigma_{\eta}^{2}}{(1-\phi^{2})}\Big{)}$
$\displaystyle=$
$\displaystyle\exp\Big{(}2\mu+\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}\Big{(}3\exp\Big{(}\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}-1\Big{)},$

and the kurtosis is

$\text{K}\big{(}r_{t_{j}}\big{)}=\dfrac{\text{E}\big{(}r_{t_{j}}^{4}\big{)}}{\text{E}\big{(}r_{t_{j}}^{2}\big{)}^{2}}=3\exp\Big{(}\frac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}.$

The kurtosis is always greater than $3$ as long as $\sigma_{\eta}^{2}>0$. For
$l=0,1,\cdots$, the autocovariance function of $r_{t_{j}}^{2}$ is

$\displaystyle\text{Cov}(r^{2}_{t_{j+l}},r^{2}_{t_{j}})$
$\displaystyle=\text{E}(r^{2}_{t_{j+l}}r^{2}_{t_{j}})-\text{E}(r^{2}_{t_{j+l}})\text{E}(r^{2}_{t_{j}})$
$\displaystyle=\text{E}(\exp(h_{t_{j}}+h_{t_{j+l}})\epsilon_{t_{j}}^{2}\epsilon_{t_{j+l}}^{2})-\exp\Big{(}2\mu+\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}$

By independence of $\epsilon_{t_{j}}$’s and $h_{t_{j}}$’s, we have

$\displaystyle\text{Cov}(r^{2}_{t_{j+l}},r^{2}_{t_{j}})$
$\displaystyle=\text{E}(\exp(h_{t_{j}}+h_{t_{j+l}}))-\exp\Big{(}2\mu+\dfrac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}.$
(10)

We note that

$\displaystyle\text{Var}(h_{t_{j}}+h_{t_{j+l}})$
$\displaystyle=\text{Var}(h_{t_{j}})+\text{Var}(h_{t_{j+l}})+2\text{Cov}(h_{t_{j}},h_{t_{j+l}})$
$\displaystyle=\dfrac{2\sigma_{\eta}^{2}(1+\phi^{l})}{1-\phi^{2}}$

By normality of $h_{t_{j}}$, it follows that $h_{t_{j}}+h_{t_{j+l}}\sim
N\Big{(}2\mu,\dfrac{2\sigma_{\eta}^{2}(1+\phi^{l})}{1-\phi^{2}}\Big{)}$.
Therefore (10) becomes

$\displaystyle\text{Cov}(r^{2}_{t_{j+l}},r^{2}_{t_{j}})$
$\displaystyle=\exp\Big{(}2\mu+\frac{\sigma_{\eta}^{2}}{1-\phi^{2}}\Big{)}\left[\exp\Big{(}\dfrac{\sigma_{\eta}^{2}\phi^{l}}{1-\phi^{2}}\Big{)}-1\right].$
(11)

Therefore, $r_{t_{j}}^{2}$ is a stationary process. ∎

###  2.2 Bayesian Estimation of the IR-SV Model

Let $\boldsymbol{r}=(r_{t_{1}},\cdots,r_{t_{T}})$ be the set of irregularly
spaced observed returns following the IR-SV model in (1)-(3). Let
$\boldsymbol{\theta}=(\sigma_{\eta}^{2},\phi,\mu)$ be the set of model
parameters. For $j=1,\cdots,T$, we can rewrite (1)-(3) as follows:

$\displaystyle r_{t_{j}}|h_{t_{j}},\mathbf{\theta}$ $\displaystyle\sim
N\big{(}0,\exp(h_{t_{j}})\big{)}$ (12) $\displaystyle h_{t_{1}}$
$\displaystyle\sim N\big{(}\mu,\dfrac{\sigma^{2}_{\eta}}{1-\phi^{2}}\big{)}$
(13) $\displaystyle h_{t_{j}}|h_{t_{j-1}}$ $\displaystyle\sim
N\Big{(}\mu+\phi^{g_{j}}(h_{t_{j-1}}-\mu),\dfrac{\sigma^{2}_{\eta}(1-\phi^{2g_{j}})}{1-\phi^{2}}\Big{)},\quad
j>1.$ (14)

We assume the following priors:

$\displaystyle\dfrac{\phi+1}{2}\sim B(20,1.5),$
$\displaystyle\dfrac{1}{\sigma_{\eta}^{2}}\sim G(2.5,0.025),$
$\displaystyle\mu\sim N(0,10),$

where $B$ denotes the beta distribution, $G$ denotes the gamma distribution
and $N$ is the normal distribution. Let
$p(\boldsymbol{\theta})=p(\phi)p(\sigma^{2}_{\eta})p(\mu)$ denote the joint
prior of all the parameters, assuming independence. The joint posterior
distribution of the hidden volatility states
$\boldsymbol{h}=(h_{t_{1}},\cdots,h_{t_{T}})$ and the parameters
$\boldsymbol{\theta}=(\sigma_{\eta}^{2},\phi,\mu)$ is obtained by Bayes’ rule
as

$p(\boldsymbol{h},\boldsymbol{\theta}|\boldsymbol{r})\propto
p(\boldsymbol{\theta})\prod_{j=1}^{T}p(r_{t_{j}}|h_{t_{j}},\boldsymbol{\theta})p(h_{t_{j}}|h_{t_{j-1}},\boldsymbol{\theta}),$
(15)

where $p(r_{t_{j}}|h_{t_{j}},\boldsymbol{\theta})$ and
$p(h_{t_{j}}|h_{t_{j-1}},\boldsymbol{\theta})$ follows from (12)-(14). We
employ a Metropolis-Hastings adaptive random-walk sampler for each of the
parameters with a univariate normal proposal distribution within R package
NIMBLE package (de Valpine et al., 2017) to generate samples from their
respective posterior distribution.

###  2.3 Simulation Study

We demonstrate the accuracy of the Bayesian estimation of the IR-SV model
parameters under a few different simulation setups. We generated $100$ sets
(replicates) of zero mean log returns data, each of length $T=5000$. The gap
times $g_{j},1\leq j\leq T$ are generated as follows:

  * •

Generate $g_{j}\sim\mathcal{P}(\lambda=3)$ and $g_{j}>0,\forall\ j$ where
$\mathcal{P}(\lambda)$ is the Poisson distribution with mean $\lambda$.

  * •

Scale the $g_{j}$’s such that $0<g_{j}\leq 1,\forall j$

In our empirical analysis we have observed the mean gap times to be around 3
seconds with $10\%$ of them greater than 5 seconds and since the persistence
parameter $|\phi|<1$, we scale the gap times such that $g_{j}\in(0,1]$ to
ensure $\phi^{g_{j}}$ is bounded away from 0, $\forall j$. To study the effect
of persistence parameter $\phi$ on the estimation accuracy, in all three
scenarios we generated zero mean log returns from the IR-SV model in (1-3)
with the true parameter values of $\mu=-9$, $\sigma=0.8$. In scenario 1, the
true value of the persistence parameter $\phi$ is 0.2, in scenario 2, true
value of $\phi$ is 0.6 and in scenario 3, true value of $\phi$ is 0.9.

Parameter | True Value | Mean | Q(2.5%) | Q(97.5%)
---|---|---|---|---
$\mu$ | -9.0000 | -8.9997 | -9.1117 | -8.8627
$\phi$ | 0.2000 | 0.2438 | 0.1565 | 0.3346
$\sigma$ | 0.8000 | 0.7918 | 0.7366 | 0.8551
Table 1: True values and posterior estimates of parameters for Scenario 1 from the IR-SV model. Parameter | True Value | Mean | Q(2.5%) | Q(97.5%)
---|---|---|---|---
$\mu$ | -9.0000 | -8.9833 | -9.2253 | -8.7326
$\phi$ | 0.6000 | 0.6253 | 0.5515 | 0.7056
$\sigma$ | 0.8000 | 0.7832 | 0.7026 | 0.8769
Table 2: True values and posterior estimates of parameters for Scenario 2 from the IR-SV model. Parameter | True Value | Mean | Q(2.5%) | Q(97.5%)
---|---|---|---|---
$\mu$ | -9.0000 | -8.9872 | -9.9697 | -8.0021
$\phi$ | 0.9000 | 0.8960 | 0.8365 | 0.9458
$\sigma$ | 0.8000 | 0.7860 | 0.7059 | 0.8818
Table 3: True values and posterior estimates of parameters for Scenario 3 from
the IR-SV model.

We run 520,000 MCMC iterations and discard the first 20,000 as burn-in and
thinned every $1000^{th}$ sample to reduce autocorrelation between MCMC
samples. Convergence for the parameters are assessed using trace and posterior
density plots, autocorrelation and thinning (not shown here). In Table 1,
Table 2, and Table 3, we report the posterior sample means of the parameters
along with their true values and $95\%$ credible intervals, averaged the 100
data sets (replicates). The true values of the parameters lie inside the
$95\%$ credible intervals for all parameters in all three scenarios.

##  3 Multivariate Stochastic Volatility Models for Irregularly Spaced
Financial Time Series (IR-MSV)

We propose a multivariate stochastic volatility model to fit irregularly
spaced synchronized intraday zero mean log returns for multiple assets by
modifying the basic multivariate stochastic volatility (BMSV) model (Chib et
al., 2009) and using the gap time idea of Nagaraja et al. (2011) for handling
the latent state.

###  3.1 Refresh Time Sampling

Why synchronization is necessary? In this section we describe refresh time
sampling which allows synchronization of high frequency prices from multiple
stocks. In high frequency trading (HFT) conducting multivariate analysis is
difficult since assets do not trade on a fixed grid, trades and quotes don’t
arrive synchronously. Hence synchronization of HFT data from multiple assets
is necessary. As in Section (2.5.3), we briefly describe the refresh sampling
procedure for synchronization of high frequency prices from multiple stocks.

Assume there are $p$ stocks, and trading time of the $i^{th}$ stock is given
by $t_{i\ell}$, $\ell=1,\cdots,n_{i}$ , $i=1,\cdots,p$. For a given time $t$,
define $N_{t}^{i}=$ the number of $t_{i\ell}\leq t$, $\ell=1,\cdots,n_{i}$,
which counts the number of distinct data points $t_{i\ell}$ available for the
$i^{th}$ asset up to time $t$. The first refresh time is defined as
$\tau_{1}=\max\left\\{t_{11},\ldots,t_{p1}\right\\}$, which is the first time
taken to trade all assets and refresh their posted prices. The subsequent
refresh times are defined as follows. Given the $j^{th}$ refresh time
$\tau_{j}$, define the $(j+1)^{th}$ refresh time

$\tau_{j+1}=\max\left\\{t_{1,N_{\tau_{j}}^{1}+1},\ldots,t_{p,N_{\tau_{j}}^{p}+1}^{p}\right\\}$
(16)

Suppose there are $m$ refresh time points $\tau_{1},\cdots,\tau_{m}$.
Intuitively, $\tau_{2}$ is the second time when all the assets are traded and
their prices are refreshed. In our empirical analysis we consider three health
sector stocks BMY, CVS and MDT traded on NYSE on $24^{th}$ June, $2016$. We
consider the high frequency prices for these three stocks and illustrate the
refresh time sampling procedure for synchronization of multivariate analysis
which we discuss later. Figure 1 illustrates the refresh time sampling idea
and in this example we have $\tau_{1}=09:45:01.00$,$\tau_{2}=09:45:02.00$,
$\tau_{3}=09:45:03.00$, $\tau_{4}=09:45:06.00$ and $\tau_{5}=09:45:07.00$ are
the first five refresh times.

Figure 1:  Refresh Time Sampling for BMY, CVS and MDT. The vertical dotted
blue lines represent the refresh time points on $24^{th}$ June, 2016

###  3.2 IR-MSV Model Formulation

Let
$\boldsymbol{r}_{t_{j}}=(r_{1t_{j}},r_{2t_{j}},\cdots,r_{pt_{j}})^{\prime}$,
$j=1,\cdots,T$ be a $p$-dimensional vector of zero mean log returns of $p$
financial assets observed at irregularly spaced time points $t_{1},\cdots
t_{T}$. Here, $t_{j}$ denotes the time of the $j^{th}$ observation, and
$g_{j}=t_{j}-t_{j-1},j>1$ denotes the known gap time between consecutive
observations. We introduce irregular multivariate stochastic volatility (IR-
MSV) model as

$\boldsymbol{r}_{t_{j}}=\mathbf{H}_{t_{j}}^{1/2}\boldsymbol{\epsilon_{t_{j}}},\quad\boldsymbol{\epsilon_{t_{j}}}\sim\text{MVN}(\boldsymbol{0},\mathbf{R}),$
(17)

where,

$\displaystyle\mathbf{H}_{t_{j}}=\text{diag}\big{(}\exp(h_{1,t_{j}}),\cdots,\exp(h_{p,t_{j}})\big{)},\quad
j=1,\cdots,T,$ (18) $\displaystyle
h_{i,t_{1}}=\mu_{i}+\eta_{i,t_{1}},\quad\eta_{i,t_{1}}\sim\text{N}\bigg{(}0,\dfrac{\sigma_{i}^{2}}{1-\phi_{i}^{2}}\bigg{)},\quad
i=1,\cdots,p,$ $\displaystyle
h_{i,t_{j}}=\mu_{i}+\phi_{i}^{g_{j}}(h_{i,t_{j-1}}-\mu_{i})+\eta_{i,t_{j}},\quad\eta_{i,t_{j}}\sim\text{N}\bigg{(}0,\dfrac{\sigma_{i}^{2}(1-\phi_{i}^{2g_{j}})}{1-\phi_{i}^{2}}\bigg{)},\quad
i=1,\cdots,p;j=1,\cdots,T,$

where $\phi_{i}$ is the persistence parameter constrained by $|\phi_{i}|<1$,
$i=1,\cdots,p$, $\mu_{i}$ is the location parameter for each asset
$i=1,\cdots,p$ and $\mathbf{R}$ is the correlation matrix of the observation
errors. Following Nagaraja et al. (2011) and our discussion in the previous
section, we observe that $\\{h_{i,t_{j}}\\}$ is a stationary process for each
$i=1,\cdots,p$. Since the error correlation matrix $\mathbf{R}$ is assumed to
remain constant over time, we refer to (17) as a constant correlation IR-MSV
model.

######  Proposition 3.1.

Let
$\boldsymbol{r}^{2}_{t_{j}}=(r^{2}_{1t_{j}},r^{2}_{2t_{j}},\cdots,r^{2}_{pt_{j}})^{\prime}$,
$j=1,\cdots,T$ be the vector of squared returns. We discuss properties of the
distribution of $\\{\boldsymbol{r}_{t_{j}}^{2}\\}$.

  1. 1.

Let
$\text{E}(\boldsymbol{r}^{2}_{t_{j}})=\boldsymbol{m}=(m_{1},\cdots,m_{p})^{\prime}$
be the $p\times 1$ vector of expected values where $m_{i}$ is

$\displaystyle
m_{i}=\text{E}\big{(}r_{i,t_{j}}^{2}\big{)}=\exp\Big{(}\mu_{i}+\frac{\sigma_{i}^{2}}{2(1-\phi_{i}^{2})}\Big{)},\quad
i=1,\cdots,p$ (19)

  2. 2.

Let $\text{Cov}(\boldsymbol{r}^{2}_{t_{j}})=\boldsymbol{\Sigma}$ be the
$p\times p$ variance-covariance matrix of squared returns vector. The
$(i,k)^{th}$ element of $\boldsymbol{\Sigma}$ is

$\text{Cov}(r_{i,t_{j}}^{2},r_{k,t_{j}}^{2})=\boldsymbol{\Sigma}_{ik}=\begin{cases}\exp\Big{(}2\mu_{i}+\dfrac{\sigma^{2}_{i}}{1-\phi_{i}^{2}}\Big{)}\Bigg{(}3\exp\Big{(}\dfrac{\sigma^{2}_{i}}{1-\phi_{i}^{2}}\Big{)}-1\Bigg{)}&i=k\\\
2\rho_{ik}^{2}\exp\left[(\mu_{i}+\mu_{k})+\dfrac{1}{2}\Big{(}\dfrac{\sigma_{i}^{2}}{1-\phi_{i}^{2}}+\dfrac{\sigma_{k}^{2}}{1-\phi_{k}^{2}}\Big{)}\right]&i\neq
k\end{cases}$ (21)

###### Proof.

∎

###  3.3 Bayesian Inference

In this section, we describe Bayesian analysis for fitting the IR-MSV model in
(17) to irregularly spaced multiple financial time series. We show the
likelihood, prior, posterior and also discuss the computational aspects using
MCMC algorithms. The likelihood function is

$\mathcal{L}(\boldsymbol{\Theta}|\boldsymbol{r}_{t_{1}},\cdots,\boldsymbol{r}_{t_{T}})=|\mathbf{\Sigma}_{\boldsymbol{r}}|^{-T/2}\prod_{j=1}^{T}\exp\Bigg{(}\dfrac{-\boldsymbol{r}_{t_{j}}\mathbf{\Sigma}_{\boldsymbol{r}}^{-1}\boldsymbol{r}_{t_{j}}^{\prime}}{2}\Bigg{)},$
(22)

where
$\mathbf{\Sigma}_{\mathbf{r}}=\mathbf{H}_{t_{j}}^{1/2}\mathbf{R}\mathbf{H}_{t_{j}}^{1/2}$
, $\boldsymbol{\sigma}=(\sigma_{1},\cdots,\sigma_{p})$ and
$\boldsymbol{\phi}=(\phi_{1},\cdots,\phi_{p})$, and
$\boldsymbol{\Theta}=(\boldsymbol{\sigma},\boldsymbol{\phi},\mathbf{R},\mu)$.

For $i=1,\cdots,p$, we assume the following priors:

$\displaystyle\mu_{i}\sim\text{N}(0,10),$
$\displaystyle\dfrac{1}{\sigma_{i}^{2}}\sim\text{G}(2.5,0.025),$
$\displaystyle\mathbf{R}\sim\text{LKJ}(\eta=1.2),\text{ where
}\text{LKJ}(\mathbf{\Sigma}|\eta)\propto\det(\mathbf{\Sigma})^{(\eta-1)},$
$\displaystyle\phi_{i}\sim\text{N}(0,0.5),$

where Lewandowski-Kurowicka-Joe (LKJ) distribution Lewandowski et al. (2009)
is a useful prior for correlation matrices, and IG denotes the inverse gamma
distribution.

The posterior distribution of $\boldsymbol{\Theta}$ is

$\displaystyle\pi(\boldsymbol{\Theta}|\boldsymbol{r}_{t_{1}},\cdots,\boldsymbol{r}_{t_{T}})\propto$
$\displaystyle\quad|\mathbf{\Sigma}_{\boldsymbol{r}}|^{-T/2}\prod_{j=1}^{T}\exp\Bigg{(}\dfrac{-\boldsymbol{r}_{t_{j}}\mathbf{\Sigma}_{\boldsymbol{r}}^{-1}\boldsymbol{r}_{t_{j}}^{\prime}}{2}\Bigg{)}$
(23)
$\displaystyle\times\prod_{j=1}^{T}\pi(\boldsymbol{h}_{t_{j}}|\boldsymbol{h}_{t_{j-1}},\boldsymbol{\sigma},\boldsymbol{\phi})\times\pi(\boldsymbol{\sigma})\pi(\boldsymbol{\phi})\pi(\mathbf{R})\pi(\mu),$

where we have assumed independence of the priors and
$\pi(\boldsymbol{h}_{t_{j}}|\boldsymbol{h}_{t_{j-1}},\boldsymbol{\sigma},\boldsymbol{\phi})$
is the joint distribution of the stochastic volatilities as in (18). We employ
a block random walk sampler for the correlation matrix $\mathbf{R}$ and for
other parameters i.e. $(\boldsymbol{\sigma},\boldsymbol{\phi},\mu)$ we employ
a Metropolis-Hastings adaptive random-walk sampler with a univariate normal
proposal distribution, which are the default samplers in NIMBLE package (de
Valpine et al., 2017) in R.

###  3.4 Simulation Study

In this simulation our goal is to examine the effect of the correlation among
the components on the accuracy of estimation. We generated 100 sets of zero
mean log returns data from the IR-MSV model represented by equations (17) and
(18) with $p=3$ and $T=5000$. The gap times $g_{j},1\leq j\leq T$ are
generated using a similar procedure as described in Section 2.3. The $3\times
3$ correlation matrix has the following representation

$\mathbf{R}=\begin{pmatrix}1&\rho_{12}&\rho_{13}\\\ \rho_{12}&1&\rho_{23}\\\
\rho_{13}&\rho_{23}&1\end{pmatrix}$

We describe three different scenarios. In scenario 1, we consider moderate
positive correlations among the components with $\rho_{12}=0.6$,
$\rho_{13}=0.4$ and $\rho_{23}=0.2$. In scenario 2, we consider a mix of high
and low positive correlation and a negative correlation with $\rho_{12}=-0.4$,
$\rho_{13}=0.7$ and $\rho_{23}=0.3$. In scenario 3, we consider all
correlations to be equal, high and positive with
$\rho_{12}=\rho_{13}=\rho_{23}=0.7$.

We run 210000 MCMC iterations and discard the first 10000 as burn-in and
thinned every $400^{th}$ sample to reduce autocorrelation between MCMC
samples. We assessed the convergence for the parameters using trace and
posterior density plots. In Table 4, Table 5, and Table 6 we report the
posterior sample mean along with their true values and $95\%$ credible
intervals averaged the 100 data sets for scenario 1, scenario 2 and scenario 3
respectively.

Parameter | True Value | Mean | Q(2.5%) | Q(97.5%)
---|---|---|---|---
$\rho_{12}$ | 0.6000 | 0.5987 | 0.5740 | 0.6158
$\rho_{13}$ | 0.4000 | 0.3987 | 0.3724 | 0.4286
$\rho_{23}$ | 0.2000 | 0.2008 | 0.1740 | 0.2254
$\sigma^{2}_{1}$ | 1.0000 | 0.9831 | 0.8639 | 1.1203
$\sigma^{2}_{1}$ | 0.8000 | 0.7807 | 0.6445 | 0.9119
$\sigma^{2}_{3}$ | 0.5000 | 0.5658 | 0.4606 | 0.6636
$\mu_{1}$ | -9.0000 | -8.9964 | -9.3681 | -8.6911
$\mu_{2}$ | -9.5000 | -9.4927 | -9.6809 | -9.3255
$\mu_{3}$ | -8.5000 | -8.5018 | -8.5973 | -8.3994
$\phi_{1}$ | 0.7000 | 0.7021 | 0.6374 | 0.7694
$\phi_{2}$ | 0.5000 | 0.4945 | 0.4119 | 0.5651
$\phi_{3}$ | 0.3000 | 0.3009 | 0.2181 | 0.3943
Table 4: True values and posterior estimates of parameters for Scenario 1 from the IR-MSV model. Parameter | True Value | Mean | Q(2.5%) | Q(97.5%)
---|---|---|---|---
$\rho_{12}$ | -0.4000 | -0.3993 | -0.4245 | -0.3749
$\rho_{13}$ | 0.7000 | 0.7003 | 0.6835 | 0.7142
$\rho_{23}$ | 0.3000 | 0.3003 | 0.2710 | 0.3268
$\sigma^{2}_{1}$ | 1.0000 | 0.9778 | 0.8918 | 1.0933
$\sigma^{2}_{2}$ | 0.8000 | 0.7867 | 0.6824 | 0.8885
$\sigma^{2}_{3}$ | 0.5000 | 0.5613 | 0.4988 | 0.6305
$\mu_{1}$ | -9.0000 | -9.0009 | -9.2870 | -8.6325
$\mu_{2}$ | -9.5000 | -9.5044 | -9.6575 | -9.3313
$\mu_{3}$ | -8.5000 | -8.4880 | -8.5960 | -8.4005
$\phi_{1}$ | 0.7000 | 0.6928 | 0.6255 | 0.7494
$\phi_{2}$ | 0.5000 | 0.4927 | 0.4120 | 0.5586
$\phi_{3}$ | 0.3000 | 0.2954 | 0.2402 | 0.3541
Table 5: True values and posterior estimates of parameters for Scenario 2 from the IR-MSV model. Parameter | True Value | Mean | Q(2.5%) | Q(97.5%)
---|---|---|---|---
$\rho_{12}$ | 0.7000 | 0.6985 | 0.6833 | 0.7147
$\rho_{13}$ | 0.7000 | 0.6993 | 0.6859 | 0.7100
$\rho_{23}$ | 0.7000 | 0.7006 | 0.6872 | 0.7139
$\sigma^{2}_{1}$ | 1.0000 | 0.9865 | 0.9045 | 1.0815
$\sigma^{2}_{2}$ | 0.8000 | 0.7811 | 0.7305 | 0.8642
$\sigma^{2}_{3}$ | 0.5000 | 0.4885 | 0.4241 | 0.5620
$\mu_{1}$ | -9.0000 | -8.9865 | -9.1761 | -8.7485
$\mu_{2}$ | -9.5000 | -9.4942 | -9.6221 | -9.3748
$\mu_{3}$ | -8.5000 | -8.4971 | -8.6012 | -8.3849
$\phi_{1}$ | 0.7000 | 0.6998 | 0.6351 | 0.7572
$\phi_{2}$ | 0.5000 | 0.4912 | 0.3955 | 0.5725
$\phi_{3}$ | 0.5000 | 0.4965 | 0.3965 | 0.5785
Table 6: True values and posterior estimates of parameters for Scenario 3 from
the IR-MSV model.

##  4 Data Analysis: Multiple Intraday Log Returns

We implement the proposed models on intraday log returns for three health
stocks traded on the NYSE: Medtronic PLC (MDT), Bristol-Myers Squibb Co (BMY),
and CVS Health Corp (CVS). In particular, we compare the forecast performance
of the IR-MSV model with the IR-SV models fit to each stock.

###  4.1 Data Description

We considered the intra-day prices for MDT, BMY and CVS traded on June 24,
2016. We aggregated the high frequency prices at the one second level by
taking the price with respect to the latest time point as in Buccheri et al.
(2021b). There were 10066, 8797 and 11330 observations (transactions) for MDT,
BMY, and CVS respectively. We synchronized these prices from three stocks
using the refresh time sampling approach (Barndorff-Nielsen et al., 2011) to
obtain a sample size of $T=6239$ synchronized observations. A brief
description of refresh sampling technique is discussed in (Section 2.5.3).
Using these irregularly spaced refreshed prices ($P_{t_{j}}$) we calculate the
log returns as

$r_{t_{j}}=\log(P_{t_{j}})-\log(P_{t_{j-1}}),\quad 2\leq j\leq T$

We have plotted the time series of log returns calculated from the refreshed
prices of MDT, BMY and CVS on $24^{th}$ June, 2016 in Figure 2. They seem to
be fairly correlated.

Figure 2: Log returns of BMY, CVS and MDT calculated from the refreshed prices
on $24^{th}$ June, 2016

###  4.2 Results

We fit IR-MSV model as in equations (17) and (18) to the irregularly spaced
multivariate intraday log returns series of three stocks. We use first 6195
observations for fitting IR-MSV model and kept 44 observations as hold out
samples. We have used exactly the same priors for the parameters as in the
simulation setup. We run 300,000 MCMC iterations and discard the first 50,000
as burn-in to obtain samples from every $500^{th}$ iteration of the last
250,000 iterations. For convergence of MCMC iterations we have used trace and
posterior density plots. In Table 7, we report the posterior sample mean along
with their $95\%$ credible intervals.

Parameter | Mean | SD | 2.50% | 97.50%
---|---|---|---|---
$\rho_{12}$ | 0.276 | 0.013 | 0.249 | 0.301
$\rho_{13}$ | 0.352 | 0.012 | 0.326 | 0.372
$\rho_{23}$ | 0.292 | 0.013 | 0.267 | 0.317
$\sigma^{2}_{1}$ | 0.388 | 0.046 | 0.302 | 0.480
$\sigma^{2}_{2}$ | 0.715 | 0.057 | 0.610 | 0.837
$\sigma^{2}_{3}$ | 0.442 | 0.048 | 0.352 | 0.542
$\mu_{1}$ | -7.286 | 0.054 | -7.390 | -7.176
$\mu_{2}$ | -8.589 | 0.051 | -8.690 | -8.481
$\mu_{3}$ | -8.286 | 0.052 | -8.395 | -8.190
$\phi_{1}$ | 0.403 | 0.062 | 0.286 | 0.529
$\phi_{2}$ | 0.143 | 0.042 | 0.077 | 0.236
$\phi_{3}$ | 0.299 | 0.059 | 0.190 | 0.409
Table 7: Table showing posterior estimates and $95\%$ credible intervals from
IR-MSV model for the intraday log returns of the three stocks.

All parameters lie within the $95\%$ credible intervals limits. The estimated
correlations are moderate indicating fair dependency among the components
which is expected. The persistence parameter for BMY is much higher than that
of CVS and MDT indicating the temporal correlation of volatility of BMY is
much higher than the temporal correlation of volatility of CVS and MDT.

We have also fitted IR-SV models to the univariate refreshed log return series
of each stock of length $T=6195$. We have used the same prior as discussed in
the simulation study of IR-SV models in Section 2.3. We run 550,000 MCMC
iterations and discard the first 50,000 as burn-in to obtain samples from
every 1000th iteration of the last 500,000 iterations. In Table 10, Table 8
and Table 9 we report the posterior sample mean along with their $95\%$
credible intervals. All the parameters lie inside the $95\%$ credible
intervals. Convergence for the parameters are assessed using trace and
posterior density plots. We observe that the persistence parameter $\phi$ for
BMY is higher than that of CVS and MDT which is consistent with the result
obtained in case of IR-MSV model fitting.

Parameter | Mean | SD | 2.50% | 97.50%
---|---|---|---|---
$\mu$ | -7.248 | 0.084 | -7.413 | -7.088
$\phi$ | 0.670 | 0.065 | 0.539 | 0.793
$\sigma$ | 0.486 | 0.042 | 0.403 | 0.568
Table 8: Table showing posterior estimates and $95\%$ credible intervals from IR-SV model for the refreshed intraday log returns of BMY. Parameter | Mean | SD | 2.50% | 97.50%
---|---|---|---|---
$\mu$ | -8.571 | 0.053 | -8.672 | -8.468
$\phi$ | 0.170 | 0.048 | 0.085 | 0.271
$\sigma$ | 0.829 | 0.034 | 0.763 | 0.898
Table 9: Table showing posterior estimates and $95\%$ credible intervals from IR-SV model for the refreshed intraday log returns of CVS. Parameter | Mean | SD | 2.50% | 97.50%
---|---|---|---|---
$\mu$ | -8.267 | 0.059 | -8.381 | -8.148
$\phi$ | 0.381 | 0.077 | 0.231 | 0.528
$\sigma$ | 0.648 | 0.039 | 0.570 | 0.724
Table 10: Table showing posterior estimates and $95\%$ credible intervals from
IR-SV model for the refreshed intraday log returns of MDT.

##  5 Summary and Discussion

In this chapter we extend gap time modeling idea of Nagaraja et al. (2011) to
propose univariate volatility models for irregular financial time series by
modifying the regularly spaced stochastic volatility models. We also extend
this approach to propose multivariate stochastic volatility (MSV) models for
multiple irregularly spaced time series which is a modification of the MSV
model of Chib et al. (2009) that was used with daily data. We used simulation
studies to show the accuracy of estimation of the proposed models. We have
used the proposed models to model intraday logarithmic returns for three
health sector NYSE stocks and compared the univariate and multivariate
forecasting performance on shorter and longer horizons.

We can construct an IR-GARCH model in a similar way. We have illustrated the
accuracy of estimation in IR-GARCH(1,1) models using simulations. This idea
can also be extended to multivariate GARCH models which has been left for
future research.

##  6 Appendix

###  6.1 Explorations with Irregular GARCH models

One of the early approaches to model the return series sampled at irregularly
spaced time intervals set by the trade arrivals is ACD-GARCH model by Ghysels
and Jasiak (1998). This is a random coefficient GARCH, or doubly stochastic
GARCH, where the stochastic durations between transactions determine the
parameter dynamics for the GARCH equation. They proposed a Generalized Method
of Moments (GMM) based two-step procedure for estimating the parameters. ACD-
GARCH is quite cumbersome to estimate and also difficult to generalize to
multiple time series. Meddahi et al. (2006) proposed a GARCH type model for
irregularly spaced data which is an exact discretization of continuous time
stochastic volatility processes observed at irregularly spaced times. Their
model combines the advantages of ACD-GARCH Ghysels and Jasiak (1998) and ACD
(Engle and Russell, 1998). A continuous version of GARCH (COGARCH) is another
way to model irregularly spaced time series data (Maller et al., 2008).
Recently, Buccheri et al. (2021a) propose to model intraday log-prices through
a multivariate local-level model with score-driven covariance matrices and to
treat asynchronicity as a missing value problem.

In this section we do some preliminary investigation of extending the gap time
autoregressive model of Nagaraja et al. (2011) to develop an irregular GARCH
(IR-GARCH) model. We introduce the IR-GARCH model and discuss its properties.
We conduct a small simulation study to assess the accuracy of estimation under
different scenarios. Future research of this topic can be useful.

Let $\\{r_{t_{j}}\\}$ be a sequence of zero mean log-returns of an asset. Here
$t_{j}$ denotes the time of the $j^{th}$ observation and
$g_{j}=t_{j}-t_{j-1}$, $j>1$ is the known $j^{th}$ observed gap time. Suppose
$\\{\epsilon_{t_{j}}\\}$ is a sequence of real valued independent and
identically distributed random variables, having mean 0 and variance 1. We
start with an IR-GARCH(1,1) model as shown below.

#### IR-GARCH(1,1)

Then, $r_{t_{j}}$ follows IR-GARCH$(1,1)$ model if

$r_{t_{j}}=\sigma_{t_{j}}\epsilon_{t_{j}}$ (24)
$\begin{gathered}\sigma_{t_{1}}^{2}=\omega(1-\alpha_{1}^{g_{1}}-\beta_{1}^{g_{1}}),\\\
\sigma_{t_{j}}^{2}=\omega(1-\alpha_{1}^{g_{j}}-\beta_{1}^{g_{j}})+\alpha_{1}^{g_{j}}r_{t_{j-1}}^{2}+\beta_{1}^{g_{j}}\sigma_{t_{j-1}}^{2},\text{
for $j>1$.}\end{gathered}$ (25)

Let $g*=\min_{j}g_{j}$; then the constraints for $\sigma_{t_{j}}^{2}$ to be
positive are $\omega>0,\alpha_{1}>0,\beta_{1}>0$ and
$\alpha_{1}^{g*}+\beta_{1}^{g*}<1$. The unconditional mean of $r_{t_{j}}$ is

$\displaystyle\text{E}(r_{t_{j}})=\text{E}\big{(}\text{E}\big{(}r_{t_{j}}|\mathcal{F}_{t_{j-1}}\big{)}\big{)}=\text{E}\big{(}\text{E}\big{(}\sigma_{t_{j}}\epsilon_{t_{j}}|\mathcal{F}_{t_{j-1}}\big{)}\big{)}=\text{E}\big{(}\sigma_{t_{j}}\text{E}\big{(}\epsilon_{t_{j}}|\mathcal{F}_{t_{j-1}}\big{)}\big{)}=0,\text{$\forall
j$.}$

The unconditional variance of $r_{t_{j}}$ is

$\text{Var}(r_{t_{1}})=\text{E}(r^{2}_{t_{1}})=\omega(1-\alpha_{1}^{g_{1}}-\beta_{1}^{g_{1}}),\text{
since $E(\epsilon_{t_{1}}^{2})=1$}$ (26)

and for $j>1$,

$\displaystyle\text{Var}(r_{t_{j}})=\text{E}(r^{2}_{t_{j}})=\text{E}[\text{E}[r^{2}_{t_{j}}|\mathcal{F}_{t_{j-1}}]]=E[\sigma^{2}_{t_{j}}]=\omega(1-\alpha_{1}^{g_{j}}-\beta_{1}^{g_{j}})+\alpha_{1}^{g_{j}}\text{E}[r_{t_{j-1}}^{2}]+\beta_{1}^{g_{j}}\text{E}[\sigma_{t_{j-1}}^{2}]$

Assuming stationarity of $\\{r_{t_{j}}\\}$ with $\text{E}(r_{t_{j}})=0$,
$\text{Var}(r_{t_{j}})=\text{Var}(r_{t_{j-1}})=\text{E}(r_{t_{j-1}}^{2})$ we
have

$\text{Var}(r_{t_{j}})=\omega,\text{ provided
$\alpha_{1}^{g*}+\beta_{1}^{g*}<1$. }$ (27)

#### IR-ARCH(1) model

The returns $r_{t_{j}}$ are said to follow an IR-ARCH(1) model if

$\begin{gathered}r_{t_{j}}=\sigma_{t_{j}}\epsilon_{t_{j}},\\\
\sigma_{t_{1}}^{2}=\omega(1-\alpha_{1}^{g_{1}}),\\\
\sigma_{t_{j}}^{2}=\omega(1-\alpha_{1}^{g_{j}})+\alpha_{1}^{g_{j}}r_{t_{j-1}}^{2}\text{,
for $j>1$},\end{gathered}$ (28)

where $g_{j}$’s are the observed gap time between two consecutive observations
and $\\{\epsilon_{t_{j}}\\}$ is a sequence of real valued independent and
identically distributed random variables with mean 0 and variance 1.

######  Proposition 6.1.

Suppose $\\{r_{t_{j}}\\}$ follows an IR-ARCH(1) model defined by (28). Let
$\eta_{t_{1}}=r_{t_{1}}^{2}-\sigma_{t_{1}}^{2}-\omega\alpha_{1}^{g_{1}}$ and
$\eta_{t_{j}}=r_{t_{j}}^{2}-\sigma_{t_{j}}^{2}$, for $j>1$. Suppose
$E[\eta_{t_{1}}^{2}]=C$ and for $j>1$,
$E[\eta_{t_{j}}^{2}]=C(1-\alpha_{1}^{2g_{j}})$, where $C>0$ is a constant.
Then $r_{t_{j}}^{2}$ is a stationary process.

###### Proof.

We have

$\displaystyle\eta_{t_{1}}=r_{t_{1}}^{2}-\sigma_{t_{1}}^{2}-\omega\alpha_{1}^{g_{1}}$
$\displaystyle\Rightarrow\sigma_{t_{1}}^{2}=r_{t_{1}}^{2}-\omega\alpha_{1}^{g_{1}}-\eta_{t_{1}}$

Substituting in equation (28) we get,

$r_{t_{1}}^{2}-\omega=\eta_{t_{1}},$

where $E(\eta_{t_{1}})=E(r_{t_{1}}^{2})-E(\sigma_{t_{1}}^{2})=0$, it follows
from equation(26). Similarly, for $j>1$,

$\displaystyle
r_{t_{j}}^{2}-\eta_{t_{j}}=\omega(1-\alpha_{1}^{g_{j}})+\alpha_{1}^{g_{j}}r_{t_{j-1}}^{2}$
$\displaystyle\Rightarrow
r_{t_{j}}^{2}-\omega=\alpha_{1}^{g_{j}}(r_{t_{j-1}}^{2}-\omega)+\eta_{t_{j}}.$

where
$E(\eta_{t_{j}})=E(r_{t_{j}}^{2}-\sigma_{t_{j}}^{2})=E(E(\sigma_{t_{j}}^{2}\epsilon_{t_{j}}^{2}|\mathcal{F}_{t_{j-1}}))-E(\sigma_{t_{j}}^{2})=E(\sigma_{t_{j}}^{2})-E(\sigma_{t_{j}}^{2})=0$,
for $j>1$. Let $x_{t_{j}}=r_{t_{j}}^{2}-\omega$, along with
$E[\eta_{t_{1}}^{2}]=C$ and for $j>1$,
$E[\eta_{t_{j}}^{2}]=C(1-\alpha_{1}^{2g_{j}})$, where $C>0$ is a constant.
Using Propostion 5.1 Nagaraja et al. (2011), it follows that $x_{t_{j}}$ and
hence $r_{t_{j}}^{2}$ is a stationary process.

∎

#### Estimation

We will use conditional maximum likelihood estimation (Tsay, 2005) by assuming
standard normal distribution for errors for estimating the parameters. Let
$\boldsymbol{\theta}=(\omega,\alpha_{1},\beta_{1})^{\prime}$ be the vector of
parameters for instance, assuming $\epsilon_{t_{j}}$ to follow a standard
normal distribution, the conditional Gaussian log-likelihood is given by

$l_{n}(\boldsymbol{\theta})=-\frac{1}{2}\sum_{j=2}^{n}\Bigg{[}\log(2\pi)+\log(\sigma_{t_{j}}^{2})+\dfrac{r_{t_{j}}^{2}}{\sigma_{t_{j}}^{2}}\Bigg{]},$
(29)

where $\sigma_{t_{j}}^{2}$ is updated according to (25) for the IR-GARCH
model. The estimated parameter vector is

$\displaystyle\hat{\boldsymbol{\theta}}=\max_{\boldsymbol{\theta}}l_{n}(\boldsymbol{\theta}).$
(30)

###  6.2 IR-GARCH(1,1)

We describe two scenarios. For each scenario, we generated 100 sets of zero
mean log returns data, each of length $T=5000$. The gap times $g_{j},1\leq
j\leq T$ are generated using a similar procedure as described in Section 2.3.
The following Table 11 and Table 12 give the true values and conditional ML
estimates of parameters for Scenario 1 and Scenario 2.

Parameter | True Value | Estimate | 2.5% | 97.5%
---|---|---|---|---
$\omega$ | 0.0100 | 0.0102 | 0.0088 | 0.0115
$\alpha$ | 0.7000 | 0.7070 | 0.6479 | 0.7662
$\beta$ | 0.2500 | 0.2427 | 0.1511 | 0.3343
Table 11: True values and conditional ML estimates of parameters for Scenario 1 from the IR-GARCH model. Parameter | True Value | Estimate | 2.5% | 97.5%
---|---|---|---|---
$\omega$ | 0.0100 | 0.0106 | 0.0073 | 0.0138
$\alpha$ | 0.9000 | 0.9053 | 0.8682 | 0.9424
$\beta$ | 0.0500 | 0.0475 | 0.0130 | 0.0820
Table 12: True values and conditional ML estimates of parameters for Scenario
2 from the IR-GARCH model.

We have used R package Rsolnp to optimize the conditional maximum likelihood
function to obtain the parameters in IR-GARCH(1,1) model. All the true values
of the parameters are inside the $95\%$ confidence interval.

## References

  * Barndorff-Nielsen et al. (2011) O. E. Barndorff-Nielsen, P. R. Hansen, A. Lunde, and N. Shephard.  Multivariate realised kernels: consistent positive semi-definite estimators of the covariation of equity prices with noise and non-synchronous trading.  _Journal of Econometrics_ , 162(2):149–169, 2011\.
  * Buccheri et al. (2021a) G. Buccheri, G. Bormetti, F. Corsi, and F. Lillo.  A score-driven conditional correlation model for noisy and asynchronous data: An application to high-frequency covariance dynamics.  _Journal of Business & Economic Statistics_, 39(4):920–936, 2021a.
  * Buccheri et al. (2021b) G. Buccheri, F. Corsi, and S. Peluso.  High-frequency lead-lag effects and cross-asset linkages: a multi-asset lagged adjustment model.  _Journal of Business & Economic Statistics_, 39(3):605–621, 2021b.
  * Chib et al. (2009) S. Chib, Y. Omori, and M. Asai.  Multivariate stochastic volatility.  In _Handbook of financial time series_ , pages 365–400. Springer, 2009.
  * de Valpine et al. (2017) P. de Valpine, D. Turek, C. J. Paciorek, C. Anderson-Bergman, D. T. Lang, and R. Bodik.  Programming with models: writing statistical algorithms for general model structures with nimble.  _Journal of Computational and Graphical Statistics_ , 26(2):403–413, 2017.
  * Engle and Russell (1998) R. F. Engle and J. R. Russell.  Autoregressive Conditional Duration: A New Model for Irregularly Spaced Transaction Data.  _Econometrica_ , 66(5):1127–1162, 1998.
  * Ghysels and Jasiak (1998) E. Ghysels and J. Jasiak.  Garch for irregularly spaced financial data: the acd-garch model.  _Studies in Nonlinear Dynamics & Econometrics_, 2(4), 1998\.
  * Lewandowski et al. (2009) D. Lewandowski, D. Kurowicka, and H. Joe.  Generating random correlation matrices based on vines and extended onion method.  _Journal of multivariate analysis_ , 100(9):1989–2001, 2009.
  * Maller et al. (2008) R. A. Maller, G. Müller, and A. Szimayer.  Garch modelling in continuous time for irregularly spaced time series data.  _Bernoulli_ , 14(2):519–542, 2008.
  * Meddahi et al. (2006) N. Meddahi, E. Renault, and B. Werker.  Garch and irregularly spaced data.  _Economics Letters_ , 90(2):200–204, 2006.
  * Nagaraja et al. (2011) C. H. Nagaraja, L. D. Brown, and L. H. Zhao.  An autoregressive approach to house price modeling.  _The Annals of Applied Statistics_ , pages 124–149, 2011.
  * Tsay (2005) R. S. Tsay.  _Analysis of financial time series_.  John wiley & sons, 2005.