# Modeling Multivariate Positive-Valued Time Series Using R-INLA

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 describe fast Bayesian statistical analysis of vector
positive-valued time series, with application to interesting financial data
streams. We discuss a flexible level correlated model (LCM) framework for
building hierarchical models for vector positive-valued time series. The LCM
allows us to combine marginal gamma distributions for the positive-valued
component responses, while accounting for association among the components at
a latent level. We use integrated nested Laplace approximation (INLA) for fast
approximate Bayesian modeling via the R-INLA package, building custom
functions to handle this setup. We use the proposed method to model
interdependencies between realized volatility measures from several stock
indexes.

Keywords: Positive time series, Level correlated models, Multivariate gamma,
Approximate Bayesian inference, INLA, Realized volatility

##  1 Introduction

Multivariate positive-valued time series are ubiquitous across diverse
application domains, including but not limited to, epidemiology, econometrics,
finance, insurance, and signal processing. In financial applications, we
increasingly observe long multivariate positive-valued time series from
multiple stocks. Examples include time series related to high-frequency
trading, such as volume, price, realized volatility measures, durations, etc.
While there has been a considerable amount of research on modeling univariate
positive-valued time series with exponential and gamma marginals (Lawrance and
Lewis (1980, 1981, 1985), Jacobs and Lewis (1977), Raftery (1982), Gaver and
Lewis (1980) and Gouriéroux and Jasiak (2006)), research on the multivariate
setup is rather limited. A few references which deal with modeling
multivariate positive-valued time series include Prékopa and Szántai (1978),
Yue (2001), Cipollini et al. (2006) and Hautsch et al. (2013). While these
static modeling approaches may be useful when the time series are stationary,
many interesting situations in practice involve nonstationary time series and
will benefit from fitting dynamic models.

Gaussian dynamic models have been widely used for several decades for modeling
univariate and multivariate time series (West and Harrison, 2006, Petris,
2010). More recently, dynamic generalized linear models (DGLM’s) have been
tailored to model non-Gaussian time series without the necessity of data
transformations, or other adaptations; see West et al. (1985), Grunwald et al.
(1993), Fahrmeir and Kaufmann (1987), Lindsey and Lambert (1995), Gamerman
(1998), Chiogna et al. (2002), Shephard and Pitt (1997), Benjamin et al.
(2003) and Godolphin and Triantafyllopoulos (2006). Fahrmeir et al. (1994) and
Durbin and Koopman (2012) are excellent books on this topic.

Non-Gaussian dynamic modeling with marginals having positive support (such as
gamma, lognormal, or Weibull) is an attractive setup for analyzing
multivariate positive-valued series exhibiting temporal nonstationarity, as
well as skewness and heavy-tails. Univariate gamma regression models for each
of the components of the response vector cannot account for the dependence
among the components, where the dependence may be due to omitted variables
which simultaneously effect the response vector (Aitchison and Ho, 1989, Chib
and Winkelmann, 2001, Ma et al., 2008). Therefore, an adequate multivariate
model for positive time series is needed.

To our knowledge, Aktekin et al. (2020) is the only paper describing a dynamic
model for multivariate positive-valued time series. It accounts for
correlations among the components through a common random environment that
evolves over time. They develop MCMC methods as well as sequential Monte Carlo
methods for parameter estimation. Tsionas (2004) described Bayesian analysis
using a multivariate gamma likelihood (Mathai and Moschopoulos, 1991), non-
conjugate prior specifications, and Gibbs sampler with data augmentation. This
is not a dynamic model, and moreover, computing the likelihood function can be
cumbersome and time consuming. Incorporating the dependence among the
components of the positive-valued vector through a different mechanism will be
useful. In general, using Markov chain Monte Carlo (MCMC) methods for Bayesian
hierarchical dynamic modeling of vector positive-valued time series under a
multivariate gamma sampling distributional assumption may be computationally
demanding, especially since there does not exist any convenient form of the
likelihood function.

The level correlated model (LCM) framework described in Serhiyenko (2015) for
multivariate time series of counts offers a computationally simpler way of
handling dependency among the components in vector-valued responses. We adapt
this framework for modeling multivariate positive-valued time series, which
assumes marginal gamma component distributions, and accounts for the
association between the components at a latent Gaussian level. This avoids the
need for computing the multivariate gamma likelihood, which can be
computationally demanding.

Another way to reduce the computational burden is to replace the MCMC approach
by the integrated nested Laplace approximation (INLA) approach. INLA enables
approximate Bayesian inference for a large class of latent Gaussian models
(LGMs) with wide-ranging applications Rue et al. (2009), and unlike simulation
intensive Bayesian approaches, it performs approximate inference using a
series of deterministic approximations that take advantage of the LGM
structure to provide fast and accurate inference. Ruiz-Cárdenas et al. (2012)
presented a general framework which enabled users to use R-INLA for a variety
of state-space models. Serhiyenko et al. (2015, 2018) used multivariate
Gaussian random effects to model level correlated effect in multivariate time
series of counts in the LCM framework thus facilitating the estimation of LCMs
using R-INLA package.

In this article, we describe the LCM framework for dynamic modeling of vector
positive-valued time series. At the observation level, the components are
conditionally independent with marginal gamma distributions, while a level
correlation random effect knits the dependence between components. The latent
state vector evolves as a stationary vector autoregression (VAR) or vector
autoregression and moving average (VARMA) process. We use R-INLA for the data
analysis, developing custom code via the rgeneric() function for implementing
the VAR or VARMA latent models. This setup is useful for financial modeling of
the joint dynamics of realized volatility measures among the stock indexes.

The format of the paper follows. Section 2 describes the general framework of
the level correlated model (LCM). Section 3 gives the description of the
approximate Bayesian analysis using R-INLA including the derivation of the
joint distribution of the latent states in subsection 3.1, followed by the
prior specification and model fitting using R-INLA in subsection 3.2. For
simulated data, Section 4 demonstrates that the R-INLA approach enables
accurate model fitting. Section 5 presents an analysis of the multivariate
realized measures of volatility. Finally, Section 6 concludes with discussion
and directions for future research.

##  2 Level Correlated Model Framework

Let $\mathbf{y}_{it}=(Y_{1,it},\dots,Y_{m,it})^{\prime}$ be a $m$-variate
vector of continuous positive responses for $i=1,\dots,n$ and $t=1,\dots,T$.
In the financial example, we would observe $m$ components of positive
responses on $n$ stocks over $T$ days. We describe the LCM setup for positive
responses similar to the description in Serhiyenko (2015). In the following,
we let $i=1\dots n$, $t=1\dots T$, and $j=1\dots m$. Assuming gamma marginals
for the components, the observation equation can be written as

$Y_{j,it}|\theta_{j,it},\tau\sim\text{Gamma}\bigg{(}\tau,\dfrac{\tau}{\theta_{j,it}}\bigg{)},$
(1)

where the precision parameter $\tau>0$ is static, and the random mean
$\theta_{j,it}>0$. Then,
$\text{E}(Y_{j,it}|\theta_{j,it},\tau)=\theta_{j,it}$, and
$\text{Var}(Y_{j,it}|\theta_{j,it},\tau)=\theta_{j,it}^{2}/\tau$. A dynamic
model for the random mean $\theta_{j,it}$ is

$\log(\theta_{j,it})=\beta_{j,i0}+x_{j,t}+\boldsymbol{S}^{\prime}_{j,it}\boldsymbol{\beta}_{j,i}+\alpha_{j,it},$
(2)

where the random effect $\beta_{j,i0}$ is a subject-specific intercept
corresponding to the $j$th component of the response, $x_{j,t}$ is a time-
varying latent state vector which depends on the $j$th component, and the
random effect $\alpha_{j,it}$ is a response type, time, and subject-specific
level correlated error component. The vector $\boldsymbol{S}_{j,it}^{\prime}$
denotes a $p_{j}$-dimensional vector of covariates, and
$\boldsymbol{\beta}_{j,i}$ is a corresponding $p_{j}$-dimensional vector of
coefficients. In general, $\boldsymbol{\beta}_{j,i}$ can also be treated as
time-varying. We have chosen to use the logarithmic link function in (2)
corresponding to the gamma marginals in (1), but other suitable link functions
may also be used. Let the level correlation vector
$\boldsymbol{\alpha}_{it}=(\alpha_{1,it},\dots,\alpha_{m,it})^{\prime}$. The
dependence between the components of $\boldsymbol{y}_{it}$ can be introduced
via $\boldsymbol{\alpha}_{it}\sim N(\boldsymbol{0},\boldsymbol{\Sigma})$,
where $\boldsymbol{\Sigma}$ is a is a variance-covariance matrix for the level
correlated error. More generally, we can allow $\boldsymbol{\Sigma}$ to be a
subject-specific parameter and assume that $\boldsymbol{\alpha}_{it}\sim
N(\boldsymbol{0},\boldsymbol{\Sigma}_{i})$.

We specify the evolution of the latent state vector
$\boldsymbol{x}_{t}=(x_{1,t},\dots,x_{m,t})^{\prime}$ via a vector
autoregression, VAR$(p)$, given by

$\boldsymbol{x}_{t}=\sum_{h=1}^{p}\boldsymbol{\Phi}_{h}\boldsymbol{x}_{t-h}+\boldsymbol{w}_{t}=\boldsymbol{\Phi}(B)\boldsymbol{x}_{t}+\boldsymbol{w}_{t},$
(3)

where $t=1\dots T$, $\boldsymbol{\Phi}_{h}$ are coefficient matrices of order
$m\times m$ and the error vector
$\boldsymbol{w}_{t}\sim\text{N}(\boldsymbol{0},\boldsymbol{W})$,
$\boldsymbol{W}$ being a $m\times m$ state error variance-covariance matrix.
In (3), $\boldsymbol{\Phi}(B)$ is the VAR($p$) operator defined as

$\boldsymbol{\Phi}(B)=\boldsymbol{I}_{m}-\boldsymbol{\Phi}_{1}B-\dots-\boldsymbol{\Phi}_{p}B^{p},$
(4)

where $B$ is the back-shift operator such that
$B^{j}\boldsymbol{x}_{t}=\boldsymbol{x}_{t-j}$, for any $j\geq 1$. The
VAR($p$) process is stationary if the roots of the determinant of
$\boldsymbol{\Phi}(z)$ are outside the unit circle, i.e., the modulus of each
root is greater than 1, (Lütkepohl, 2005). We refer to the model defined by
(1)-(3) as an LCM-VAR, with independent gamma marginals at the observation
level. Let
$\boldsymbol{\Theta}^{*}=(\beta_{j,i0},\boldsymbol{\beta}^{\prime}_{j,i},\text{vec}^{\prime}(\boldsymbol{\Sigma}),\text{vec}^{\prime}(\boldsymbol{\Phi}_{1}),\dots,\text{vec}^{\prime}(\boldsymbol{\Phi}_{p}),\text{vec}^{\prime}(\boldsymbol{W}))^{\prime}$
denote the set of all LCM-VAR parameters, where vec is the vectorization
operator.

If the $\boldsymbol{\Phi}_{h}$ matrices in (3) are all assumed to be diagonal,
the LCM-VAR model simplifies to $m$ separate LCM-AR models, again with
independent gamma marginals. Here, each component of the latent state vector
$\boldsymbol{x}_{t}$ evolves according to an autoregressive (AR) process of
order $p$, and the state equation (3) becomes

$x_{j,t}=\sum_{h=1}^{p}\phi_{j,h}x_{j,(t-h)}+w_{j,t},~{}j=1,\ldots,m,$ (5)

where the error $w_{j,t}\sim\text{N}(0,W_{j})$, for scalar $W_{j}$.

##  3 Approximate Bayesian Analysis using INLA

In the LCM framework, it is possible to use R-INLA templates when the state
equation follows (5), but not when it follows the VAR($p$) model in (3). To
implement a model with latent VAR effects, we can use the rgeneric() function
to code the joint p.d.f. of the state vector $\boldsymbol{x}_{t}$ which is
described in the next section.

###  3.1 Joint Distribution of the Latent States

Suppose $\boldsymbol{x}_{t}=(X_{t,1},\dots,X_{t,m})^{\prime}$ is an
$m$-dimensional latent process in (3). Let $\boldsymbol{\Theta}$ denote the
vector of all the parameters. The latent Gaussian Markov field representation
is that of a multivariate Gaussian vector with a sparse precision matrix as in
(Rue and Held, 2005, Palmí-Perales et al., 2021); i.e.,

$vec(\boldsymbol{\Theta})\sim N(\boldsymbol{0},\boldsymbol{\Sigma}^{*}),$ (6)

where $\boldsymbol{0}$ denotes vector of zeros and $\boldsymbol{\Sigma}^{*}$
is the variance-covariance matrix. Since R-INLA works with the precision
matrix, we derive the form of the joint p.d.f. in terms of
$\boldsymbol{\Sigma}^{*-1}$. We first show the form of the joint p.d.f. for
the VAR(1) latent process.

####  3.1.1 Joint p.d.f. for a VAR(1) Process

When $p=1$, (3) simplifies to a VAR(1) process:

$\boldsymbol{x}_{t}=\boldsymbol{\Phi}\boldsymbol{x}_{t-1}+\boldsymbol{w}_{t}.$
(7)

To derive the joint distribution of
$(\boldsymbol{x}_{1}^{\prime},\dots,\boldsymbol{x}_{T}^{\prime})^{\prime}$, we
see that

$\log\pi(\boldsymbol{x}_{t}|\boldsymbol{x}_{t-1})=-\frac{1}{2}\boldsymbol{w}_{t}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{w}_{t}+C,$

where
$\boldsymbol{w}_{t}=\boldsymbol{x}_{t}-\boldsymbol{\Phi}\boldsymbol{x}_{t-1}$,
$\pi(.)$ denotes the probability density function and C is a constant. We
assume diffuse initial conditions for $\boldsymbol{x}_{1}\sim
N(\boldsymbol{0},\frac{1}{\kappa}\boldsymbol{I}_{m})$, where $\kappa=0.01$.
The joint distribution of
$(\boldsymbol{x}_{1}^{\prime},\dots,\boldsymbol{x}_{T}^{\prime})^{\prime}$ is
given by

$\displaystyle\log\pi(\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots\boldsymbol{x}_{n})=\sum_{t=2}^{n}\log\pi(\boldsymbol{x}_{t}|\boldsymbol{x}_{t-1})+\log\pi(\boldsymbol{x}_{1})$
$\displaystyle=-\dfrac{1}{2}\big{[}(\boldsymbol{x}_{n}-\boldsymbol{\Phi}\boldsymbol{x}_{n-1})^{\prime}\boldsymbol{W}^{-1}(\boldsymbol{x}_{n}-\boldsymbol{\Phi}\boldsymbol{x}_{n-1})+(\boldsymbol{x}_{n-1}-\boldsymbol{\Phi}\boldsymbol{x}_{n-2})^{\prime}\boldsymbol{W}^{-1}(\boldsymbol{x}_{n-1}-\boldsymbol{\Phi}\boldsymbol{x}_{n-2})$
$\displaystyle\hskip
85.35826pt+\dots+(\boldsymbol{x}_{2}-\boldsymbol{\Phi}\boldsymbol{x}_{1})^{\prime}\boldsymbol{W}^{-1}(\boldsymbol{x}_{2}-\boldsymbol{\Phi}\boldsymbol{x}_{1})+\boldsymbol{x}_{1}^{\prime}\kappa\boldsymbol{I}_{m}\boldsymbol{x}_{1}\big{]}+C^{\prime}$
$\displaystyle=-\dfrac{1}{2}\big{[}\boldsymbol{x}_{n}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{x}_{n}-\boldsymbol{x}_{n}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{\Phi}\boldsymbol{x}_{n-1}-\boldsymbol{x}_{n-1}^{\prime}\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{x}_{n}+\boldsymbol{x}_{n-1}^{\prime}(\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}+\boldsymbol{W}^{-1})\boldsymbol{x}_{n-1}$
$\displaystyle\hskip
85.35826pt+\dots+\boldsymbol{x}_{2}^{\prime}(\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{\Phi}\boldsymbol{x}_{2}+\boldsymbol{W}^{-1})\boldsymbol{x}_{2}-\boldsymbol{x}_{2}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{\Phi}^{\prime}\boldsymbol{x}_{1}-\boldsymbol{x}_{1}^{\prime}\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{x}_{2}$
$\displaystyle\hskip
85.35826pt+\boldsymbol{x}_{1}^{\prime}(\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{\Phi}+\kappa\boldsymbol{I}_{m})\boldsymbol{x}_{1}\big{]}+C^{\prime}$

Let
$\boldsymbol{A}=\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}\boldsymbol{\Phi}$,
$\boldsymbol{B}=-\boldsymbol{\Phi}^{\prime}\boldsymbol{W}^{-1}$ and
$\boldsymbol{C}=\boldsymbol{W}^{-1}$ be the $m\times m$ matrices and let
$\boldsymbol{0}_{m}$ be the $m\times m$ matrix of 0’s. Hence the precision
matrix for the joint distribution of
$(\boldsymbol{x}_{1}^{\prime},\dots,\boldsymbol{x}_{T}^{\prime})^{\prime}$ is
a $mT\times mT$ matrix given by

$\boldsymbol{\Sigma}^{*-1}=\begin{bmatrix}\boldsymbol{A}+\kappa\boldsymbol{I}_{m}&\boldsymbol{B}&\cdots&\cdots&\cdots&\cdots&\cdots&\boldsymbol{0}_{m}\\\
\boldsymbol{B}^{\prime}&\boldsymbol{A+C}&\boldsymbol{B}&\ddots&&&&\vdots\\\
\vdots&\boldsymbol{B}^{\prime}&\boldsymbol{A+C}&\boldsymbol{B}&&&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&&\vdots\\\
\vdots&&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&&\ddots&\ddots&\ddots&\ddots&\boldsymbol{0}_{m}\\\
\vdots&&&&\ddots&\ddots&\ddots&\boldsymbol{B}\\\
\boldsymbol{0}_{m}&\cdots&\cdots&\cdots&\cdots&\boldsymbol{0}_{m}&\boldsymbol{B}^{\prime}&\boldsymbol{C}\\\
\end{bmatrix}$ (8)

####  3.1.2 Joint p.d.f. for a VAR($p$) Process

Following Lütkepohl (2005), we can represent a VAR($p$) process of an
$m$-dimensional time series (in (3)) as a VAR(1) process of a $pm$-dimensional
time series, by stacking
$\boldsymbol{x}_{t},\boldsymbol{x}_{t-1},\dots\boldsymbol{x}_{t-p+1}$:

$\displaystyle\begin{pmatrix}\boldsymbol{x}_{t}\\\ \boldsymbol{x}_{t-1}\\\
\boldsymbol{x}_{t-2}\\\ \vdots\\\ \vdots\\\ \boldsymbol{x}_{t-p+1}\\\
\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\boldsymbol{\Phi}_{1}&\boldsymbol{\Phi}_{2}&\cdots&\cdots&\boldsymbol{\Phi}_{p-1}&\boldsymbol{\Phi}_{p}\\\
\boldsymbol{I}_{m}&\boldsymbol{0}_{m}&\ddots&&\boldsymbol{0}_{m}&\boldsymbol{0}_{m}\\\
\boldsymbol{0}_{m}&\boldsymbol{I}_{m}&\ddots&\ddots&\boldsymbol{0}_{m}&\boldsymbol{0}_{m}\\\
\vdots&\vdots&\ddots&\ddots&\ddots&\vdots\\\
\boldsymbol{0}_{m}&\boldsymbol{0}_{m}&\cdots&\cdots&\boldsymbol{I}_{m}&\boldsymbol{0}_{m}\\\
\end{pmatrix}\begin{pmatrix}\boldsymbol{x}_{t-1}\\\ \boldsymbol{x}_{t-2}\\\
\boldsymbol{x}_{t-3}\\\ \vdots\\\ \vdots\\\ \boldsymbol{x}_{t-p}\\\
\end{pmatrix}+\begin{pmatrix}\boldsymbol{w}_{t}\\\ \boldsymbol{0}\\\
\boldsymbol{0}\\\ \vdots\\\ \vdots\\\ \boldsymbol{0}\\\ \end{pmatrix},\text{
i.e.,}$ $\displaystyle\boldsymbol{x}^{*}_{t}$ $\displaystyle=$
$\displaystyle\boldsymbol{\Phi}^{*}\boldsymbol{x}^{*}_{t-1}+\boldsymbol{w}^{*}_{t},$
(9)

where $\boldsymbol{w}^{*}_{t}\sim SN(\boldsymbol{0},\boldsymbol{W}^{*})$,
$\boldsymbol{W}^{*}=\begin{pmatrix}\boldsymbol{W}&\boldsymbol{0}_{m}&\cdots&\boldsymbol{0}_{m}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{0}_{m}&\boldsymbol{0}_{m}&\cdots&\boldsymbol{0}_{m}\\\
\end{pmatrix}$ is a $pm\times pm$ p.s.d. matrix and $SN$ denotes a singular
normal distribution.

Following the setup for the VAR$(1)$ process in Section 3.1.1, we can then
obtain the $pmT\times pmT$ dimensional precision matrix for the joint
distribution of $T$ realizations from a VAR($p$) process as

$\boldsymbol{\Sigma}^{*-1}=\begin{bmatrix}\boldsymbol{A}^{*}+\kappa^{*}\boldsymbol{I}_{m}&\boldsymbol{B}^{*}&\cdots&\cdots&\cdots&\cdots&\cdots&\boldsymbol{0}_{m}\\\
\boldsymbol{B}^{*\prime}&\boldsymbol{A}^{*}+\boldsymbol{C}^{*}&\boldsymbol{B}^{*}&\ddots&&&&\vdots\\\
\vdots&\boldsymbol{B}^{*\prime}&\boldsymbol{A}^{*}+\boldsymbol{C}^{*}&\boldsymbol{B}^{*}&&&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&&\vdots\\\
\vdots&&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&&\ddots&\ddots&\ddots&\ddots&\boldsymbol{0}_{m}\\\
\vdots&&&&\ddots&\ddots&\ddots&\boldsymbol{B}^{*}\\\
\boldsymbol{0}_{m}&\cdots&\cdots&\cdots&\cdots&\boldsymbol{0}_{m}&\boldsymbol{B}^{*\prime}&\boldsymbol{C}^{*}\\\
\end{bmatrix},$ (10)

where,
$\boldsymbol{A}^{*}=\boldsymbol{\Phi}^{*\prime}\boldsymbol{W}^{*+}\boldsymbol{\Phi}^{*}$,
$\boldsymbol{B}^{*}=-\boldsymbol{\Phi}^{*\prime}\boldsymbol{W}^{*+}$ and
$\boldsymbol{C}^{*}=\boldsymbol{W}^{*+}$ are $pm\times pm$ matrices and
$\boldsymbol{W}^{*+}$ denotes the Moore-Penrose inverse. Let
$\boldsymbol{\Psi}=(\text{vec}(\boldsymbol{\Phi}_{1})^{\prime},\dots,\text{vec}(\boldsymbol{\Phi}_{p})^{\prime},\text{vec}(\boldsymbol{W})^{\prime})$
be the vector of parameters and let
$\boldsymbol{x}=(\boldsymbol{x}_{1}^{*^{\prime}},\dots,\boldsymbol{x}_{T}^{*^{\prime}})^{\prime}$
be the observations from a VAR$(p)$ process. The joint p.d.f. is given by

$f(\boldsymbol{x}|\boldsymbol{\Psi})=(2\pi)^{-pmT/2}|\boldsymbol{\Sigma}^{*}|^{-1/2}\exp\bigg{(}{-\frac{1}{2}\boldsymbol{x}^{\prime}\boldsymbol{\Sigma}^{*-1}\boldsymbol{x}}\bigg{)}.$
(11)

####  3.1.3 Joint p.d.f. for a VARMA($p,q$) Process

In the literature, we rarely see applications where the state evolution in a
dynamic model follows a VARMA process. Nevertheless, for completeness, we
derive the precision matrix for the joint p.d.f. of an $m$-dimensional time
series following a stationary VARMA($p,q$) process defined by

$\boldsymbol{\Phi}(B)\boldsymbol{x}_{t}=\boldsymbol{\Theta}(B)\boldsymbol{w}_{t},$
(12)

where $\boldsymbol{\Phi}(B)$ is the VAR operator as defined in (4) and the
moving average (MA) operator $\boldsymbol{\Theta}(B)$ is defined as

$\boldsymbol{\Theta}(B)=\boldsymbol{I}_{m}-\boldsymbol{\Theta}_{1}B-\dots-\boldsymbol{\Theta}_{q}B^{q}.$

In (12), $\boldsymbol{\Phi}_{i}$ and $\boldsymbol{\Theta}_{i}$ are $m\times m$
unknown parameter matrices and the $m$-dimensional white noise process is
represented by $\boldsymbol{w}_{t}=(w_{t,1},\dots,w_{t,m})$ with
$E(\boldsymbol{w}_{t})=\boldsymbol{0}$ and
$Cov(\boldsymbol{w}_{t})=\boldsymbol{\Sigma}_{w}$, which a p.d. variance-
covariance matrix. Let
$\boldsymbol{\Psi}^{*}=(\text{vec}(\boldsymbol{\Phi}_{1})^{\prime},\dots,\text{vec}(\boldsymbol{\Phi}_{p})^{\prime},\text{vec}(\boldsymbol{\Theta}_{1})^{\prime},\dots,\text{vec}(\boldsymbol{\Theta}_{q})^{\prime},\text{vec}(\boldsymbol{\Sigma}_{w})^{\prime})^{\prime}$
be the vector of parameters and assuming Gaussian errors, the joint p.d.f. is
given by

$L(\boldsymbol{\Psi}^{*}|\boldsymbol{x})=(2\pi)^{-mT/2}|\boldsymbol{\Sigma}|^{-1/2}\exp\bigg{(}-\frac{1}{2}\boldsymbol{x}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}\bigg{)},$
(13)

where $T$ is the number of observations,
$\boldsymbol{x}=(\boldsymbol{x}^{\prime}_{1},\dots,\boldsymbol{x}^{\prime}_{T})^{\prime}$
is $mT$-dimensional vector and $\boldsymbol{\Sigma}$ is the $mT\times mT$
variance-covariance matrix of the observations. Using the closed form
expressions in Gallego (2009), we have

$L(\boldsymbol{\Psi}^{*}|\boldsymbol{x})=L(\boldsymbol{\Psi}^{*}|\hat{\boldsymbol{w}}_{0})=(2\pi)^{-mT/2}|\boldsymbol{\Sigma}_{0}|^{-1/2}\exp\bigg{(}-\frac{1}{2}\hat{\boldsymbol{w}}^{\prime}_{0}\boldsymbol{\Sigma}^{-1}_{0}\hat{\boldsymbol{w}}_{0}\bigg{)},$
(14)

where $\boldsymbol{\Sigma}_{0}$ is the variance-covariance matrix of
$\hat{\boldsymbol{w}}_{0}$. Let
$\boldsymbol{u}=\left(\boldsymbol{x}_{-p+1},\ldots,\boldsymbol{x}_{0},\boldsymbol{w}_{-q+1},\ldots,\boldsymbol{w}_{0}\right)$
is a $m(p+q)\times 1$ vector of pre-sample values and covariance of pre-sample
values be given by $\boldsymbol{\Sigma}_{u}$. The block matrix
$\boldsymbol{\Phi}=\left[\boldsymbol{\Phi}_{ij}\right](i=1,\ldots,T;j=1,\ldots,T)$
is defined as $\boldsymbol{\Phi}_{ij}=-\boldsymbol{\Phi}_{i-j}$ for $0\leq
i-j\leq p$ and $\boldsymbol{\Phi}_{ij}=\mathbf{0}$ otherwise, where the
$m\times m$ matrix $\boldsymbol{\Phi}_{ij}$ is the typical block of
$\boldsymbol{\Phi}$ and $\boldsymbol{\Phi}_{0}=-\mathbf{I}_{m}$. Similarly,
the block matrix
$\boldsymbol{\Theta}=\left[\boldsymbol{\Theta}_{ij}\right](i=1,\ldots,T;j=1,\ldots,T)$
is such that $\boldsymbol{\Theta}_{ij}=-\boldsymbol{\Theta}_{i-j}$ for $0\leq
i-j\leq q$ and $\boldsymbol{\Theta}_{ij}=\mathbf{0}$ otherwise, where
$\boldsymbol{\Theta}_{0}=-\mathbf{I}_{m}$. Let $r=\max(p,q)$ and the block
matrix $\mathbf{F}=\left[\mathbf{F}_{ij}\right](i=1,\ldots,r;j=1,\ldots,p+q)$
is such that $\mathbf{F}_{ij}=\boldsymbol{\Phi}_{p-(j-i)}$ for $0\leq
j-i<p,\mathbf{F}_{ij}=-\boldsymbol{\Theta}_{p+q-(j-i)}$ for $p\leq j-i<p+q$,
and $\mathbf{F}_{ij}=\mathbf{0}$ otherwise. The $kT\times kr$ matrix
$\mathbf{G}$ has 1s along the main diagonal and 0s otherwise. Hence the
precision matrix as in Gallego (2009) is given by

$\boldsymbol{\Sigma}_{0}^{-1}=\left(\mathbf{I}_{T}\otimes\boldsymbol{\Sigma}_{w}^{-1}\right)-\left(\mathbf{I}_{T}\otimes\boldsymbol{\Sigma}_{w}^{-1}\right)\boldsymbol{\Theta}^{-1}\mathbf{GF}\left(\mathbf{Z}^{\prime}\mathbf{Z}\right)^{-1}\mathbf{F}^{\prime}\mathbf{G}^{\prime}\Theta^{\prime-1}\left(\mathbf{I}_{T}\otimes\boldsymbol{\Sigma}_{w}^{-1}\right),$
(15)

where

$\mathbf{Z}^{\prime}\mathbf{Z}=\boldsymbol{\Sigma}_{u}^{-1}+\mathbf{F}^{\prime}\mathbf{G}^{\prime}\mathbf{\Theta}^{\prime-1}\left(\mathbf{I}_{T}\otimes\boldsymbol{\Sigma}_{w}^{-1}\right)\boldsymbol{\Theta}^{-1}\mathbf{GF}$

is a square matrix of order $m(p+q)$.

###  3.2 Approximate Posterior Analysis

Let $\boldsymbol{Y}$ and $\boldsymbol{S}$ respectively denote the responses
$\boldsymbol{y}_{it}$ and the static predictors
$\boldsymbol{S}^{\prime}_{j,it}$ for $i=1,\dots,n$, $j=1,\dots,m$ and
$t=1,\dots,T$. Let
$\boldsymbol{b}=(\beta_{ji0},\boldsymbol{\beta}_{j,i}^{\prime})^{\prime}$
denote all the subject-specific intercepts $\beta_{ji0}$ and the fixed effects
$\boldsymbol{\beta}_{j,i}$, for $i=1,\dots,n$ and $j=1,\dots,m$,
$\boldsymbol{\Phi}=(\text{vec}^{\prime}(\boldsymbol{\Phi}_{1}),\dots,\text{vec}^{\prime}(\boldsymbol{\Phi}_{p}))^{\prime}$
denote all the parameters related to the coefficient matrices in the latent
states and let
$\boldsymbol{x}=(\boldsymbol{x}^{\prime}_{1},\dots,\boldsymbol{x}^{\prime}_{T})^{\prime}$
denotes the latent states. Let $\tau$ be the precision parameter of the
observation equation and
$\boldsymbol{\theta}_{it}=(\theta_{1,it},\dots,\theta_{m,it})^{\prime}$, where
$\theta_{j,it}$ evolves as in (2) for $j=1,\dots,m$. The likelihood function
under the model described by equations (1) - (3) is:

$\displaystyle
L(\boldsymbol{b},\boldsymbol{\Phi},\boldsymbol{W},\boldsymbol{\Sigma},\boldsymbol{x},\tau|\boldsymbol{Y},\boldsymbol{S})=\int\prod_{i=1}^{n}\prod_{t=1}^{T}p(\boldsymbol{y}_{it}|\boldsymbol{\theta}_{it},\tau)\times
p(\boldsymbol{x}_{t}|\boldsymbol{x}_{t-1},\dots,\boldsymbol{x}_{t-p},\boldsymbol{W})\times
p(\boldsymbol{\alpha}_{it}|\boldsymbol{\Sigma})\,d\boldsymbol{\alpha}_{it},$
(16)

where

$p(\boldsymbol{y}_{it}|\boldsymbol{\theta}_{it},\tau)=\prod_{j=1}^{m}\bigg{(}\dfrac{\tau^{\tau}}{\Gamma(\tau)\theta_{j,it}^{\tau}}y_{j,it}^{\tau-1}\exp\Big{(}-\frac{\tau}{\theta_{j,it}}y_{j,it}\Big{)}\bigg{)},$
$p(\boldsymbol{\alpha}_{it}|\boldsymbol{\Sigma})=\dfrac{1}{\sqrt{(2\pi)^{m}\det(\boldsymbol{\Sigma})}}\exp\bigg{(}-\dfrac{1}{2}\Big{(}\boldsymbol{\alpha}_{it}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{\alpha}_{it}\Big{)}\bigg{)}$

and
$\prod_{t=1}^{T}p(\boldsymbol{x}_{t}|\boldsymbol{x}_{t-1},\dots,\boldsymbol{x}_{t-p},\boldsymbol{W})$
is the joint p.d.f. of VAR$(p)$ process and hence can be computed as in (11).

The LCM model parameters include the subject-specific intercepts
$\beta_{j,i0}$, the fixed effects $\boldsymbol{\beta}_{j,i}$ – a $p_{j}$
dimensional vector corresponding to $\boldsymbol{S}^{\prime}_{j,it}$, the
$m\times m$ variance-covariance matrix of the random effects
$\boldsymbol{\Sigma}$ and the $m\times m$ variance-covariance matrix of the
state errors $\boldsymbol{W}$. We have mostly used the default prior
specifications of R-INLA which seem to work well. A Wishart distribution is
considered as a prior specification for the $m\times m$ covariance matrix
$\boldsymbol{\Sigma}$. In particular, the precision matrix
$\boldsymbol{\Sigma}^{-1}\sim\text{Wishart}_{m}(r,\boldsymbol{I}_{m})$, $r$ is
the degree of freedom, whose default value is $2m+1$ representing the default
prior specification for $m$-variate Gaussian random effects in R-INLA. For
modeling purposes we have used a diagonal matrix specification for
$\boldsymbol{W}$ = $\text{diag}(\sigma^{2}_{w1},\dots,\sigma^{2}_{wm})$, where
$\sigma^{2}_{wj}$ are marginal variances for $j=1,\dots,m$. Let the
corresponding log-precision be given by
$\theta_{wj}=\log\bigg{(}\dfrac{1}{\sigma^{2}_{wj}}\bigg{)}$. Then for each
$j$, the logarithm of the precision $\theta_{wj}$ is assumed to follow a log-
gamma distribution, whose shape and inverse scale parameters are 1 and 0.00005
respectively. This is also the default prior specification for R-INLA. The
prior for the fixed effects is chosen to be the default Gaussian prior with
mean 0 and precision 0.001.

Given the data and model parameters, the Bayesian formulation requires
specification of the likelihood and prior, from which, using Bayes’ Theorem,
we obtain the posterior density as a normalized product of the likelihood and
prior. Since the likelihood (16) is not available in closed form, its exact
analysis is cumbersome. The high computational requirements of a fully
Bayesian approach has prompted us to use R-INLA for model fitting (Rue et al.,
2009, Lindgren and Rue, 2015). For approximate Bayesian inference, INLA uses
nested Laplace approximations and numerical integration, thereby significantly
reducing the computational time compared to a fully Bayesian inference via
MCMC.

One of the main advantages of INLA is that templates for many useful models in
the literature are available in the R-INLA package. Other cases can be
programmed and implemented using generic functions. We provide more details
about the implementation of VAR($1$) as a latent process using rgeneric in
R-INLA in the Appendix (Appendix). In the following section we show accurate
estimation of parameters using simulated data from trivariate LCM-AR and
trivariate LCM-VAR and then we fit trivariate LCMs to model interdependencies
between the realized volatility measures for several stock indexes.

##  4 Simulated Data

Data from the trivariate LCM-AR model with $Y_{j,it}$ ($i=1,\dots,30$;
$j=1,2,3$ and $t=1,\dots,500$) is simulated according to the following model:

$\displaystyle
Y_{j,it}|\theta_{j,it},\tau\sim\text{Gamma}\bigg{(}\tau,\frac{\tau}{\theta_{j,it}}\bigg{)}$
(17) $\displaystyle\log(\theta_{j,it})=x_{j,t}+\alpha_{j,it}+\beta S_{j,it},$

where $S_{j,it}$ represents a static predictor which is simulated from
$N(0,1)$ and $\beta=0.2$. The state equations are given by (5) with the order
of autoregression $p=1$, $\phi_{1,1}=\phi_{2,1}=\phi_{3,1}=0.8$ and the state
errors $w_{1,t},w_{2,t}$ and $w_{3,t}$ are simulated from
$N(0,W_{1}),N(0,W_{2})\text{ and }N(0,W_{3})$ respectively, where $W_{1}=0.3$,
$W_{2}=0.2$ and $W_{3}=0.5$. The level correlated term
$\boldsymbol{\alpha}_{it}$ is simulated from
$N_{3}(\boldsymbol{0},\boldsymbol{\Sigma})$, where

$\boldsymbol{\Sigma}=\begin{pmatrix}\sigma^{2}_{1}&\rho_{12}\sigma_{1}\sigma_{2}&\rho_{13}\sigma_{1}\sigma_{3}\\\
\rho_{12}\sigma_{1}\sigma_{2}&\sigma^{2}_{2}&\rho_{23}\sigma_{2}\sigma_{3}\\\
\rho_{13}\sigma_{1}\sigma_{3}&\rho_{23}\sigma_{2}\sigma_{3}&\sigma^{2}_{3}\\\
\end{pmatrix},$ (18)

is the variance-covariance matrix with
$\sigma^{2}_{1}=\sigma^{2}_{2}=\sigma^{2}_{3}=0.5$ and $\rho_{12}=0.8$,
$\rho_{13}=0.7$ and $\rho_{23}=0.5$. The responses are simulated from gamma
distribution with $\tau=300$. For the model estimation, we used the R-INLA
package Rue et al. (2009), Martins et al. (2013) with the default prior
specification for the parameters as discussed in subsection 3.2. We assume a
normal prior for $\beta$, trivariate inverse Wishart prior for the variance-
covariance matrix $\boldsymbol{\Sigma}$ and gamma prior for the precisions. We
explicitly assumed $\tau\sim Gamma(0.01,0.01)$ i.e. E[$\tau$] = 1 and
Var[$\tau$]=100.

Simulations were conducted on a personal laptop (11th Gen Intel(R) Core(TM)
i7-1165G7 @ 2.80GHz) with 16Gb of RAM using 64-bit version of Windows 11
operating system. It took about 7 minutes to run LCM-AR on the simulated data.
Most of the estimated parameters in LCM-AR are close to their true values,
except for $\tau$, the precision parameter of the observational equation,
whose estimate is a bit off. Table 1 gives the posterior summaries for the
estimated parameters of LCM-AR.

Table 1: Posterior summaries for LCM-AR simulated data Parameter | True value | Posterior Mean | Posterior Std. Dev.
---|---|---|---
$\tau$ | 300 | 455.422 | 109.925
$\beta$ | 0.200 | 0.206 | 0.003
$1/\sigma^{2}_{1}$ | 2.000 | 1.962 | 0.046
$1/\sigma^{2}_{2}$ | 2.000 | 1.909 | 0.037
$1/\sigma^{2}_{3}$ | 2.000 | 1.875 | 0.049
$\rho_{12}$ | 0.800 | 0.795 | 0.009
$\rho_{13}$ | 0.700 | 0.728 | 0.011
$\rho_{23}$ | 0.500 | 0.515 | 0.011
$\phi_{11}$ | 0.800 | 0.730 | 0.010
$\phi_{22}$ | 0.800 | 0.819 | 0.008
$\phi_{33}$ | 0.800 | 0.807 | 0.007
$1/\sigma^{2}_{w1}$ | 3.333 | 3.312 | 0.044
$1/\sigma^{2}_{w2}$ | 5.000 | 5.438 | 0.059
$1/\sigma^{2}_{w3}$ | 2.000 | 2.149 | 0.021

Data from the trivariate LCM-VAR model with $Y_{j,it}$ ($i=1,\dots,30$;
$j=1,2,3$ and $t=1,\dots,500$) is simulated according to observation equation
as in (17) and the state equation is described by

$\boldsymbol{x}_{t}=\boldsymbol{\Phi}\boldsymbol{x}_{t-1}+\boldsymbol{w}_{t},$

where

$\boldsymbol{\Phi}=\begin{pmatrix}0.5&0&0.3\\\ 0.6&0.1&0.5\\\ 0.1&0&0.8\\\
\end{pmatrix}$

and $\boldsymbol{w}_{t}$ is simulated from
$N_{3}(\boldsymbol{0},\boldsymbol{W})$, where
$\boldsymbol{W}=\text{diag}(0.5,0.25,0.5)$, i.e., $W_{1}=0.5$, $W_{2}=0.25$
and $W_{3}=0.5$. The static predictor is simulated from $N(0,1)$ and the
coefficient $\beta=0.2$. The level correlated term $\boldsymbol{\alpha}_{it}$
is simulated from $N_{3}(\boldsymbol{0},\boldsymbol{\Sigma})$, where
$\boldsymbol{\Sigma}$ is defined as in (18) with $\sigma^{2}_{1}=1/3$,
$\sigma^{2}_{2}=1/2$,$\sigma^{2}_{3}=1/3$ and $\rho_{12}=0.6$, $\rho_{13}=0.8$
and $\rho_{23}=0.8$. The responses are simulated from gamma distribution with
$\tau=100$.

It took about 38 minutes to run LCM-VAR on the simulated data. Table 2 gives
the posterior summaries for the estimated parameters of LCM-VAR. Similar to
the LCM-AR simulation setup, most of the estimated parameters are close to
their respective true values. The posterior standard deviation of the
precisions of the state errors in LCM-VAR seems to be more than that in LCM-
AR.

Table 2: Posterior summaries for LCM-VAR simulated data Parameter | True value | Posterior Mean | Posterior Std. Dev.
---|---|---|---
$\tau$ | 100 | 83.321 | 10.360
$\beta$ | 0.200 | 0.201 | 0.002
$1/\sigma_{1}^{2}$ | 3.000 | 3.109 | 0.038
$1/\sigma_{2}^{2}$ | 2.000 | 2.025 | 0.023
$1/\sigma_{3}^{2}$ | 3.000 | 3.115 | 0.038
$\rho_{12}$ | 0.600 | 0.644 | 0.007
$\rho_{13}$ | 0.800 | 0.847 | 0.006
$\rho_{23}$ | 0.800 | 0.843 | 0.005
$\phi_{11}$ | 0.500 | 0.481 | 0.049
$\phi_{21}$ | 0.600 | 0.569 | 0.030
$\phi_{31}$ | 0.100 | 0.013 | 0.036
$\phi_{12}$ | 0.000 | -0.006 | 0.042
$\phi_{22}$ | 0.100 | 0.090 | 0.031
$\phi_{32}$ | 0.000 | -0.041 | 0.039
$\phi_{13}$ | 0.300 | 0.289 | 0.049
$\phi_{23}$ | 0.500 | 0.517 | 0.025
$\phi_{33}$ | 0.800 | 0.855 | 0.028
$1/\sigma^{2}_{w1}$ | 2.000 | 1.783 | 0.095
$1/\sigma^{2}_{w2}$ | 4.000 | 3.557 | 0.207
$1/\sigma^{2}_{w3}$ | 2.000 | 2.314 | 0.124

##  5 Modeling Multivariate Realized Volatility Measures

Volatility measurement is central to asset pricing, asset allocation,
portfolio optimization and risk management. A first attempt to measure
volatility using intraday data was made by Parkinson (1980) for the estimation
of the daily range. Since then the literature has expanded significantly from
the realized volatility of Andersen and Bollerslev (1998) to bipower variation
(Barndorff-Nielsen and Shephard, 2004), realized kernels of Barndorff-Nielsen
et al. (2009) and realized covariance matrices (Aït-Sahalia et al., 2010,
Peluso et al., 2014). As these intraday volatility measures evolved, there has
been a natural complementary effort to build adequate models to describe their
joint dynamics (Engle and Gallo, 2006, Donelli et al., 2021, Cipollini et al.,
2013).

###  5.1 Data description

The focus of this study is to model three robust intraday volatility measures:
median realized variance ($MedRV$), realized kernel variance based on Parzen
kernel ($RK$) and bipower variation ($BPV$) using the different level
correlated model specifications for the $30$ stock indexes data (as shown in
Table 3) available on Oxford-Man Institute’s Quantitative Finance Realized
Library (Heber et al., 2009). These volatility measures are based on
$5$-minute sampling frequency. The sample period spans from January $2$,
$2015$, to December $29$, $2017$. The number of observations for each stock
index is 513. Thus we have $m=3$, $n=30$ and $T=513$. In Table (3), we list
all the stock indexes we have used in our analysis along with their associated
regions.

Table 3: List of stock indexes along with their associated regions Symbol | Name | Region
---|---|---
AEX | AEX index | Netherlands
AORD | All Ordinaries | Australia
BFX | Bell 20 Index | Belgium
BSESN | S&P BSE Sensex | India
BVSP | BVSP BOVESPA Index | United States
DJI | Dow Jones Industrial Average | United States
FCHI | CAC 40 | France
FTMIB | FTSE MIB | Italy
FTSE | FTSE 100 | United Kingdom
GDAXI | DAX | Germany
GSPTSE | S&P/TSX Composite index | Canada
HSI | HANG SENG Index | Hong Kong
IBEX | IBEX 35 Index | United States
IXIC | Nasdaq 100 | United States
KS11 | Korea Composite Stock Price Index | South Korea
KSE | Karachi SE 100 Index | Pakistan
MXX | IPC Mexico | Mexico
N225 | Nikkei 225 | Japan
NSEI | NIFTY 50 | India
OMXC20 | OMX Copenhagen 20 Index | Denmark
OMXHPI | OMX Helsinki All Share Index | Finland
OMXSPI | OMX Stockholm All Share Index | Sweden
OSEAX | Oslo Exchange All-share Index | Norway
RUT | Russel 2000 | United States
SMSI | Madrid General Index | Spain
SPX | S&P 500 Index | United States
SSEC | Shanghai Composite Index | China
SSMI | Swiss Stock Market Index | Switzerland
STI | Straits Times Index | Singapore
STOXX50E | EURO STOXX 50 | Europe

The mean, standard deviation (St.Dev.), skewness, kurtosis and autocorrelation
of lag 1 of $\sqrt{MedRV}$, $\sqrt{RK}$ and $\sqrt{BPV}$ are given in Table
(4). The values of skewness indicates that all of the realized measures are
positively skewed and high values of kurtosis indicate the presence of heavy
tails. All of the realized measures are persistent according to the high
values of their first order autocorrelation.

Table 4: Descriptive Statistics of $\sqrt{MedRV}$, $\sqrt{RK}$ and $\sqrt{BPV}$ | Mean | St.Dev. | Skewness | Kurtosis | AutoCorr(lag1)
---|---|---|---|---|---
$\sqrt{MedRV}$ | 0.004 | 0.003 | 2.314 | 13.094 | 0.750
$\sqrt{RK}$ | 0.007 | 0.005 | 3.487 | 30.008 | 0.583
$\sqrt{BPV}$ | 0.007 | 0.004 | 3.699 | 32.181 | 0.706

We compute the kendall’s correlation matrix for the components of the response
vector in Table 5. The off-diagonal elements indicate high degree of positive
correlation between the realized volatility measures.

Table 5: Kendall’s correlation matrix of $\sqrt{MedRV}$, $\sqrt{RK}$ and $\sqrt{BPV}$ | $\sqrt{MedRV}$ | $\sqrt{RK}$ | $\sqrt{BPV}$
---|---|---|---
$\sqrt{MedRV}$ | 1.00 | 0.62 | 0.77
$\sqrt{RK}$ | 0.62 | 1.00 | 0.72
$\sqrt{BPV}$ | 0.77 | 0.72 | 1.00

###  5.2 Realized Volatility Measures

Let $\\{P_{t_{j}}\\}$ be a sequence of observed intraday prices of an asset.
Consider a time interval $[a,b]$ and the time index $t_{j}$ for
$t_{j}\in[a,b]$. The time is partitioned into $M$ equidistant points in time
$j=1,\dots,M$. At each time point $t_{j}$, the price of an asset is observed.
Then $\\{P_{t_{j}}\\}_{j=1}^{M}$ is the observed sequence of prices of an
asset sampled at frequency $\tau=\frac{b-a}{M-1}$. Hence the intraday returns
over the sub-interval $[t_{j-1},t_{j})$ is defined as
$r_{j}=\log(P_{t_{j}})-\log(P_{t_{j-1}})$. To obtain realized variances at
5-minutes sampling frequency, $\tau=300$ seconds and similarly for 1-minute
sampling frequency, $\tau=60$ seconds.

Let us denote $r_{j,t}$ as the $j^{th}$ intraday return for $t^{th}$ day, for
$t=1,\dots,T$. We define some popular measures of intraday volatility
considered in the literature. Median realized variance ($MedRV$) is an
estimator introduced by Andersen et al. (2012), which is a scaled square of
the median of three consecutive intraday absolute returns, which is calculated
as

$\text{MedRV}_{t}=\dfrac{\pi}{6-4\sqrt{3}+\pi}\dfrac{M}{M-2}\sum_{j=2}^{M-1}\text{median}\big{(}|r_{j-1,t}|,|r_{j,t}|,|r_{j+1,t}|\big{)}^{2}.$

Barndorff-Nielsen et al. (2009) proposed the realized kernels estimator, which
uses kernel methods to combat market microstructure noise. The realized kernel
estimator is

$\displaystyle\text{RK}_{t}=\sum_{h=-H}^{H}k\Big{(}\frac{h}{H+1}\Big{)}\gamma_{h},\quad\text{where}$
$\displaystyle\gamma_{h}=\sum_{i=|h|+1}^{n}r_{i,t}r_{i-|h|,t},$

where $k(.)$ is the Parzen kernel function defined as

$k(x)=\begin{cases}1-6x^{2}+6x^{3}&0\leq x\leq 1/2\\\ 2(1-x)^{3}&1/2\leq x\leq
1\\\ 0&x>1.\end{cases}$

A simple and comprehensible measure, which isolates the variation of the pure
price process is the Bipower variation (BPV). It is calculated as the sum of
products of absolute values of two consecutive returns,

$\text{BPV}_{t}=\sum_{j=2}^{M}|r_{j,t}||r_{j-1,t}|.$

For each day $t$, we consider the $m=3$-dimensional time series composed of
standard deviation versions of median realized variance, realized kernel and
bipower variation,
$\boldsymbol{y}_{t}=(\sqrt{RV}_{t},\sqrt{RK}_{t},\sqrt{BPV}_{t})^{\prime}$,
where $t=1,\dots,T.$ We have plotted the time series of $\sqrt{MedRV}$,
$\sqrt{RK}$ and $\sqrt{BPV}$ for two stock indexes S&P 500 and NASDAQ-100 for
illustration.

Figure 1: Plots of $\sqrt{MedRV}$, $\sqrt{RK}$ and $\sqrt{BPV}$ for SP500 and
NASDAQ-100 from January 2, 2015, to December 29, 2017

###  5.3 Model Specification

Let $\boldsymbol{y}_{it}=(y_{1,it},\dots,y_{m,it})^{\prime}$ be a $m$-variate
vector of positive-valued volatility measures over time, for $i=1,\ldots,n$
and $t=1,\ldots,T$. That is, we observe $m$ components of different volatility
measures on $n$ stock market indexes over $T$ regularly spaced time points.
High values of the lag 1 autocorrelation has prompted us to use
$\log(Y_{j,i,t-1})$ as the predictor indicating the daily level component
driven by short term traders as in Corsi (2009). We define the daily jump
component as $J_{t}=\max(RV_{t}-BPV_{t},0)$ and the daily continuous component
as $C_{t}=RV_{t}-J_{t}$, where $RV_{t}=\sum_{j=1}^{M}r_{j,t}^{2}$, is the
realized variance at time $t$. The daily jump and continuous components
account for the jumps in the price process, and are useful predictors in
modeling realized volatility (Andersen et al., 2007). With these predictors,
(1) – (2) simplify as

$\displaystyle
Y_{j,it}|\theta_{j,it},\tau\sim\text{Gamma}\bigg{(}\tau,\frac{\tau}{\theta_{j,it}}\bigg{)},$
(19) $\displaystyle\log(\theta_{j,it})=\alpha_{j,i}+\beta_{j,t,0}+$
$\displaystyle\beta_{1,j}\log(Y_{j,i,t-1})+\beta_{2,j}\log(1+J_{j,i,t-1})+\beta_{3,j}\log(1+C_{j,i,t-1})+\xi_{j,it},$

where the level $\alpha_{j,i}$ is assumed to be stock-specific, and the level
correlated effect
$\boldsymbol{\xi}=(\xi_{1,it},\xi_{2,it},\xi_{3,it})^{\prime}\sim
N(\boldsymbol{0},\boldsymbol{\Sigma})$, with

$\boldsymbol{\Sigma}=\begin{pmatrix}\sigma^{2}_{1}&\rho_{12}\sigma_{1}\sigma_{2}&\rho_{13}\sigma_{1}\sigma_{3}\\\
\rho_{12}\sigma_{1}\sigma_{2}&\sigma^{2}_{2}&\rho_{23}\sigma_{2}\sigma_{3}\\\
\rho_{13}\sigma_{1}\sigma_{3}&\rho_{23}\sigma_{2}\sigma_{3}&\sigma^{2}_{3}\\\
\end{pmatrix}.$

We refer to this model as LCM1, and consider two different specifications for
the time-varying intercept $\beta_{j,t,0}$, similar to (3) and (5). First, an
assumption of VAR evolution leads to

$\boldsymbol{\beta}_{t,0}=\boldsymbol{\Phi}\boldsymbol{\beta}_{t-1,0}+\boldsymbol{w}_{t},\text{
where }\boldsymbol{w}_{t}\sim N(\boldsymbol{0},\boldsymbol{W}),$ (20)

while an AR specification for each component $j=1,\ldots,m$ leads to

$\beta_{j,t,0}=\phi_{j}\beta_{j,t-1,0}+w_{j,t},\text{ where }w_{j,t}\sim
N(0,W_{j}).$ (21)

We refer to (19) and (20) as LCM1-VAR, and to (19) and (21) as LCM1-AR.

We also investigate a setup where the level $\alpha_{j,i}$ is component-
specific but not stock-specific, and define the corresponding model LCM2 as

$\displaystyle Y_{j,it}|\theta_{j,it}\sim$
$\displaystyle\quad\text{Gamma}\bigg{(}\tau,\frac{\tau}{\theta_{j,it}}\bigg{)},$
(22) $\displaystyle\log(\theta_{j,it})=\alpha_{j}+\beta_{j,t,0}+$
$\displaystyle\beta_{1,j}\log(Y_{j,i,t-1})+\beta_{2,j}\log(1+J_{j,i,t-1})+\beta_{3,j}\log(1+C_{j,i,t-1})+\xi_{j,it},$

In this case, we refer to equations (22) and (20) as LCM2-VAR, and to (22) and
(21) as LCM2-AR. The difference between LCM1 and LCM2 lies in the
specification of the level $\alpha$.

###  5.4 Results

Using R-INLA, we fit the four models discussed in Section (5.3): LCM1-VAR,
LCM1-AR, LCM2-VAR, and LCM2-AR, using the same prior specification as
mentioned in Section (3.2). We have assumed $\tau\sim\text{Gamma}(1,0.1)$ such
that $\text{E}[\tau]=10$ and $\text{Var}[\tau]=100$.

In Table 6, we show the mean of the posterior means of the levels from Model
LCM1. We also show the posterior mean of the levels from Model LCM2.

Table 6: Posterior summaries of component-specific levels | LCM1-VAR | LCM1-AR | LCM2-VAR | LCM2-AR
---|---|---|---|---
Parameter | Mean of Posterior Mean | Mean of Posterior Mean | Posterior Mean | Posterior Mean
$\alpha_{1}$ | -4.401 | -4.898 | -3.312 | -3.869
$\alpha_{2}$ | -4.443 | -4.633 | -3.657 | -3.923
$\alpha_{3}$ | -3.864 | -4.450 | -3.062 | -3.583

In Table 7, we show posterior summaries of the covariate coefficients and
hyperparameters for Models LCM1 and LCM2. We observe that the $\beta$’s
corresponding to the daily effects are higher under the AR specification than
those obtained by assuming a VAR(1) state evolution, whereas the estimated
state precisions are higher for the VAR(1) state evolution than the AR(1)
state evolution. The estimated level correlations among the three volatility
measures, i.e., $\rho_{12},\rho_{13}$ and $\rho_{23}$, across all the models
are strong and significant, which empirically confirms dependency between
them. The estimated level correlations from the four models are close to the
Kendall’s tau as discussed in Table (5). Most of the estimated parameters are
significant except for the jump components.

Table 7: Posterior summaries of covariate coefficients and hyperparameters | LCM1-VAR | LCM1-AR | LCM2-VAR | LCM2-AR
---|---|---|---|---
Parameter | Posterior Mean | Posterior Mean | Posterior Mean | Posterior Mean
Covariate coefficients |  |  |  |
$\beta_{1,1}$ | 0.206 | 0.120 | 0.405 | 0.306
$\beta_{1,2}$ | 0.122 | 0.088 | 0.281 | 0.229
$\beta_{1,3}$ | 0.244 | 0.133 | 0.404 | 0.303
$\beta_{2,1}$ | 1.661 | 42.829 | 5.444 | 36.690
$\beta_{2,2}$ | 14.860 | 31.900 | -23.530 | 8.912
$\beta_{2,3}$ | 21.604 | 57.976 | 6.873 | 37.788
$\beta_{3,1}$ | 263.115 | 307.014 | 220.804 | 295.162
$\beta_{3,2}$ | 227.402 | 273.007 | 249.234 | 291.896
$\beta_{3,3}$ | 183.556 | 245.000 | 186.188 | 268.420
Hyperparameters |  |  |  |
$1/\sigma_{1}^{2}$ | 13.034 | 11.177 | 11.472 | 9.030
$1/\sigma_{2}^{2}$ | 9.128 | 7.431 | 7.672 | 6.649
$1/\sigma_{3}^{2}$ | 17.758 | 13.290 | 15.152 | 11.019
$\rho_{12}$ | 0.621 | 0.677 | 0.605 | 0.676
$\rho_{13}$ | 0.811 | 0.854 | 0.802 | 0.859
$\rho_{23}$ | 0.807 | 0.837 | 0.798 | 0.841
$\phi_{11}$ | 2.419 | 0.996 | 0.134 | 0.997
$\phi_{21}$ | 2.564 |  | 0.309 |
$\phi_{31}$ | 2.464 |  | 0.444 |
$\phi_{12}$ | 0.592 |  | 2.482 |
$\phi_{22}$ | 0.623 | 0.996 | 2.496 | 0.998
$\phi_{32}$ | 0.526 |  | 2.362 |
$\phi_{13}$ | -2.318 |  | -2.205 |
$\phi_{23}$ | -2.382 |  | -2.246 |
$\phi_{33}$ | -2.359 | 0.997 | -2.410 | 0.998
$1/W_{1}$ | 428.255 | 20.633 | 768.665 | 32.634
$1/W_{2}$ | 173.920 | 21.040 | 458.162 | 29.387
$1/W_{3}$ | 667.143 | 25.496 | 868.763 | 41.081

###  5.5 Model Comparisons

In volatility modeling, we are mostly interested in out-of-sample forecasting
accuracy. We compare the four models discussed in section 5.4 using popular
measures such as mean absolute percentage error (MAPE) and mean absolute error
(MAE). Let $y_{n_{t}}(l)$ be the $l^{th}$ step-ahead forecast from origin
$n_{t}$ for $y_{t+l}$, $l=1,\ldots,n_{h}$, and let
$e_{n_{t}}(l)=y_{n_{t}+l}-y_{n_{t}}(l)$ be the $l^{th}$ step-ahead forecast
error. MAPE is defined as

$\text{MAPE}=\frac{1}{n_{h}}\sum_{l=1}^{n_{h}}|e_{n_{t}}(l)/y_{n_{t}+l}|.$

MAE is defined as

$\text{MAE}=\frac{1}{n_{h}}\sum_{l=1}^{n_{h}}|e_{n_{t}}(l)|.$

Small values of MAPE and MAE indicate better model. To evaluate forecast
accuracy, for each stock index and each realized volatility measure we
calculated both MAPE and MAE based on 1-step, 5-step and 10-step ahead
forecasts i.e. for $n_{h}=1,5$ and $10$. In Table (8), we compute MAPE and MAE
as the average across stock indexes and across the realized volatility
measures.

Table 8: MAPE and MAE averaged across the 30 stock indexes and also across 3 realized volatility measures Horizon ($n_{h}$) | LCM1-AR | LCM1-VAR | LCM2-AR | LCM2-VAR
---|---|---|---|---
MAPE |  |  |  |
1 | 0.302 | 0.320 | 0.310 | $\boldsymbol{0.282}$
5 | 0.307 | 0.322 | 0.318 | $\boldsymbol{0.252}$
10 | 0.317 | 0.367 | 0.324 | $\boldsymbol{0.280}$
MAE |  |  |  |
1 | 0.0012 | 0.0013 | 0.0013 | $\boldsymbol{0.0011}$
5 | 0.0012 | 0.0014 | 0.0012 | $\boldsymbol{0.0010}$
10 | 0.0013 | 0.0016 | 0.0012 | $\boldsymbol{0.0011}$

Thus, based on MAPE and MAE, LCM2-VAR outperforms all other models. For in-
sample comparison we have used Deviance Information Criterion (DIC)
(Spiegelhalter et al., 2002) as the model selection criterion. In that case,
LCM1-AR was outperforming the rest due to smaller values of DIC, but this
model does not perform well in terms of forecasting accuracy as in Table 8.

##  6 Summary and discussion

This article describes the use of the rgeneric() function in R-INLA for
fitting an LCM with VAR$(p)$ latent effects to vector positive-valued time
series. We demonstrate the method using simulated data and realized financial
volatilities. The proposed LCM-VAR framework is not limited to financial
applications, or to temporal models with exogenous predictors. We can easily
adapt the model to incorporate spatial dependence between different locations,
with applications in modeling climate variables.

Another useful direction is a hierarchical multivariate dynamic model (HMDM)
setup where the observed data vector
$\boldsymbol{y}_{it}=(Y_{1,it},\dots,Y_{m,it})^{\prime}$ follows a
multivariate Gamma distribution. Several multivariate extensions of the
univariate gamma distributions exist in the literature (see Cherian (1941),
Ramabhadran (1951), Mathai and Moschopoulos (1991), and Kotz et al. (2004) and
references therein). Here we define multivariate gamma distribution as defined
in (Tsionas, 2004) which is in the same spirit as in (Mathai and Moschopoulos,
1991). Following (Mathai and Moschopoulos, 1991) and (Tsionas, 2004) to define
the probability distribution function (p.d.f.) of an $m$-variate gamma
distribution, we assume that

$Y_{j,it}=V_{j,it}+Z,$

where $V_{j,it}$ and $Z$ are independent random variables distributed as
Gamma$(\lambda_{j,it},\theta_{j,it})$, $i=1,\dots,n$, $j=1,\dots,m$,
$t=1,\dots,T$ and Gamma$(\alpha,\gamma)$ respectively. Then for $i=1,\dots,n$,
$j=1,\dots,m$ and $t=1,\dots,T$,

$p(y_{j,it}|z,\lambda_{j,it},\theta_{j,it})=\dfrac{\theta_{j,it}^{\lambda_{j,it}}}{\Gamma(\lambda_{j,it})}(y_{j,it}-z)^{\lambda_{j,it}-1}\exp[-\theta_{j,it}(y_{j,it}-z)],\hskip
5.69046pty_{j,it}\geq z,$

and the p.d.f. of $z$ is given by

$p(z|\alpha,\gamma)=\dfrac{\gamma^{\alpha}}{\Gamma(\alpha)}z^{\alpha-1}\exp(-\gamma
z),\hskip 5.69046ptz>0,\alpha>0,\gamma>0.$

The sampling distribution is given by

$p(\boldsymbol{y}_{it}|\alpha,\gamma,\boldsymbol{\lambda}_{it},\boldsymbol{\theta}_{it})=\dfrac{\gamma^{\alpha}}{\Gamma(\alpha)}\prod_{j=1}^{m}\dfrac{\theta_{j,it}^{\lambda_{j,it}}}{\Gamma(\lambda_{j,it})}\int_{0}^{B}(y_{j,it}-z)^{\lambda_{j,it}-1}\exp(-\theta_{j,it}(y_{j,it}-z))z^{\alpha-1}\exp(-\gamma
z)\,dz,$

where $B=\min\\{{y_{j,it}:i=1,\dots,n;j=1,\dots,m;t=1\dots,T\\}}$,
$\boldsymbol{\lambda}_{it}=(\lambda_{1,it},\dots,\lambda_{m,it})^{\prime}$,
$\boldsymbol{\theta}_{it}=(\theta_{1,it},\dots,\theta_{m,it})^{\prime}$ and
$\boldsymbol{\Theta}_{it}=(\boldsymbol{\lambda}_{it},\boldsymbol{\theta}_{it})^{\prime}$
and $\alpha>0$,$\gamma>0$. The second-level observational equation for
$i=1,\dots,n$, $j=1,\dots,m$ and $t=1,\dots,T$ is given by

$\displaystyle\log(\boldsymbol{\Theta}_{j,it})=\boldsymbol{D}^{\prime}_{j,it}\boldsymbol{\delta}_{j,it}+\boldsymbol{S}^{\prime}_{j,it}\boldsymbol{\eta}_{j,it},$

where
$\boldsymbol{D}_{j,it}=\left(D_{j,it,1},\dots,D_{j,it,a_{j}}\right)^{\prime}$
is an $a_{j}$-dimensional vector of exogenous predictors with stock-time
varying (dynamic) coefficients $\boldsymbol{\delta}_{j,it}=$
$\left(\delta_{j,it,1},\cdots,\delta_{j,it,a_{j}}\right)^{\prime}$ and
$\boldsymbol{S}_{j,it}=\left(S_{j,it,1},\cdots,S_{j,it,b_{j}}\right)^{\prime}$
is a $b_{j}$-dimensional vector of exogenous predictors with static
coefficients
$\boldsymbol{\eta}_{j}=\left(\eta_{j,1},\cdots,\eta_{j,it,b_{j}}\right)^{\prime}$.
Let $\boldsymbol{\xi}_{it}$ be a $p_{d}$-dimensional vector constructed by
stacking the $a_{j}$ coefficients $\boldsymbol{\delta}_{j,it}$ for
$j=1,\cdots,m$. The structural equation of the HMDM relates the stock-time
varying parameter $\boldsymbol{\xi}_{it}$ to an aggregate (pooled) state
parameter $\boldsymbol{x}_{t}$ :

$\boldsymbol{\xi}_{it}=\boldsymbol{x}_{t}+\boldsymbol{v}_{it}$

where the errors $\boldsymbol{v}_{it}$ are assumed to be i.i.d.
$N_{p_{d}}\left(\mathbf{0},\boldsymbol{V}_{i}\right)$. The state (or system)
equation of the HMDM is:

$\boldsymbol{x}_{t}=\boldsymbol{G}\boldsymbol{x}_{t-1}+\boldsymbol{w}_{t}$

where $\boldsymbol{G}$ is a $p_{d}\times p_{d}$ state transition matrix and
the state errors $\boldsymbol{w}_{t}$ are assumed to be i.i.d.
$N_{p_{d}}(\mathbf{0},\boldsymbol{W})$. This model could be explored in the
context of analyzing the joint dynamics of stock indexes. Although we have not
discussed about incorporating dependence between the components of
$\boldsymbol{Y}$ through copulas, this could be an interesting topic for
further exploration.

## Appendix

###  rgeneric() construction for VAR(1)

In this section we briefly explain how to build an rmodel to fit latent VAR(1)
using rgeneric() in R-INLA. We have to define rmodel that we label
‘inla.rgeneric.VAR1.model’. Through the inla.rgeneric.define() function, we
create an inla.rgeneric object that we call model.VAR.

model.VAR <\- inla.rgeneric.define(inla.rgeneric.VAR1.model, …).

Then, we can embed this “inla.rgeneric” object in the usual INLA syntax to fit
latent VAR(1).

formula.inla <\- $\dots+\text{f(idx, model = model.VAR, \ldots)}+$ …,

where idx refers to the indices in R-INLA setup.

Let $\boldsymbol{x}_{t}$ be a $m$-dimensional time series following VAR($1$)
process as in equation (7). Assuming
$\boldsymbol{W}=\text{diag}(\sigma^{2}_{w1},\dots,\sigma^{2}_{wm})$, be the
$m\times m$ covariance matrix of the latent states. The set of parameters in
latent VAR($1$) consists of $\\{\phi_{ij},i=1\dots,m;j=1,\dots,m\\}$ and the
marginal variances $\sigma^{2}_{wj}$, for $j=1,\dots,m$. For the internal
representation of the parameters in INLA, we transform the marginal variances
into log-precisions such that
$\boldsymbol{\theta}=(\phi_{11},\phi_{21},\dots,\phi_{mm},\log(1/\sigma^{2}_{w1}),\dots,\log(1/\sigma^{2}_{wm}))^{\prime}$
such that each of its components has $\mathbb{R}$ as its support, as required
by R-INLA. For more details about how to build an rgeneric() function,
vignette("rgeneric", package= "INLA") line in R provides documentation.

In order to define latent VAR$(1)$, we need to construct
‘inla.rgeneric.VAR1.model’, which returns several functions such as Q(),
graph(), mu(), log.prior(), initial(), log.norm.const(), etc.

  * •

The Q() function defines the precision matrix of the VAR$(1)$ process which is
derived in (8).

  * •

The graph() function defines the 0/1 representation of the precision matrix
indicating the zero and non-zero entries of the precision matrix as returned
by Q().

  * •

The mu() function returns $\boldsymbol{0}$.

  * •

We assign a Gaussian prior
$\mathcal{N}(\text{mean}=\mu,\text{variance}=\sigma^{2})$ on $\phi_{ij}$, for
$i=1,\dots,m$, $j=1,\dots,m$ and Gamma prior $\Gamma(.,a;b)$ (with mean $a/b$
and variance $a/b^{2}$) on the precisions, $1/\sigma^{2}_{wj}$, $j=1,\dots,m$.
Hence the joint prior for $\boldsymbol{\theta}$ becomes

$\pi(\boldsymbol{\theta})=\prod_{i=1}^{m^{2}}\mathcal{N}(\theta_{i},1;1)\times\prod_{j=m^{2}+1}^{m^{2}+m}\Gamma(\theta_{j};1,1)\times\exp(\theta_{j}),$
(23)

where for each $j=m^{2}+1,\dots,m^{2}+m$, $\exp(\theta_{j})$ is the Jacobian
for the change of variable from $1/\sigma^{2}_{j}$ to
$\theta_{j}=\log(1/\sigma^{2}_{wj})$. Hence the log.prior() returns the
logarithm of the joint prior as in (23).

  * •

We have defined the initial values of $\phi_{ij}=0.1$ and precisions
$1/\sigma^{2}_{wj}$ as $1$ for $i=1,\dots,m$ and $j=1,\dots,m$. The function
initial() contains the initial values of the parameter $\boldsymbol{\theta}$
in the internal scale.

  * •

The log.norm.const() denotes the logarithm of the normalizing constant from a
multivariate Gaussian distribution.

## References

  * Aït-Sahalia et al. (2010) Y. Aït-Sahalia, J. Fan, and D. Xiu.  High-frequency covariance estimates with noisy and asynchronous financial data.  _Journal of the American Statistical Association_ , 105(492):1504–1517, 2010.
  * Aitchison and Ho (1989) J. Aitchison and C. Ho.  The multivariate poisson-log normal distribution.  _Biometrika_ , 76(4):643–653, 1989.
  * Aktekin et al. (2020) T. Aktekin, N. G. Polson, and R. Soyer.  A family of multivariate non-gaussian time series models.  _Journal of Time Series Analysis_ , 41(5):691–721, 2020.
  * Andersen and Bollerslev (1998) T. G. Andersen and T. Bollerslev.  Answering the skeptics: Yes, standard volatility models do provide accurate forecasts.  _International Economic Review_ , pages 885–905, 1998.
  * Andersen et al. (2007) T. G. Andersen, T. Bollerslev, and F. X. Diebold.  Roughing it up: Including jump components in the measurement, modeling, and forecasting of return volatility.  _The Review of Economics and Statistics_ , 89(4):701–720, 2007.
  * Andersen et al. (2012) T. G. Andersen, D. Dobrev, and E. Schaumburg.  Jump-robust volatility estimation using nearest neighbor truncation.  _Journal of Econometrics_ , 169(1):75–93, 2012\.
  * Barndorff-Nielsen and Shephard (2004) O. E. Barndorff-Nielsen and N. Shephard.  Power and bipower variation with stochastic volatility and jumps.  _Journal of Financial Econometrics_ , 2(1):1–37, 2004.
  * Barndorff-Nielsen et al. (2009) O. E. Barndorff-Nielsen, P. R. Hansen, A. Lunde, and N. Shephard.  Realized kernels in practice: Trades and quotes.  _The Econometrics Journal_ , 2009.
  * Benjamin et al. (2003) M. A. Benjamin, R. A. Rigby, and D. M. Stasinopoulos.  Generalized autoregressive moving average models.  _Journal of the American Statistical Association_ , 98(461):214–223, 2003.
  * Cherian (1941) K. Cherian.  A bi-variate correlated gamma-type distribution function.  _The Journal of the Indian Mathematical Society_ , 5:133–144, 1941.
  * Chib and Winkelmann (2001) S. Chib and R. Winkelmann.  Markov chain monte carlo analysis of correlated count data.  _Journal of Business & Economic Statistics_, 19(4):428–435, 2001.
  * Chiogna et al. (2002) Chiogna, M. Carlo Gaetan, and C. Gaetan.  Dynamic generalized linear models with application to environmental epidemiology.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 51(4):453–468, 2002.
  * Cipollini et al. (2006) F. Cipollini, R. F. Engle III, and G. M. Gallo.  Vector multiplicative error models: representation and inference.  _Working Paper Series, NBER_ , 2006.
  * Cipollini et al. (2013) F. Cipollini, R. F. Engle, and G. M. Gallo.  Semiparametric vector MEM.  _Journal of Applied Econometrics_ , 28(7):1067–1086, 2013.
  * Corsi (2009) F. Corsi.  A simple approximate long-memory model of realized volatility.  _Journal of Financial Econometrics_ , 7(2):174–196, 2009.
  * Donelli et al. (2021) N. Donelli, S. Peluso, and A. Mira.  A Bayesian semiparametric vector multiplicative error model.  _Computational Statistics & Data Analysis_, 161:107242, 2021.
  * Durbin and Koopman (2012) J. Durbin and S. J. Koopman.  _Time Series Analysis by State Space Methods_ , volume 38.  OUP Oxford, 2012.
  * Engle and Gallo (2006) R. F. Engle and G. M. Gallo.  A multiple indicators model for volatility using intra-daily data.  _Journal of Econometrics_ , 131(1-2):3–27, 2006\.
  * Fahrmeir and Kaufmann (1987) L. Fahrmeir and H. Kaufmann.  Regression models for non-stationary categorical time series.  _Journal of Time Series Analysis_ , 8(2):147–160, 1987.
  * Fahrmeir et al. (1994) L. Fahrmeir, G. Tutz, W. Hennevogl, and E. Salem.  _Multivariate Statistical Modelling Based on Generalized Linear Models_ , volume 425.  Springer, 1994.
  * Gallego (2009) J. L. Gallego.  The exact likelihood function of a vector autoregressive moving average process.  _Statistics & Probability Letters_, 79(6):711–714, 2009.
  * Gamerman (1998) D. Gamerman.  Markov chain monte carlo for dynamic generalised linear models.  _Biometrika_ , 85(1):215–227, 1998.
  * Gaver and Lewis (1980) D. P. Gaver and P. Lewis.  First-order autoregressive gamma sequences and point processes.  _Advances in Applied Probability_ , 12(3):727–745, 1980.
  * Godolphin and Triantafyllopoulos (2006) E. Godolphin and K. Triantafyllopoulos.  Decomposition of time series models in state-space form.  _Computational Statistics & Data Analysis_, 50(9):2232–2246, 2006.
  * Gouriéroux and Jasiak (2006) C. Gouriéroux and J. Jasiak.  Autoregressive gamma processes.  _Journal of Forecasting_ , 25(2):129–152, 2006\.
  * Grunwald et al. (1993) G. K. Grunwald, P. Guttorp, and A. E. Raftery.  Prediction rules for exponential family state space models.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 55(4):937–943, 1993.
  * Hautsch et al. (2013) N. Hautsch, O. Okhrin, and A. Ristig.  Modeling time-varying dependencies between positive-valued high-frequency time series.  In _Copulae in Mathematical and Quantitative Finance_ , pages 115–127. Springer, 2013.
  * Heber et al. (2009) G. Heber, A. Lunde, N. Shephard, and K. Sheppard.  Oxford-Man Institute’s realized library, version 0.1, 2009.
  * Jacobs and Lewis (1977) P. Jacobs and P. Lewis.  A mixed autoregressive-moving average exponential sequence and point process (EARMA (1, 1)).  _Advances in Applied Probability_ , 9(1):87–104, 1977.
  * Kotz et al. (2004) S. Kotz, N. Balakrishnan, and N. L. Johnson.  _Continuous Multivariate Distributions, Volume 1: Models and Applications_ , volume 1.  John Wiley & Sons, 2004.
  * Lawrance and Lewis (1980) A. Lawrance and P. Lewis.  The exponential autoregressive-moving average earma (p, q) process.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 42(2):150–161, 1980.
  * Lawrance and Lewis (1981) A. Lawrance and P. Lewis.  A new autoregressive time series model in exponential variables (near (1)).  _Advances in Applied Probability_ , 13(4):826–845, 1981.
  * Lawrance and Lewis (1985) A. Lawrance and P. Lewis.  Modelling and residual analysis of nonlinear autoregressive time series in exponential variables.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 47(2):165–183, 1985.
  * Lindgren and Rue (2015) F. Lindgren and H. Rue.  Bayesian spatial modelling with R-INLA.  _Journal of statistical software_ , 63:1–25, 2015.
  * Lindsey and Lambert (1995) J. K. Lindsey and P. Lambert.  Dynamic generalized linear models and repeated measurements.  _Journal of Statistical Planning and Inference_ , 47(1-2):129–139, 1995.
  * Lütkepohl (2005) H. Lütkepohl.  _New Introduction to Multiple Time Series Analysis_.  Springer Science & Business Media, 2005.
  * Ma et al. (2008) J. Ma, K. M. Kockelman, and P. Damien.  A multivariate poisson-lognormal regression model for prediction of crash counts by severity, using Bayesian methods.  _Accident Analysis & Prevention_, 40(3):964–975, 2008.
  * Martins et al. (2013) T. G. Martins, D. Simpson, F. Lindgren, and H. Rue.  Bayesian computing with INLA: new features.  _Computational Statistics & Data Analysis_, 67:68–83, 2013.
  * Mathai and Moschopoulos (1991) A. M. Mathai and P. G. Moschopoulos.  On a multivariate gamma.  _Journal of Multivariate Analysis_ , 39(1):135–153, 1991.
  * Palmí-Perales et al. (2021) F. Palmí-Perales, V. Gómez-Rubio, and M. A. Martinez-Beneito.  Bayesian multivariate spatial models for lattice data with inla.  _Journal of Statistical Software_ , 98(2):1–29, 2021.  doi: 10.18637/jss.v098.i02.  URL https://www.jstatsoft.org/index.php/jss/article/view/v098i02.
  * Parkinson (1980) M. Parkinson.  The extreme value method for estimating the variance of the rate of return.  _Journal of Business_ , pages 61–65, 1980.
  * Peluso et al. (2014) S. Peluso, F. Corsi, and A. Mira.  A Bayesian high-frequency estimator of the multivariate covariance of noisy and asynchronous returns.  _Journal of Financial Econometrics_ , 13(3):665–697, 2014.
  * Petris (2010) G. Petris.  An R package for Dynamic Linear Models.  _Journal of Statistical Software_ , 36:1–16, 2010.
  * Prékopa and Szántai (1978) A. Prékopa and T. Szántai.  A new multivariate gamma distribution and its fitting to empirical streamflow data.  _Water Resources Research_ , 14(1):19–24, 1978\.
  * Raftery (1982) A. Raftery.  Generalized non-normal time series models.  _Time series analysis: Theory and practice_ , 1(0):621–640, 1982.
  * Ramabhadran (1951) V. Ramabhadran.  A multivariate gamma-type distribution.  _Sankhyā: The Indian Journal of Statistics_ , pages 45–46, 1951\.
  * Rue and Held (2005) H. Rue and L. Held.  _Gaussian Markov Random Fields: Theory and Applications_.  Chapman and Hall/CRC, 2005.
  * Rue et al. (2009) H. Rue, S. Martino, and N. Chopin.  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 71(2):319–392, 2009.
  * Ruiz-Cárdenas et al. (2012) R. Ruiz-Cárdenas, E. T. Krainski, and H. Rue.  Direct fitting of dynamic models using integrated nested laplace approximations—INLA.  _Computational Statistics & Data Analysis_, 56(6):1808–1828, 2012.
  * Serhiyenko (2015) V. Serhiyenko.  Dynamic modeling of multivariate counts-fitting, diagnostics, and applications.  _PhD Thesis, University of Connecticut_ , 2015.
  * Serhiyenko et al. (2015) V. Serhiyenko, N. Ravishanker, and R. Venkatesan.  Approximate Bayesian estimation for multivariate count time series models.  In _Ordered Data Analysis, Modeling and Health Research Methods_ , pages 155–167. Springer, 2015.
  * Serhiyenko et al. (2018) V. Serhiyenko, N. Ravishanker, and R. Venkatesan.  Multi-stage multivariate modeling of temporal patterns in prescription counts for competing drugs in a therapeutic category.  _Applied Stochastic Models in Business and Industry_ , 34(1):61–78, 2018.
  * Shephard and Pitt (1997) N. Shephard and M. K. Pitt.  Likelihood analysis of non-gaussian measurement time series.  _Biometrika_ , 84(3):653–667, 1997.
  * Spiegelhalter et al. (2002) D. J. Spiegelhalter, N. G. Best, B. P. Carlin, and A. Van Der Linde.  Bayesian measures of model complexity and fit.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 64(4):583–639, 2002.
  * Tsionas (2004) E. G. Tsionas.  Bayesian inference for multivariate gamma distributions.  _Statistics and Computing_ , 14(3):223–233, 2004\.
  * West and Harrison (2006) M. West and J. Harrison.  _Bayesian Forecasting and Dynamic Models_.  Springer Science & Business Media, 2006.
  * West et al. (1985) M. West, P. J. Harrison, and H. S. Migon.  Dynamic generalized linear models and bayesian forecasting.  _Journal of the American Statistical Association_ , 80(389):73–83, 1985.
  * Yue (2001) S. Yue.  A bivariate gamma distribution for use in multivariate flood frequency analysis.  _Hydrological Processes_ , 15(6):1033–1045, 2001.