# Moment based estimation for the multivariate COGARCH(1,1) process

Thiago do Rêgo Sousa Center for Mathematical Sciences, Technical University of
Munich, 85748 Garching, Boltzmannstr. 3, Germany, email<EMAIL_ADDRESS>Robert Stelzer Institute of Mathematical Finance, Ulm University, 89075 Ulm,
Helmholtzstr. 18, Germany, email<EMAIL_ADDRESS>
###### Abstract

For the multivariate COGARCH process, we obtain explicit expressions for the
second-order structure of the “squared returns” process observed on an
equidistant grid. Based on this, we present a generalized method of moments
estimator for its parameters. Under appropriate moment and strong mixing
conditions, we show that the resulting estimator is consistent and
asymptotically normal. Sufficient conditions for strong mixing, stationarity
and identifiability of the model parameters are discussed in detail. We
investigate the finite sample behavior of the estimator in a simulation study.

AMS 2010 Subject Classifications: primary: 62M05, 62M10, secondary: 60G51,
91B84.
Keywords: multivariate continuous time GARCH, second-order moment structure,
estimation, generalized method of moments, model identification, Lévy process

##  1 Introduction

The modeling of financial data has received much attention over the last
decades, where several models have been proposed for capturing its “stylized
facts”. Prominent models are the class of ARCH (autoregressive conditionally
heteroskedastic) and GARCH (generalized ARCH) processes introduced in Engle
(1982); Bollerslev (1986). They are able to capture most of these stylized
facts of financial data (see Cont (2001); Guillaume et al. (1997)). A special
feature of GARCH like processes is that they usually exhibit heavy tails even
if the driving noise is light tailed, a feature most other stochastic
volatility models do not have (Fasen et al. (2006)).

In many financial applications, it is most natural to model the price
evolution in continuous time, especially when dealing with high-frequency
data. The COGARCH process is a natural generalization of the discrete time
GARCH process to continuous time. It exhibits many “stylized features” of
financial time series and is well suited for modeling high-frequency data (see
Bayracı and Ünal (2014); Bibbona and Negri (2015); Haug et al. (2007);
Klüppelberg et al. (2011); Maller et al. (2008); Müller (2010)).

In many cases one needs to model the joint price of several financial assets
which exhibit a non-trivial dependence structure and therefore, multivariate
models are needed. The MUCOGARCH process introduced in Stelzer (2010) is a
multivariate extension of the COGARCH process. It combines the features of the
continuous time GARCH processes with the ones of the multivariate BEKK GARCH
process of Engle and Kroner (1995). It is a $d-$dimensional stochastic process
and it is defined as

$G_{t}=\int_{0}^{t}V_{s-}^{1/2}{\rm d}L_{s},\quad t\geq 0,$ (1.1)

where $L$ is an ${\mathbb{R}}^{d}$-valued Lévy process with non-zero Lévy
measure and càdlàg sample paths. The matrix-valued volatility process
$(V_{s})_{s\in{\mathbb{R}}^{+}}$ depends on a parameter
$\theta\in\Theta\subset{\mathbb{R}}^{q}$, it is predictable and its randomness
depends only on $L$. We assume that we have a sample of size $n$ of the log-
price process (1.1) with true parameter $\theta_{0}\in\Theta$ observed on a
fixed grid of size $\Delta>0$, and compute the log returns

${\bm{G}}_{i}=\int_{(i-1)\Delta}^{i\Delta}V_{s-}^{1/2}{\rm d}L_{s},\quad
i=1,\dots,n.$ (1.2)

Therefore, an important question is how to estimate the true parameter
$\theta_{0}$ based on observations $({\bm{G}}_{i})_{i=1}^{n}$. In the
univariate case, several methods have been proposed to estimate the parameters
of the COGARCH process (Bayracı and Ünal (2014); Bibbona and Negri (2015); do
Rêgo Sousa et al. (2019); Haug et al. (2007); Maller et al. (2008)). All these
methods rely on the fact that the COGARCH process is, under certain regularity
conditions, ergodic and strongly mixing.

In the univariate case, Fasen (2010) proved geometric ergodicity results for
the COGARCH process (in fact, their results apply to a wider class of Lévy
driven models). Recently, Stelzer and Vestweber (2019) derived sufficient
conditions for the existence of a unique stationary distribution, for the
geometric ergodicity, and for the finiteness of moments of the stationary
distribution in the MUCOGARCH process. These results imply ergodicity and
strong mixing of the log-price process $({\bm{G}}_{i})_{i=1}^{\infty}$, thus
paving the way for statistical inference. We will use their results to apply
the generalized method of moments (GMM) for estimating the parameters of the
MUCOGARCH process. To this end we compute the second-order structure of the
squared returns in closed form, under appropriate assumptions.

Consistency and asymptotic normality of the GMM estimator is obtained under
standard assumptions of strong mixing, existence of moments of the MUCOGARCH
volatility process and model identifiability. Thus we discuss sufficient
conditions, easily checkable for given parameter spaces ensuring strong mixing
and existence of relevant moments.

The identifiability question is rather delicate, since the formulae for the
second-order structure of the log-price returns involve operators which are
not invertible and, therefore, the strategy used for showing identifiability
as used in the one-dimensional COGARCH process cannot be generalised. In the
end we can establish identifiability conditions that are not overly
restrictive and easy to use.

Our paper is organized as follows. In Section 2, we fix the notation and
briefly introduce Lévy processes. In Section 3 we define the MUCOGARCH
process, and obtain in Section 4 its second-order structure. Section 5
introduces the GMM estimator and discusses sufficient conditions for
stationarity, strong mixing and identifiability of the model. In Section 6, we
study the finite sample behavior of the estimators in a simulation study.
Finally, Section A presents the proofs for the results of Sections 3 and 4.

##  2 Preliminaries

###  2.1 Notation

Denote the set of non-negative real numbers by ${\mathbb{R}}^{+}$. For
$z\in{\mathbb{C}},\Re(z)$ and $\Im(z)$ denote the real and imaginary part,
respectively. We denote by $M_{m,d}({\mathbb{R}})$, the set of real $m\times
d$ matrices and write $M_{d}({\mathbb{R}})$ for $M_{d,d}({\mathbb{R}})$. The
group of invertible $d\times d$ matrices is denoted by $GL_{d}({\mathbb{R}})$,
the linear subspace of symmetric matrices by ${\mathbb{S}}_{d}$, the (closed)
positive semidefinite cone by ${\mathbb{S}}^{+}_{d}$ and the (open) positive
definite cone by ${\mathbb{S}}^{++}_{d}$. We write $I_{d}$ for the $d\times d$
identity matrix. The tensor (Kronecker) product of two matrices $A,B$ is
written as $A\otimes B$. The $\operatorname{vec}$ operator denotes the well-
known vectorization operator that maps the set of $d\times d$ matrices to
${\mathbb{R}}^{d^{2}}$ by stacking the columns of the matrices below one
another. Similarly, $\operatorname{vech}$ stacks the entries on and below the
main diagonal of a square matrix. For more information regarding the tensor
product, $\operatorname{vec}$ and $\operatorname{vech}$ operators we refer to
Bernstein (2009); Horn and Johnson (1991). The spectrum of a square matrix is
denoted by $\sigma(\cdot)$. Finally, $A^{*}$ denotes the transpose of a matrix
$A\in M_{m,d}({\mathbb{R}})$ and $A_{(i,j)}$ denotes the entry in the $i$th
line and $j$th column of $A$. Arbitrary norms of vectors or matrices are
denoted by $\|\cdot\|$ in which case it is irrelevant which particular norm is
used. The norm $\|\cdot\|_{2}$ denotes the operator norm on
$M_{d^{2}}(\mathbb{R})$ associated with the usual Euclidean norm. The symbol
$c$ stands for any positive constant, whose value may change from line to
line, but is not of particular interest.

Additionally, we employ an intuitive notation with respect to (stochastic)
integration with matrix-valued integrators, referring to any of the standard
texts (for example, Protter (2005)) for a comprehensive treatment of the
theory of stochastic integration. Let $(A_{t})_{t\in{\mathbb{R}}^{+}}$ in
$M_{m,d}({\mathbb{R}})$ and $(B_{t})_{t\in{\mathbb{R}}^{+}}$ in
$M_{r,u}({\mathbb{R}})$ be càdlàg and adapted processes and
$(L_{t})_{t\in{\mathbb{R}}^{+}}$ in $M_{d,r}({\mathbb{R}})$ be a
semimartingale. We then denote by $\int_{0}^{t}A_{s-}{\rm d}L_{s}B_{s-}$ the
matrix $C_{t}\in M_{m,u}({\mathbb{R}})$ which has $ij$-th entry
$\sum_{k=1}^{d}\sum_{l=1}^{r}\int_{0}^{t}A_{ik,s-}B_{lj,s-}\mathrm{d}L_{kl,s}$.
If $(X_{t})_{t\in{\mathbb{R}}^{+}}$ is a semimartingale in ${\mathbb{R}}^{m}$
and $(Y_{t})_{t\in{\mathbb{R}}^{+}}$ one in ${\mathbb{R}}^{d}$, then the
quadratic variation $([X,Y]_{t})_{t\in{\mathbb{R}}^{+}}$ is defined as the
finite variation process in $M_{m,d}({\mathbb{R}})$ with $ij$-th entry
$[X_{i},Y_{j}]_{t}$ for $t\in{\mathbb{R}}^{+}$, $i=1,\dots,m$ and
$j=1,\dots,d$. We also refer to Lemma 2.2 in Behme (2012) for a collection of
basic properties related to integration with matrix-valued integrators.
Lastly, let ${\bm{Q}}:M_{d^{2}}({\mathbb{R}})\mapsto M_{d^{2}}({\mathbb{R}})$
be the linear map defined by

$({\bm{Q}}X)_{(k-1)d+l,(p-1)d+q}=X_{(k-1)d+p,(l-1)d+q}\quad\text{ for all
}k,l,p,q={1,\dots,d},$

which has the property that
${\bm{Q}}(\operatorname{vec}(X)\operatorname{vec}(Z)^{T})=X\otimes Z$ for all
$X,Z\in{\mathbb{S}}_{d}$ ((Pigorsch and Stelzer, 2009b, Theorem 4.3)). Let
$K_{d}$ be the commutation matrix characterized by
$K_{d}\operatorname{vec}(A)=\operatorname{vec}(A^{*})$ for all $A\in
M_{d}({\mathbb{R}})$ (see Magnus and Neudecker (1979) for more details).
Define ${\mathcal{Q}}\in M_{d^{4}}({\mathbb{R}})$ as the matrix associated
with the linear map
$\operatorname{vec}\circ{\bm{Q}}\circ\operatorname{vec}^{-1}$ on
${\mathbb{R}}^{d^{4}}$, and ${\mathcal{K}}_{d}\in M_{d^{4}}({\mathbb{R}})$ as
the matrix associated with the linear map
$\operatorname{vec}(K_{d}\operatorname{vec}^{-1}(x))$ for
$x\in{\mathbb{R}}^{d^{4}}$.

###  2.2 Lévy processes

A Lévy process $L=(L_{t})_{t\in{\mathbb{R}}^{+}}$ in ${\mathbb{R}}^{d}$ is
characterized by its characteristic function in Lévy-Khintchine form
${\mathbb{E}}e^{i\langle u,L_{t}\rangle}$ = $\exp\\{t\psi_{L}(u)\\}$ for
$t\in{\mathbb{R}}^{+}$ with

$\psi_{L}(u)=i\langle\gamma_{L},u\rangle-\frac{1}{2}\langle
u,\Gamma_{L}u\rangle+\int_{\mathbb{R}^{d}}\big{(}e^{i\langle
u,x\rangle}-1-i\langle u,x\rangle I_{[0,1]}(\|x\|)\big{)}\nu_{L}(dx),\quad
u\in\mathbb{R}^{d},$

where $\gamma_{L}\in\mathbb{R}^{d},\Gamma_{L}\in\mathbb{S}_{d}^{+}$ and the
Lévy measure $\nu_{L}$ is a non-zero measure on $\mathbb{R}^{d}$ satisfying
$\nu_{L}(\\{0\\})=0$ and $\int_{\mathbb{R}^{d}}\left(\|x\|^{2}\wedge
1\right)\nu_{L}(dx)<\infty$. We assume w.l.o.g. $L$ to have càdlàg paths. The
discontinuous part of the quadratic variation of $L$ is denoted by
$([L,L]_{t}^{{\mathfrak{d}}})_{t\in{\mathbb{R}}^{+}}$ and it is also a Lévy
process. It has finite variation, zero drift and Lévy measure
$\nu_{[L,L]^{{\mathfrak{d}}}}(B)=\int_{\mathbb{R}^{d}}I_{B}\left(xx^{*}\right)\nu_{L}(dx)$
for all Borel sets $B\subseteq\mathbb{S}_{d}$. For more details on Lévy
processes we refer to Applebaum (2009); Sato (1999).

##  3 The MUCOGARCH process

Throughout, we assume that all random variables and processes are defined on a
given filtered probability space
$(\Omega,{\mathcal{F}},P,({\mathcal{F}}_{t})_{t\in T})$, with $T={\mathbb{N}}$
in the discrete-time case and $T={\mathbb{R}}^{+}$ in the continuous-time one.
In the continuous-time setting, we assume the usual conditions (complete,
right-continuous filtration) to be satisfied. We can now recall the definition
of the MUCOGARCH process.

######  Definition 3.1 (MUCOGARCH(1,1) - (Stelzer, 2010, Definition 3.1)).

Let $L$ be an ${\mathbb{R}}^{d}$-valued Lévy process, $A,B\in
M_{d}({\mathbb{R}})$ and $C\in{\mathbb{S}}_{d}^{++}$. The process
$G=(G_{t})_{t\in{\mathbb{R}}^{+}}$ solving

$\displaystyle{\rm d}G_{t}$ $\displaystyle=$ $\displaystyle V_{t-}^{1/2}{\rm
d}L_{t}$ (3.1) $\displaystyle V_{t}$ $\displaystyle=$ $\displaystyle C+Y_{t}$
(3.2) $\displaystyle{\rm d}Y_{t}$ $\displaystyle=$
$\displaystyle(BY_{t-}+Y_{t-}B^{\ast}){\rm d}t+AV_{t-}^{1/2}{\rm
d}[L,L]_{t}^{{\mathfrak{d}}}V_{t-}^{1/2}A^{\ast}$ (3.3)

with initial values $G_{0}$ in ${\mathbb{R}}^{d}$ and $Y_{0}$ in
$\mathbb{S}_{d}^{+}({\mathbb{R}})$ is called a MUCOGARCH(1,1) process. The
process $Y=(Y_{t})_{t\in{\mathbb{R}}^{+}}$ is called a MUCOGARCH(1,1)
volatility process. Hereafter we will always write MUCOGARCH for short.

The interpretation of the model parameters $B$ and $C$ is the following. If
$\sigma(B)\in\\{z\in{\mathbb{C}}:\Re(z)<0\\}$, the process $V$, as long as no
jump occurs, “mean reverts” to the level $C$ at matrix exponential rate given
by $B$. Since all jumps are positive semidefinite, $C$ is not a mean level
but, instead, a lower bound for $V$.

By (Stelzer, 2010, Theorems 3.2 and 4.4), the MUCOGARCH process is well-
defined, the solution $(Y_{t})_{t\in{\mathbb{R}}^{+}}$ is locally bounded and
of finite variation. Additionally, the process
$(G_{t},Y_{t})_{t\in{\mathbb{R}}^{+}}$ and its volatility process
$(Y_{t})_{t\in{\mathbb{R}}^{+}}$ are time homogeneous strong Markov processes
on ${\mathbb{R}}^{d}\times{\mathbb{S}}^{+}_{d}$ and ${\mathbb{S}}^{+}_{d}$,
respectively.

Since the price process $(G_{t})_{t\in{\mathbb{R}}^{+}}$ in (3.1) is defined
in terms of the Lévy process $L$ and $(Y_{t})_{t\in{\mathbb{R}}_{+}}$, the
existence of its moments is closely related to the existence of moments of $L$
and the stationary distribution of $(Y_{t})_{t\in{\mathbb{R}}_{+}}$.

######  Lemma 3.2.

Suppose that ${\mathbb{E}}\|Y_{0}\|^{p}<\infty$ and
${\mathbb{E}}\|L_{1}\|^{2p}<\infty$ for some $p\geq 1$. Then:

  * (a)

${\mathbb{E}}\|Y_{t}\|^{p}<\infty$ for all $t\in{\mathbb{R}}^{+}$ and
$t\mapsto{\mathbb{E}}\|Y_{t}\|^{p}$ is locally bounded.

  * (b)

${\mathbb{E}}\|G_{t}\|^{2p}<\infty$ for all $t\in{\mathbb{R}}^{+}$ and
$t\mapsto{\mathbb{E}}\|G_{t}\|^{2p}$ is locally bounded.

##  4 Second-order structure of “squared returns”

In this section, we derive the second-order structure of the MUCOGARCH
“squared returns” process
$({\bm{G}}_{i}{\bm{G}}_{i}^{\ast})_{i\in{\mathbb{N}}}$ defined in terms of
(1.2), which will be used in Section 5 to estimate the parameters $A,B$ and
$C$ of the MUCOGARCH process. The proofs are postponed to Section A. We group
the needed assumptions as follows.

######  Assumptions a (Lévy process).

  1. $(a.1)$

${\mathbb{E}}L_{1}=0$.

  2. $(a.2)$

${\rm var}(L_{1})=(\sigma_{W}+\sigma_{L})I_{d}$, with $\sigma_{W}\geq 0$ and
$\sigma_{L}>0$.

  3. $(a.3)$

$\int_{{\mathbb{R}}^{d}}x_{i}x_{j}x_{k}\,\nu_{L}(dx)=0,\quad\text{for all
}i,j,k\in\\{1,\dots,d\\}.$

  4. $(a.4)$

${\mathbb{E}}\|L_{1}\|^{4}<\infty$.

  5. $(a.5)$

There exists a constant $\rho_{L}>0$ such that

${\mathbb{E}}[\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}),\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}})^{\ast}]_{1}^{\mathfrak{d}}=\rho_{L}(I_{d^{2}}+K_{d}+\operatorname{vec}(I_{d})\operatorname{vec}(I_{d})^{\ast}).$

  6. $(a.6)$

${\mathbb{E}}\|L_{1}\|^{8}<\infty$.

######  Assumptions b (Parameters).

  1. $(b.1)$

$A\in GL_{d}({\mathbb{R}})$.

  2. $(b.2)$

The matrices ${\mathcal{B}}$ and ${\mathcal{C}}$ defined below satisfy
$\sigma({\mathcal{B}}),\sigma({\mathcal{C}})\in\\{z\in{\mathbb{C}}:\Re(z)<0\\}$.

$\displaystyle{\mathcal{B}}$ $\displaystyle:=$ $\displaystyle B\otimes
I+I\otimes B+\sigma_{L}(A\otimes A)$ (4.1) $\displaystyle{\mathcal{C}}$
$\displaystyle:=$ $\displaystyle{\mathcal{B}}\otimes
I_{d^{2}}+I_{d^{2}}\otimes{\mathcal{B}}+{\mathcal{A}}{\mathcal{R}},$

where ${\mathcal{A}}=(A\otimes A)\otimes(A\otimes A)$,
${\mathcal{R}}=\rho_{L}({\mathcal{Q}}+{\mathcal{K}}_{d}{\mathcal{Q}}+I_{d^{4}})$,
and ${\mathcal{K}}_{d}$ and ${\mathcal{Q}}$ as in Section 2.1.

######  Assumption c (MUCOGARCH volatility).

  1. $(c.1)$

$(Y_{t})_{t\in{\mathbb{R}}_{+}}$ is a second-order stationary MUCOGARCH
volatility process.

  2. $(c.2)$

$(Y_{t})_{t\in{\mathbb{R}}_{+}}$ is a stationary MUCOGARCH volatility process
and its stationary distribution satisfies ${\mathbb{E}}\|Y_{0}\|^{4}<\infty$.

Sufficient conditions for Assumption c are given in (Stelzer, 2010, Theorem
4.5). Note that $(c.2)$ implies $(c.1)$. We recall now the expressions for the
second-order structure of the process $Y$ and of the log-price returns process
$({\bm{G}}_{i})_{i\in{\mathbb{N}}}$. First, for a second-order stationary
${\mathbb{R}}^{d}$-valued process, its autocovariance function ${\rm
acov}_{X}:{\mathbb{R}}\mapsto M_{d}({\mathbb{R}})$ is denoted by
$\operatorname{acov}_{X}(h)=\operatorname{cov}\left(X_{h},X_{0}\right)={\mathbb{E}}\left(X_{h}X_{0}^{*}\right)-{\mathbb{E}}\left(X_{0}\right)E\left(X_{0}\right)^{*}$
for $h\geq 0$ and by
$\operatorname{acov}_{X}(h)=\left(\operatorname{acov}_{X}(-h)\right)^{*}$ for
$h<0$. For matrix-valued processes $(Z_{t})_{t\in{\mathbb{R}}}$, we set ${\rm
acov}_{Z}={\rm acov}_{\operatorname{vec}(Z)}$.

######  Proposition 4.1 ((Stelzer, 2010, Theorems 4.8, 4.11, Corollary 4.19
and Proposition 5.2)).

If Assumptions $(a.1)$-$(a.5)$, $(b.2)$ and $(c.1)$ hold, then

$\displaystyle{\mathbb{E}}(\operatorname{vec}(Y_{0}))$ $\displaystyle=$
$\displaystyle-\sigma_{L}{\mathcal{B}}^{-1}(A\otimes A)\operatorname{vec}(C)$
(4.2) $\displaystyle{\rm var}(\operatorname{vec}(Y_{0}))$ $\displaystyle=$
$\displaystyle{\rm
var}(\operatorname{vec}(V_{0}))=-{\mathcal{C}}^{-1}\big{[}\big{(}\sigma_{L}^{2}{\mathcal{C}}({\mathcal{B}}^{-1}\otimes{\mathcal{B}}^{-1}){\mathcal{A}}+{\mathcal{A}}{\mathcal{R}}\big{)}(\operatorname{vec}(C)\otimes\operatorname{vec}(C))$
$\displaystyle\quad\quad+\big{(}\sigma_{L}(A\otimes A)\otimes
I_{d^{2}}+{\mathcal{A}}{\mathcal{R}}\big{)}\operatorname{vec}(C)\otimes{\mathbb{E}}(\operatorname{vec}(Y_{0}))$
$\displaystyle\quad\quad+\big{(}\sigma_{L}I_{d^{2}}\otimes(A\otimes
A)+{\mathcal{A}}{\mathcal{R}}\big{)}{\mathbb{E}}(\operatorname{vec}(Y_{0}))\otimes\operatorname{vec}(C)\big{]}$
$\displaystyle{\rm acov}_{Y}(h)$ $\displaystyle=$ $\displaystyle{\rm
acov}_{V}(h)=e^{{\mathcal{B}}h}{\rm var}(\operatorname{vec}(Y_{0}))$
$\displaystyle{\mathbb{E}}({\bm{G}}_{1})$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle{\rm var}({\bm{G}}_{1})$ $\displaystyle=$
$\displaystyle(\sigma_{L}+\sigma_{W})\Delta{\mathbb{E}}(C+Y_{0})$ (4.3)
$\displaystyle{\rm acov}_{\bm{G}}(h)$ $\displaystyle=$ $\displaystyle
0\quad\text{for all }h\in{\mathbb{Z}}\backslash\\{0\\}.$

Based on Lemma 3.2 and Proposition 4.1, we obtain now the second-order
properties of the MUCOGARCH process.

######  Lemma 4.2.

If Assumptions a,b and c hold, then

$\begin{split}{\rm acov}_{\bm{G}\bm{G}^{\ast}}(h)&=e^{{\mathcal{B}}\Delta
h}{\mathcal{B}}^{-1}(I_{d^{2}}-e^{-{\mathcal{B}}\Delta})(\sigma_{L}+\sigma_{W}){\rm
var}(\operatorname{vec}(V_{0}))\\\
&\quad\quad\times(e^{{\mathcal{B}}^{\ast}\Delta}-I_{d^{2}})[(\sigma_{W}+\sigma_{L})({\mathcal{B}}^{\ast})^{-1}-2((A\otimes
A)^{\ast})^{-1}],\quad h\in{\mathbb{N}},\end{split}$ (4.4)
$\displaystyle{\mathbb{E}}\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})^{\ast}$
(4.5)
$\displaystyle=\Delta\rho_{L}\big{(}({\bm{Q}}+K_{d}{\bm{Q}}+I_{d^{2}})({\mathbb{E}}\operatorname{vec}(V_{0})\operatorname{vec}(V_{0})^{\ast})\big{)}\times$
$\displaystyle\,\,\,\,\,\,(I_{d^{2}}+K_{d}){\bm{Q}}(D^{\ast})(I_{d^{2}}+K_{d})+D+D^{\ast},$

with,

$D:=(\sigma_{L}+\sigma_{W})\Big{(}\frac{1}{2}(\sigma_{L}+\sigma_{W})\Delta^{2}{\mathbb{E}}\operatorname{vec}(V_{0}){\mathbb{E}}\operatorname{vec}(V_{0})^{\ast}+{\rm
var}(\operatorname{vec}(V_{0}))\tilde{{\mathcal{B}}}\Big{)}$ (4.6)
$\tilde{{\mathcal{B}}}:=\big{[}({\mathcal{B}}^{\ast})^{-1}(e^{{\mathcal{B}}^{\ast}\Delta}-I_{d^{2}})-I_{d^{2}}\Delta\big{]}\big{[}(\sigma_{W}+\sigma_{L})({\mathcal{B}}^{\ast})^{-1}-2((A\otimes
A)^{\ast})^{-1}\big{]}$ (4.7)

######  Remark 4.3.

If the Lévy process $L$ has paths of finite variation, then Lemma 4.2 holds
without the moment assumptions $(a.6)$ and $(c.2)$. This is because
expectations involving stochastic integrals with finite variation Lévy
integrators can be computed by using the compensation formula (see Remark
A.2). In the following, we will define the moment based estimator for
MUCOGARCH processes driven by general Lévy processes (without path
restrictions). Only in Section 5.5 we will give a consistency result that
distinguishes between Lévy process with paths of finite and infinite
variation.

Next, we define an estimator for the parameters $A,B$ and $C$, which basically
consists of comparing the sample moments to the model moments.

##  5 Moment based estimation of the MUCOGARCH process

In this section, we consider the matrices $A_{\theta},B_{\theta}\in
M_{d}({\mathbb{R}})$ and $C_{\theta}\in\mathbb{S}_{d}^{++}$ from Definition
3.1 as depending on a parameter $\theta\in\Theta\subset{\mathbb{R}}^{q}$ for
$q\in{\mathbb{N}}$.

The data used for estimation is an equidistant sample of $d$-dimensional log-
prices $({\bm{G}}_{i})_{i=1}^{n}$ as defined in (1.2) with true parameter
$\theta_{0}\in\Theta$. We assume that the true $\sigma_{L},\sigma_{W}$ and
$\rho_{L}$ as used in Assumptions $(a.2)$ and $(a.5)$ are known. These
assumptions are not very restrictive and are comparable to assuming iid
standard normal noise in the discrete time multivariate GARCH process, which
is very common (Francq and Zakoïan, 2019, eq. (10.6)).

###  5.1 Generalized Method of Moments (GMM) estimator

In order to estimate the parameter $\theta_{0}\in\Theta$, we compare the
sample moments (based on a sample of log-prices) to the model moments (based
on the expressions (4.3), (4.4) and (4.5), provided they are well defined).
More specifically, based on the observations $({\bm{G}}_{i})_{i=1}^{n}$ and a
fixed $r<n$, the sample moments are defined as

$\hat{k}_{n,r}=\frac{1}{n}\sum_{i=1}^{n-r}D_{i}=\frac{1}{n}\sum_{i=1}^{n-r}\begin{pmatrix}\operatorname{vec}(\bm{G}_{i}\bm{G}_{i}^{\ast})\\\
\operatorname{vec}(\operatorname{vec}(\bm{G}_{i}\bm{G}_{i}^{\ast})\operatorname{vec}(\bm{G}_{i}\bm{G}_{i}^{\ast})^{*})\\\
\vdots\\\
\operatorname{vec}(\operatorname{vec}(\bm{G}_{i}\bm{G}_{i}^{\ast})\operatorname{vec}(\bm{G}_{i+r}\bm{G}_{i+r}^{\ast})^{*})\end{pmatrix}.$
(5.1)

The used number of lags of the true autocovariance function $r$ needs to be
chosen in such a way that the model parameters are identifiable and also to
ensure a good fit of the autocovariance structure to the data. For each
$\theta\in\Theta$, let

$k_{\theta,r}=\begin{pmatrix}{\mathbb{E}}_{\theta}\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\\\
{\mathbb{E}}_{\theta}\operatorname{vec}(\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})^{*})\\\
\vdots\\\
{\mathbb{E}}_{\theta}\operatorname{vec}(\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1+r}\bm{G}_{1+r}^{\ast})^{*})\end{pmatrix},$
(5.2)

where the expectations are explicitly given by (4.3), (4.4) and (4.5) by
replacing $A,B$ and $C$ by $A_{\theta},B_{\theta}$ and $C_{\theta}$,
respectively. Then, the GMM estimator of $\theta_{0}$ is given by

$\hat{\theta}_{n}=\operatorname*{arg\,min}_{\theta\in\Theta}\Big{\\{}(\hat{k}_{n,r}-k_{\theta,r})^{T}\Omega(\hat{k}_{n,r}-k_{\theta,r})\Big{\\}},$
(5.3)

where $\Omega$ is a positive definite weight matrix. The matrix $\Omega$ may
be depend on the data but should converge in probability to a positive
definite matrix of constants.

###  5.2 Asymptotic properties: general case

Additionally to Assumptions a, b and c we need assumptions for proving
consistency and asymptotic normality of $\hat{\theta}_{n}$. These are mainly
related to identifiability of the model parameters, stationarity, strong
mixing and existence of certain moments of
$({\bm{G}}_{i})_{i\in{\mathbb{N}}}$.

######  Assumptions d (Parameter space and log-price process).

  1. $(d.1)$

The parameter space $\Theta$ is a compact subset of ${\mathbb{R}}^{q}$.

  2. $(d.2)$

The true parameter $\theta_{0}$ lies in the interior of $\Theta$.

  3. $(d.3)$

[Identifiability]. Let $r>1$ be fixed. For any
$\theta\not=\tilde{\theta}\in\Theta$ we have
$k_{\theta,r}\not=k_{\tilde{\theta},r}$.

  4. $(d.4)$

The map $\theta\mapsto(A_{\theta},B_{\theta},C_{\theta})$ is continuously
differentiable.

  5. $(d.5)$

The sequence $({\bm{G}}_{i})_{i\in{\mathbb{N}}}$ is strictly stationary and
exponentially $\alpha$-mixing.

######  Assumption e (Moments).

  1. $(e.1)$

There exists a positive constant $\delta>0$ such that
${\mathbb{E}}\|{\bm{G}}_{1}\|^{8+\delta}<\infty$.

Assumption e can be written in terms of moments of $L$ and $Y_{0}$ (see Lemma
3.2). We are now ready to state the strong consistency of the empirical
moments in (5.1).

######  Lemma 5.1.

If Assumptions a, b, c and $(d.5)$ hold, then
$\hat{k}_{n,r}\stackrel{{\scriptstyle\rm a.s.}}{{\rightarrow}}k_{\theta_{0}}$
as ${n\to\infty}$.

###### Proof.

It follows from $(d.5)$ that the log-price process
$({\bm{G}}_{i})_{i\in{\mathbb{N}}}$ is ergodic and since both
${\mathbb{E}}\|\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\|$ and
${\mathbb{E}}\|\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1+h}\bm{G}_{1+h}^{\ast})^{*}\|$
are finite (Lemma 3.2 with $p=2$ under $(a.4)$ and $(c.1)$), we can apply
Birkhoff’s ergodic theorem ((Krengel, 1985, Theorem 4.4)) to conclude the
result. ∎

Next, we state the weak consistency property of the GMM estimator.

######  Theorem 5.2.

If Assumptions a,b, c, $(d.1)$, $(d.3)$-$(d.5)$ hold, then the GMM estimator
defined in (5.3) is weakly consistent.

###### Proof.

We check Assumptions 1.1-1.3 in Mátyás (1999) that ensure weak consistency of
the GMM estimator in (5.3). Assumption 1.1 is satisfied due to our
identifiability condition $(d.3)$. It follows from (5.3) combined with Lemma
5.1 that

$\sup_{\theta\in\Theta}\|\hat{k}_{n,r}-k_{\theta,r}-(k_{\theta_{0},r}-k_{\theta,r})\|=\|\hat{k}_{n,r}-k_{\theta_{0},r}\|\stackrel{{\scriptstyle\rm
a.s.}}{{\rightarrow}}0,\quad{n\to\infty},$

which is Assumption 1.2 of Mátyás (1999). Since the weight matrix $\Omega$ in
(5.3) is non-random, their Assumption 1.3 is automatically satisfied,
completing the proof. ∎

In order to prove asymptotic normality of the GMM estimator, we need some
auxiliary results.

######  Lemma 5.3.

If Assumptions a, b, c, $(d.1)$ and $(d.4)$ hold, then the map $\Theta\mapsto
k_{\theta,r}$ in (5.2) is continuously differentiable.

###### Proof.

The the map $\Theta\mapsto k_{\theta,r}$ depends on the moments given in
(4.3), (4.4) and (4.5). These moments are given in terms of products and
Kronecker products involving the quantities $A_{\theta}$, $A_{\theta}^{-1}$,
${\mathcal{B}}_{\theta}$, ${\mathcal{B}}_{\theta}^{-1}$,
$e^{-\alpha{\mathcal{B}}_{\theta}}$, $\alpha>0$, $C_{\theta}$,
${\mathcal{C}}_{\theta}$ and ${\mathcal{C}}_{\theta}^{-1}$. From $(d.4)$ we
obtain the continuous differentiability of
${\mathcal{B}}_{\theta}$,${\mathcal{B}}_{\theta}^{-1}$,
${\mathcal{C}}_{\theta}$,${\mathcal{C}}_{\theta}^{-1}$ and $A_{\theta}^{-1}$
on $\Theta$. Let $i\in\\{1,\dots,q\\}$ be fixed. According to (2.1) in Wilcox
(1967), the matrix exponential is differentiable and

$\frac{\partial}{\partial\theta_{i}}e^{-\alpha{\mathcal{B}}_{\theta}}=-\int_{0}^{\alpha}e^{-(\alpha-u){\mathcal{B}}_{\theta}}\bigg{(}\frac{\partial}{\partial\theta_{i}}{\mathcal{B}}_{\theta}\bigg{)}e^{-u{\mathcal{B}}_{\theta}}{\rm
d}u.$ (5.4)

Using the definition of ${\mathcal{B}}_{\theta}$ in (4.1) combined with
$(d.1)$ and $(d.4)$ gives

$\sup_{\theta\in\Theta}\|{\mathcal{B}}_{\theta}\|\leq
2\bigg{(}\sup_{\theta\in\Theta}\|B_{\theta}\|\bigg{)}\|I_{d}\|+\sigma_{L}\bigg{(}\sup_{\theta\in\Theta}\|A_{\theta}\|^{2}\bigg{)}<\infty.$

Additionally, an application of the chain rule to
$\frac{\partial}{\partial\theta_{i}}{\mathcal{B}}_{\theta}$ combined with
$(d.1)$ and $(d.4)$ gives
$\sup_{\theta\in\Theta}\|\frac{\partial}{\partial\theta_{i}}{\mathcal{B}}_{\theta}\|<\infty$
and, therefore,

$\sup_{\theta\in\Theta}\bigg{\|}e^{-(\alpha-u){\mathcal{B}}_{\theta}}\bigg{(}\frac{\partial}{\partial\theta_{i}}{\mathcal{B}}_{\theta}\bigg{)}e^{-u{\mathcal{B}}_{\theta}}\bigg{\|}\leq\sup_{\theta\in\Theta}e^{(|\alpha-u|+|u|)\|{\mathcal{B}}_{\theta}\|}\bigg{(}\sup_{\theta\in\Theta}\bigg{\|}\frac{\partial}{\partial\theta_{i}}{\mathcal{B}}_{\theta}\bigg{\|}\bigg{)},\quad
u\in[0,\alpha].$ (5.5)

Thus, the continuous differentiability of the map in (5.4) follows by
dominated convergence with dominating function as in (5.5). Another
application of the chain rule shows that the map $\theta\mapsto k_{\theta,r}$
is continuously differentiable on $\Theta$. ∎

######  Lemma 5.4.

Assume that Assumptions a,b, c, $(d.5)$ and $(e.1)$ hold and let

$\Sigma_{\theta_{0}}={\mathbb{E}}(F_{1}F_{1}^{\ast})+\sum_{i=1}^{\infty}{\mathbb{E}}\big{\\{}(F_{1}F_{1+i}^{\ast})+{\mathbb{E}}(F_{1}F_{1+i}^{\ast})^{\ast}\big{\\}}$
(5.6)

with $F_{i}=D_{i}-k_{\theta_{0},r}$ and $D_{i}$ as defined in (5.1). Then for
$r\in{\mathbb{N}}_{0}$

$\sqrt{n}(\hat{k}_{n,r}-k_{\theta_{0},r})\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\Sigma_{\theta_{0}}),\quad
n\rightarrow\infty.$

###### Proof.

For the asymptotic normality of (5.1) we use the Cramér-Wold device and show
that

$\sqrt{n}\Big{(}\frac{1}{n}\sum_{i=1}^{n-r}\lambda^{\ast}F_{i}\Big{)}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\lambda^{\ast}\Sigma_{\theta_{0}}\lambda),\quad{n\to\infty},$

for all vectors $\lambda\in{\mathbb{R}}^{d^{2}+(r+1)d^{4}}$. Denote by
$\alpha_{{\bm{G}}}$ the mixing coefficients of
$({\bm{G}}_{i})_{i\in{\mathbb{N}}}$. Since each $F_{i}$ is a measurable
function of ${\bm{G}}_{i},\dots,{\bm{G}}_{i+r}$ it follows from $(d.5)$ and
Remark 1.8 of Bradley (2007) that $(\lambda^{*}F_{i})_{i\in{\mathbb{N}}}$ is
$\alpha$-mixing with mixing coefficients satisfying
$\alpha_{F}(n)\leq\alpha_{{\bm{G}}}(n-(r+1))$ for all $n\geq r+2$. Therefore,
$\sum_{n=0}^{\infty}(\alpha_{F}(n))^{\frac{\epsilon}{2+\epsilon}}<\infty$ for
all $\epsilon>0$. From $(e.1)$ we obtain
${\mathbb{E}}\|\lambda^{\ast}F_{1}\|^{2+\epsilon/4}<\infty$ for some
$\epsilon>0$. Thus, the CLT for $\alpha$-mixing sequences applies, see e.g.
(Ibragimov and Linnik, 1971, Theorem 18.5.3), so that

$\sqrt{n}\Big{(}\frac{1}{n}\sum_{i=1}^{n-r}\lambda^{\ast}F_{i}\Big{)}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\zeta),\quad{n\to\infty},$

where

$\zeta={\mathbb{E}}\lambda^{\ast}F_{1}F_{1}^{\ast}\lambda+2\sum_{i=1}^{\infty}{\mathbb{E}}\lambda^{\ast}F_{1}F_{1+i}^{\ast}\lambda.$

Since
$\lambda^{\ast}F_{1}F_{1+i}^{\ast}\lambda=\lambda^{\ast}(F_{1}F_{1+i}^{\ast})^{\ast}\lambda$
we get (5.6) after rearranging the above equation. ∎

######  Theorem 5.5.

Assume that Assumptions a,b, c, d and $(e.1)$ hold and that the matrix
$\Sigma$ in (5.6) is positive definite. Then the GMM estimator defined in
(5.3) is asymptotically normal with covariance matrix

$(\mathcal{J}_{\theta_{0}})^{-1}\mathcal{I}_{\theta_{0}}(\mathcal{J}_{\theta_{0}})^{-1},$
(5.7)

where
$\mathcal{J}_{\theta_{0}}=(\nabla_{\theta}k_{\theta_{0},r})^{\top}\Omega(\nabla_{\theta}k_{\theta_{0},r})$
and
$\mathcal{I}_{\theta_{0}}=(\nabla_{\theta}k_{\theta_{0},r})^{\top}\Omega\Sigma_{\theta_{0}}\Omega(\nabla_{\theta}k_{\theta_{0},r})$.

###### Proof.

We check Assumptions 1.7-1.9 of Theorem 1.2 in Mátyás (1999). Since by Lemma
5.3 the map $\theta\mapsto k_{\theta,r}$ is continuously differentiable, their
Assumption 1.7 is valid. Now, for any sequence $\tilde{\theta}_{n}$ such that
$\tilde{\theta}_{n}\stackrel{{\scriptstyle P}}{{\rightarrow}}\theta_{0}$ as
${n\to\infty}$, it follows from the continuous mapping theorem by the
continuity of the map
$\Theta\mapsto\frac{\partial}{\partial\theta}k_{\theta,r}$ in Lemma 5.3 that
$\frac{\partial}{\partial\theta}(\hat{k}_{n,r}-k_{\theta_{n}})\stackrel{{\scriptstyle
P}}{{\rightarrow}}(k_{\theta_{0}}-\frac{\partial}{\partial\theta}k_{\theta_{0}})$
as ${n\to\infty}$. Therefore, Assumption 1.8 in Mátyás (1999) is also
satisfied. Since Lemma 5.4 implies Assumption 1.9, we conclude the result. ∎

######  Remark 5.6.

In order to apply the results of Section 5.2 we need to check Assumption c,
model identifiability $(d.3)$, strong mixing of the log-price returns sequence
$(d.5)$ and existence of certain moments of its stationary distribution
(Assumption e) . In Sections 5.3 and 5.4 we give sufficient conditions for
identifiability of the model parameters, strict stationarity and strong
mixing. Then we use these results to derive in Section 5.5 more palpable
conditions under which Theorems 5.2 and 5.5 can be applied.

###  5.3 Sufficient conditions for strict stationarity and strong mixing

Sufficient conditions for the existence of a unique stationary distribution of
$(Y_{t})_{t\in{\mathbb{R}}^{+}}$, geometric ergodicity and for the finiteness
of moments of order $p$ of the stationary distribution have recently been
given in Stelzer and Vestweber (2019). We state these conditions in the next
theorem, which are conditions (i), (iv) and (v) of Theorem 4.3 in Stelzer and
Vestweber (2019).

######  Theorem 5.7 (Geometric Ergodicity - (Stelzer and Vestweber, 2019,
Theorem 4.3)).

Let $Y$ be a MUCOGARCH volatility process which is $\mu$-irreducible with the
support of $\mu$ having non-empty interior and aperiodic. Assume that one of
the following conditions is satisfied:

  * (i)

setting $p=1$ there exists $\Xi\in\mathbb{S}_{d}^{++}$ such that

$\Xi B+B^{\top}\Xi+\sigma_{L}A^{\top}\Xi A\in-\mathbb{S}_{d}^{++},$ (5.8)

  * (ii)

there exist $p\in[1,\infty)$ and $\Xi\in\mathbb{S}_{d}^{++}$ such that

$\int_{\mathbb{R}^{d}}\left(2^{p-1}\left(1+K_{\Xi,A}\|y\|_{2}^{2}\right)^{p}-1\right)\nu_{L}(dy)+pK_{\Xi,B}<0,$
(5.9)

where

$K_{\Xi,B}=\max_{X\in\mathbb{S}_{d}^{+},\operatorname{tr}(X)=1}\frac{\operatorname{tr}\left(\left(\Xi
B+B^{\top}\Xi\right)X\right)}{\operatorname{tr}(\Xi X)}\text{ and
}K_{\Xi,A}=\max_{X\in\mathbb{S}_{d}^{+},\operatorname{tr}(X)=1}\frac{\operatorname{tr}\left(A^{\top}\Xi
AX\right)}{\operatorname{tr}(\Xi X)},$

  * (iii)

there exist $p\in[1,\infty)$ and $\Xi\in\mathbb{S}_{d}^{++}$ such that

$\max\left\\{2^{p-2},1\right\\}K_{\Xi,A}\int_{\mathbb{R}^{d}}\|y\|_{2}^{2}\left(1+\|y\|_{2}^{2}K_{\Xi,A}\right)^{p-1}\nu_{L}(dy)+K_{\Xi,B}<0$
(5.10)

where $K_{\Xi,B},K_{\Xi,A}$ are as in (ii).

Then a unique stationary distribution for the MUCOGARCH volatility process $Y$
exists, $Y$ is positive Harris recurrent, geometrically ergodic and its
stationary distribution has a finite $p$-th moment.

A consequence of Theorem 5.7 is that the process $Y$ is exponentially
$\beta$-mixing. This implies $\alpha$-mixing of the log-price process as we
state next. For more details on mixing conditions we refer to Bradley (2007).

######  Corollary 5.8.

If $Y$ is strictly stationary and exponentially $\beta$-mixing, then the log-
price process $({\bm{G}}_{i})_{i\in{\mathbb{N}}}$ is stationary, exponentially
$\alpha$-mixing, and as a consequence also ergodic.

###### Proof.

Since $Y$ is an exponentially $\beta$-mixing, homogeneous strong Markov
process (Stelzer, 2010, Theorem 4.4), and driven only by the discrete part of
the quadratic variation of $L$, the proof follows by the same arguments as for
Theorem 3.4 in Haug et al. (2007). ∎

Next, we state a result which gives sufficient conditions for the
irreducibility of the MUCOGARCH volatility process $Y$ process, which is one
of the sufficient conditions for the geometric ergodicity result in Theorem
5.7.

######  Theorem 5.9 (Irreducibility and Aperiodicity - (Stelzer, 2010, Theorem
5.1 and Corollary 5.2)).

Let $Y$ be a MUCOGARCH volatility process driven by a Lévy process whose
discrete part is a compound Poisson process $L$ with $A\in GL_{d}(\mathbb{R})$
and $\Re(\sigma(B))<0$. If the jump distribution of L has a non-trivial
absolutely continuous component equivalent to the Lebesgue measure on
${\mathbb{R}}^{d}$ restricted to an open neighborhood of zero, then $Y$ is
irreducible w.r.t. the Lebesgue measure restricted to an open neighborhood of
zero in $\mathbb{S}_{d}^{+}$ and aperiodic.

###  5.4 Sufficient conditions for identifiability

In this Section we investigate the identifiability of the model parameters
from the model moments, i.e., we investigate the injectivity of the map
$\theta\mapsto k_{\theta,r}$ on an appropriate compact set $\Theta$. Recall
that we can divide the this map into the composition of
$\theta\mapsto(A_{\theta},B_{\theta},C_{\theta})\mapsto k_{\theta,r}$.
Injectivity of $\theta\mapsto(A_{\theta},B_{\theta},C_{\theta})$ holds if e.g.
it simply maps the entries of $\theta$ to the entries of the matrices
$(A_{\theta},B_{\theta},C_{\theta})$. Thus, we only need to investigate the
injectivity of the map $(A_{\theta},B_{\theta},C_{\theta})\mapsto
k_{\theta,r}$. As we will see, there will appear some restrictions on the
matrices $A_{\theta},B_{\theta}$, which are related to the fact that we need
to take the logarithm of a matrix exponential, and we need to ensure this is
well defined. We will omit $\theta$ from the notation, except when explicitly
needed. We start with the identifiability of the matrix $C$.

######  Lemma 5.10.

Assume that Assumptions $(a.1)$-$(a.5)$, $(b.2)$ and $(c.1)$ hold and that
$\sigma(B)\subset\\{z\in{\mathbb{C}}:\Re(z)<0\\}$. If the matrices $A$ and $B$
are known, then ${\mathbb{E}}({\bm{G}}_{1}{\bm{G}}_{1}^{\ast})$ uniquely
determines $C$.

###### Proof.

Since $\sigma(B\otimes I+I\otimes
B)=\sigma(B)+\sigma(B)\subset\\{z\in{\mathbb{C}}:\Re(z)<0\\}$, the matrix
$B\otimes I+I\otimes B$ is invertible. The rest of the proof follows by noting
that from (4.2) and (4.3) it follows that

$\operatorname{vec}(C)=(\sigma_{L}+\sigma_{W})^{-1}\Delta^{-1}(B\otimes
I+I\otimes
B)^{-1}{\mathcal{B}}\operatorname{vec}({\mathbb{E}}({\bm{G}}_{1}{\bm{G}}_{1}^{\ast})).$

∎

For the identification of the matrices $A$ and $B$ we need to use the second-
order structure of the squared returns process in Lemma 4.2. We first state
three auxiliary results, which provide conditions such that we can identify
the components of the autocovariance function in (4.4).

######  Lemma 5.11.

Assume that $B\in M_{d}({\mathbb{R}})$ is diagonalizable with $S\in
GL_{d}({\mathbb{C}})$ such that $S^{-1}BS$ is diagonal. If

$\bigg{|}\frac{\sigma_{L}-\sigma_{W}}{2\sigma_{L}}\bigg{|}\|A\otimes
A\|_{S}<-2\max\\{\Re(\sigma(B))\\},$ (5.11)

with

$\|X\|_{S}=\|(S^{-1}\otimes S^{-1})X(S\otimes S))\|_{2},\quad X\in
M_{d^{2}}({\mathbb{R}}),S\in GL_{d}({\mathbb{C}}),$ (5.12)

then the matrix

$(\sigma_{W}+\sigma_{L})({\mathcal{B}}^{\ast})^{-1}-2((A\otimes
A)^{\ast})^{-1}$ (5.13)

is invertible.

###### Proof.

From (Bernstein, 2009, fact 2.16.14), $X^{-1}+Y^{-1}$ is non-singular if and
only if $X+Y$ is non-singular and $X,Y$ are non-singular. Setting
$X=\frac{{\mathcal{B}}}{(\sigma_{L}+\sigma_{W})},Y=-\frac{1}{2}(A\otimes A)$
and using the definition of ${\mathcal{B}}$ in (4.1) we get

$X+Y=\frac{1}{(\sigma_{L}+\sigma_{W})}\bigg{(}(B\otimes I+I\otimes
B)+\frac{(\sigma_{L}-\sigma_{W})}{2}(A\otimes A)\bigg{)}.$

Since $B$ is diagonalizable, we can use (Bernstein, 2009, Proposition 7.1.6)
to obtain

$B\otimes I+I\otimes B=(S\otimes S)(S^{-1}BS\otimes I)(S^{-1}\otimes S^{-1}),$

which guarantees that $B\otimes I+I\otimes B$ is also diagonalizable. Now we
rewrite the first equation on p. 106 in Stelzer (2010) with the matrix
${\mathcal{B}}$ replaced by $(B\otimes I+I\otimes
B)+\frac{(\sigma_{L}-\sigma_{W})}{2}(A\otimes A)$ and apply the Bauer-Fike
Theorem (Horn and Johnson, 1991, Theorem 6.3.2) to see that (5.11) implies
that all eigenvalues of $(X+Y)(\sigma_{L}+\sigma_{W})$ are in
$\\{z\in{\mathbb{C}}:\Re(z)<0\\}$ and, therefore, $X+Y$ is invertible. ∎

######  Lemma 5.12.

If $A\in M_{d}({\mathbb{R}})$ is such that $A_{(1,1)},\dots,A_{(1,{j-1})}=0$
and $A_{(1,j)}>0$ for some $j\in\\{1,\dots,d\\}$, then the map $X\mapsto
AXA^{T}$ for $X\in{\mathbb{S}}_{d}$ identifies $A$.

###### Proof.

Assume first that $A_{(1,1)}>0$. For each $i\in\\{1,\dots,d\\}$, let $e_{i}$
be the $i$th column unit vector in ${\mathbb{R}}^{d}$ and define the matrix
$E^{(i,j)}=e_{i}e_{j}^{T}$. The first line of the matrix $AE^{(1,1)}A^{T}$ is

$(A_{(1,1)}^{2},A_{(1,1)}A_{(2,1)},\dots,A_{(1,1)}A_{(d,1)}).$ (5.14)

Since $A_{(1,1)}>0$, (5.14) allows us to identify first $A_{(1,1)}$ and then
$A_{(2,1)},\dots,A_{(d,1)}$. Now, for each $k\in\\{2,\dots,d\\}$, note that
$E^{(1,k)}+E^{(k,1)}$ is symmetric. Simple calculations reveal that the first
line of the matrix $A(E^{(1,k)}+E^{(k,1)})A^{T}$ is

$(2A_{(1,1)}A_{(1,k)},A_{(1,1)}A_{(2,k)}+A_{(1,k)}A_{(2,1)},\dots,A_{(1,1)}A_{(d,k)}+A_{(1,k)}A_{(d,1)}).$
(5.15)

Since $A_{(1,1)}>0$, we identify $A_{(1,k)}$ from the first entry of (5.15).
Now, since also $A_{(2,1)},\dots,A_{(d,1)}$ are already known, we can identify
$A_{(2,k)},\dots,A_{(d,k)}$. Thus, all entries of A can be identified. The
cases $A_{(1,j)}>0$ for some $j>1$ follow similarly. ∎

######  Lemma 5.13.

Assume that the Assumptions a,b and c and the conditions of Lemma 5.12 hold,
that the matrix in (5.13) is invertible, that
$\sigma({\mathcal{B}})\subset\\{z\in{\mathbb{C}}:-\pi<\Im(z)\Delta<\pi,\Re(z)<0\\}$
and that ${\rm var}(\operatorname{vech}(V_{0}))$ is invertible. Define
$M=(e^{{\mathcal{B}}\Delta})^{-1}{\rm acov}_{\theta,\bm{G}\bm{G}^{\ast}}(1)$.
Then, ${\rm acov}_{\theta,\bm{G}\bm{G}^{\ast}}(1)$ and ${\rm
acov}_{\theta,\bm{G}\bm{G}^{\ast}}(2)$ uniquely identify ${\mathcal{B}}$ and
$M$.

###### Proof.

Since $M$ is given in terms of ${\mathcal{B}}$ and ${\rm
acov}_{\theta,\bm{G}\bm{G}^{\ast}}(1)$, we only need to identify
${\mathcal{B}}$. Observe that we are using the $\operatorname{vec}$ operator
only for convenience, as it interacts nicely with tensor products of matrices
and thus gives nicely looking formulae. However, the volatility and “squared
returns” processes take values in $\mathbb{S}_{d}$ which is a
$d(d+1)/2$-dimensional vector space, whereas the $\operatorname{vec}$ operator
assumes values in a $d^{2}$-dimensional vector space. Instead of using the
$\operatorname{vech}$ operator and cumbersome notation, we take an abstract
point of view. The variance of a random element of $\mathbb{S}_{d}$ is a
symmetric positive semi-definite linear operator from $\mathbb{S}_{d}$ to
itself. Likewise, the autocovariance of $\bm{G}_{1}\bm{G}_{1}^{\ast}$ and
$\bm{G}_{1+h}\bm{G}_{1+h}^{\ast}$ is a linear operator from $\mathbb{S}_{d}$
to itself. The condition that ${\rm var}(\operatorname{vech}(V_{0}))$ is
invertible is equivalent to the invertibility of the linear operator, which is
the variance of $V_{0}$. Similarly all other $d^{2}\times d^{2}$ matrices in

$e^{{\mathcal{B}}\Delta
h}{\mathcal{B}}^{-1}(I_{d^{2}}-e^{-{\mathcal{B}}\Delta})(\sigma_{L}+\sigma_{W}){\rm
var}(\operatorname{vec}(V_{0}))(e^{{\mathcal{B}}^{\ast}\Delta}-I_{d^{2}})[(\sigma_{W}+\sigma_{L})({\mathcal{B}}^{\ast})^{-1}-2((A\otimes
A)^{\ast})^{-1}]$

are representing linear operators from $\mathbb{S}_{d}$ to itself. Under the
assumptions made, the above product involves only invertible linear operators.
Hence ${\rm acov}_{\theta,\bm{G}\bm{G}^{\ast}}(h)$ is invertible (over
$\mathbb{S}_{d}$) for every $h>0$. Thus,

$e^{{\mathcal{B}}\Delta}={\rm acov}_{\theta,\bm{G}\bm{G}^{\ast}}(2)[{\rm
acov}_{\theta,\bm{G}\bm{G}^{\ast}}(1)]^{-1}.$

By the assumptions on the eigenvalues of ${\mathcal{B}}$ there is a unique
logarithm for $e^{{\mathcal{B}}\Delta}$ (see (Horn and Johnson, 1991, Section
6.4) or (Schlemm and Stelzer, 2012, Lemma 3.11)), so ${\mathcal{B}}\Delta$ and
thus ${\mathcal{B}}$ is identified. Finally, note that the matrices in the
$\operatorname{vec}$ representations are uniquely identified by the employed
linear operators on $\mathbb{S}_{d}$ due to (Pigorsch and Stelzer, 2009a,
Proposition 3.1) and Lemma 5.12. ∎

######  Lemma 5.14 (Identifiability of $A$, $B$ and $C$).

For all $\theta\in\Theta$, assume the conditions of Lemma 5.13,
$\sigma(B_{\theta})\subset\\{z\in{\mathbb{C}}:\Re(z)<0\\}$ and that the
entries of the matrices $A_{\theta}$ and $B_{\theta}$ satisfy: for some
$k\not=l\in\\{1,\dots,d\\}$, $A_{(k,l),\theta}>0$,
$A_{(k,l),\theta}\not=A_{(l,k),\theta}$ and
$B_{(k,l),\theta}=B_{(l,k),\theta}$. Then $k_{\theta,2}$ uniquely identifies
$A_{\theta}$, $B_{\theta}$ and $C_{\theta}$.

###### Proof.

Recall that we omit $\theta$ in the notation. Assume w.l.o.g that
$\sigma_{L}=1$. Because of Lemma 5.10, we only need to show the identification
of $A$ and $B$.
Assume first that $d=2$. Then the matrix ${\mathcal{B}}$ from (4.1) is

$\begin{pmatrix}[r]2B_{(1,1)}+A_{(1,1)}^{2}&B_{(1,2)}+A_{(1,1)}A_{(1,2)}&B_{(1,2)}+A_{(1,1)}A_{(1,2)}&A_{(1,2)}^{2}\\\
B_{(2,1)}+A_{(1,1)}A_{(2,1)}&B_{(1,1)}+B_{(2,2)}+A_{(1,1)}A_{(2,2)}&A_{(1,2)}A_{(2,1)}&B_{(1,2)}+A_{(1,2)}A_{(2,2)}\\\
B_{(2,1)}+A_{(1,1)}A_{(2,1)}&A_{(1,2)}A_{(2,1)}&B_{(1,1)}+B_{(2,2)}+A_{(1,1)}A_{(2,2)}&B_{(1,2)}+A_{(1,2)}A_{(2,2)}\\\
A_{(2,1)}^{2}&B_{(2,1)}+A_{(2,1)}A_{(2,2)}&B_{(2,1)}+A_{(2,1)}A_{(2,2)}&2B_{(2,2)}+A_{(2,2)}^{2}\end{pmatrix}.$
(5.16)

Using the entry at position $(1,4)$ and the fact that $A_{(1,2)}>0$ allow us
to identify $A_{(1,2)}$. Then, we use the entry at position $(2,3)$ to
identify $A_{(2,1)}$. Now, we use the entries at positions $(1,2)$ and $(2,1)$
together with the fact that $A_{(1,2)}\not=A_{(2,1)}$ and
$B_{(1,2)}=B_{(2,1)}$ to write
$A_{(1,1)}=({\mathcal{B}}_{(1,2)}-{\mathcal{B}}_{(2,1)})/(A_{(1,2)}-A_{(2,1)})$.
Similarly we use the entries at positions $(3,4),(4,3)$ to get
$A_{(2,2)}=({\mathcal{B}}_{(3,4)}-{\mathcal{B}}_{(4,3)})/(A_{(1,2)}-A_{(2,1)})$.
Now, since all the entries of $A$ are known, we can use the entries at
positions $(1,1),(1,2)$ and $(2,2)$ to identify the entries of $B$.
Now assume that $d>2$. We assume w.l.o.g. that $A_{(1,2)}>0$,
$A_{(1,2)}\not=A_{(2,1)}$ and $B_{(1,2)}=A_{(2,1)}$. Write the matrix
${\mathcal{B}}$ from (4.1) in the following block form:

${\mathcal{B}}=B\otimes I+I\otimes B+A\otimes
A=\begin{pmatrix}{\mathcal{B}}^{(1,1)}&\cdots&{\mathcal{B}}^{(1,d)}\\\
\vdots&\ddots&\vdots\\\
{\mathcal{B}}^{(d,1)}&\cdots&{\mathcal{B}}^{(d,d)}\end{pmatrix},$ (5.17)

where ${\mathcal{B}}^{(i,j)}\in M_{d}({\mathbb{R}})$ for all $i,j=1,\dots,d$.
First, we have that

${\mathcal{B}}^{(1,2)}=\begin{pmatrix}[r]B_{(1,2)}+A_{(1,2)}A_{(1,1)}&A_{(1,2)}A_{(1,2)}&A_{(1,2)}A_{(1,3)}&\cdots&A_{(1,2)}A_{(1,d)}\\\
A_{(1,2)}A_{(2,1)}&B_{(1,2)}+A_{(1,2)}A_{(2,2)}&A_{(1,2)}A_{(2,3)}&\cdots&A_{(1,2)}A_{(2,d)}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
A_{(1,2)}A_{(d,1)}&A_{(1,2)}A_{(d,2)}&A_{(1,2)}A_{(d,3)}&\cdots&B_{(1,2)}+A_{(1,2)}A_{(d,d)}\end{pmatrix},$
(5.18)

Since $A_{(1,2)}>0$ we can identify it from (5.18), because
${\mathcal{B}}^{(1,2)}_{(1,2)}=A_{(1,2)}^{2}$. Then we use the off-diagonal
entries of the matrix ${\mathcal{B}}^{(1,2)}$ in (5.18) together with
$A_{(1,2)}$ to identify all the off-diagonal entries of the matrix $A$. Next
we identify the diagonal entries of $A$. It follows from (5.17) that

$\systeme*{{\mathcal{B}}^{(k,k)}_{(1,2)}=B_{(1,2)}+A_{(k,k)}A_{(1,2)},{\mathcal{B}}^{(k,k)}_{(2,1)}=B_{(2,1)}+A_{(k,k)}A_{(2,1)}},\quad
k=1,\dots,d.$ (5.19)

Since $A_{(1,2)}-A_{(2,1)}\not=0$ and $B_{(1,2)}=B_{(2,1)}$, the system of
equations (5.19) gives

$A_{(k,k)}=({\mathcal{B}}^{(k,k)}_{(1,2)}-{\mathcal{B}}^{(k,k)}_{(2,1)})/(A_{(1,2)}-A_{(2,1)}),\quad
k=1,\dots,d.$

Finally, since the matrix $A$ is now completely known, we can use (5.17) to
identify all entries of $B$. ∎

In Lemma 5.14 we identify the matrices $A$ and $B$ only from ${\mathcal{B}}$
and, therefore, some mild restrictions on the off-diagonal entries of $B$
appear. In order to avoid those restrictions, we could to take the structure
of
${\mathbb{E}}\operatorname{vec}(\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})^{*})$
in (4.5) into account when proving identifiability and we expect that one can
improve the identification results since more moment conditions are used.
However, already in the $2$-dimensional case the results on identification
conditions are quite involved, and this has mainly to do with the fact that
the linear operator $({\bm{Q}}+K_{d}{\bm{Q}}+I_{d^{2}})$ at the right hand
side of (4.5) is not one-to-one in the space of matrices of the form
${\mathbb{E}}\operatorname{vec}(V_{0})\operatorname{vec}(V_{0})^{\ast}$. In
the end, in order to use the moment conditions
${\mathbb{E}}\operatorname{vec}(\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})^{*})$,
we need to assume that the matrices $B\otimes I+I\otimes B$ and $A\otimes A$
commute (see (Do Rêgo Sousa, 2019, Lemma 3.5.18)). Since commutativity is a
quite strong condition, it seems highly preferable to work with the class of
MUCOGARCH processes, which are identifiable by Lemma 5.14. The exponential
decay of the autocovariance function of the model is still quite flexible,
because of the interplay between the matrices $A$ and $B$ (see (5.16), for
instance).

###  5.5 Asymptotic properties: general case revisited

Here, we combine the results of Sections 5.2-5.4 to give easily verifiable
conditions under which the GMM estimator $\hat{\theta}_{n}$ will be consistent
and asymptotically normal. We assume that the parameter $\theta$ contains the
entries of the matrices $(A_{\theta},B_{\theta},C_{\theta})$ so that the map
$\theta\mapsto(A_{\theta},B_{\theta},C_{\theta})$ is automatically injective
and continuously differentiable on $\Theta$.

First, we define

$\displaystyle\|x\|_{S}$ $\displaystyle=$ $\displaystyle\|(S^{-1}\otimes
S^{-1})x\|_{2},\quad x\in{\mathbb{R}}^{d^{2}},S\in GL_{d}({\mathbb{C}})$
(5.20) $\displaystyle K_{2,S}$ $\displaystyle=$
$\displaystyle\max_{X\in\mathbb{S}^{+}_{d},\|X\|_{2}=1}\bigg{(}\frac{\|X\|_{2}}{\|\operatorname{vec}(X)\|_{S}}\bigg{)},\quad
S\in GL_{d}({\mathbb{C}}).$ (5.21)

Consider now the following group of assumptions:

######  Assumptions f (Parameter space).

For all $\theta\in\Theta$ it holds:

  1. $(f.1)$

The matrices $B_{\theta}$ satisfy
$\sigma(B_{\theta})\subset\\{z\in{\mathbb{C}}:\Re(z)<0\\}$.

  2. $(f.2)$

The matrix ${\mathcal{B}}_{\theta}$ satisfy
$\sigma({\mathcal{B}}_{\theta})\subset\\{z\in{\mathbb{C}}:-\pi<\Im(z)\Delta<\pi,\Re(z)<0\\}$.

  3. $(f.3)$

The matrix $B_{\theta}\in M_{d}({\mathbb{R}})$ is diagonalizable with
$S_{\theta}\in GL_{d}({\mathbb{C}})$ such that
$S_{\theta}^{-1}B_{\theta}S_{\theta}$ is diagonal.

  4. $(f.4)$

The entries of the matrices $A_{\theta}$ and $B_{\theta}$ satisfy: for some
$k\not=l\in\\{1,\dots,d\\}$, $A_{(k,l),\theta}>0$,
$A_{(k,l),\theta}\not=A_{(l,k),\theta}$ and
$B_{(k,l),\theta}=B_{(l,k),\theta}$.

  5. $(f.5)$

The matrix ${\rm var}_{\theta}(\operatorname{vech}(V_{0}))$ is invertible.

  6. $(f.6)$

$|\frac{\sigma_{L}-\sigma_{W}}{2\sigma_{L}}|\|A_{\theta}\otimes
A_{\theta}\|_{S_{\theta}}<-2\max\\{\Re(\sigma(B_{\theta}))\\}$ with
$S_{\theta}$ as in $(f.3)$ and $\|A_{\theta}\otimes A_{\theta}\|_{S_{\theta}}$
as in (5.12).

  7. $(f.7)$

There exists $\Xi_{\theta}\in\mathbb{S}_{d}^{++}$ such that, condition (5.8)
holds with $A,B$ replaced by $A_{\theta},B_{\theta}$.

  8. $(f.8)$

$m(4,\theta)<0$ where

$m(p,\theta):=\int_{{\mathbb{R}}^{d}}((1+\alpha_{\theta}\|\operatorname{vec}(yy^{\ast})\|_{S_{\theta}})^{p}-1)\nu_{L}({\rm
d}y)+2p\max\\{\mathfrak{R}(\sigma(B_{\theta}))\\},$ (5.22)

$\alpha_{\theta}=\|S_{\theta}\|_{2}^{2}\|S_{\theta}^{-1}\|_{2}^{2}K_{2,B_{\theta}}\|A_{\theta}\otimes
A_{\theta}\|_{S_{\theta}}$ with $K_{2,B_{\theta}}$ as in (5.21),
$\|\operatorname{vec}(yy^{\ast})\|_{S_{\theta}}$ as in (5.20) and $S_{\theta}$
as in $(f.3)$.

######  Assumptions g (MUCOGARCH process at $\theta_{0}$).

  1. $(g.1)$

The MUCOGARCH volatility process $Y$ is stationary, $\mu$-irreducible with the
support of $\mu$ having non-empty interior and aperiodic.

  2. $(g.2)$

$m(p,\theta_{0})<0$ for some $p>4$.

Assumption $(f.1)$-$(f.6)$ collect the needed identifiability assumptions from
Section 5.4. Assumption $(f.7)$ is a sufficient condition under which we have
uniqueness of the stationary distribution of $Y$ and geometric ergodicity
(Section 5.3). For the asymptotic results of the GMM estimator in Theorems 5.2
and 5.5, we need to ensure ${\mathbb{E}}\|Y_{0}\|^{p}<\infty$ for appropriate
$p>1$, and this would require checking Assumptions 5.9 or 5.10 with $p>1$.
However, this imposes strong conditions on the Lévy process (Stelzer and
Vestweber, 2019, Remark 4.4). Instead, we require diagonalizability of the
matrix $B_{\theta}$ ($(f.3)$), which is not a very restrictive assumption, and
check (5.22) to ensure ${\mathbb{E}}\|Y_{0}\|^{p}<\infty$ (Stelzer, 2010,
Theorem 4.5), which is less restrictive.

In view of the above assumptions and the results of Sections 5.1-5.4 we have
the following consistency result.

######  Corollary 5.15 (Consistency of the GMM estimator - $L$ has paths of
infinite variation).

Suppose that assumptions a, b, $(d.1)$,$(f.1)$-$(f.8)$ and $(g.1)$ hold. Then
the GMM estimator defined in (5.3) is weakly consistent.

If the paths of the driving Lévy process are of finite variation, we can relax
even more the conditions from Corollary 5.15. Before we state this result, we
give the definition of asymptotic second-order stationarity which will be used
in its proof. A stochastic process $X\in{\mathbb{S}}_{d}$ is said to be
asymptotically second-order stationary with mean $\mu\in{\mathbb{R}}^{d^{2}}$,
variance $\Sigma\in{\mathbb{S}}_{d^{2}}^{+}$ and autocovariance function
$f:{\mathbb{R}}^{+}\mapsto M_{d^{2}}({\mathbb{R}})$ if it has finite second
moments and

$\begin{split}&\lim_{t\rightarrow\infty}{\mathbb{E}}\left(X_{t}\right)=\mu,\quad\quad\lim_{t\rightarrow\infty}{\rm
var}\left(\operatorname{vec}\left(X_{t}\right)\right)=\Sigma\\\
&\lim_{t\rightarrow\infty}\sup_{h\in\mathbb{R}^{+}}\left\\{\left\|{\rm
cov}\left(\operatorname{vec}\left(X_{t+h}\right),\operatorname{vec}\left(X_{t}\right)\right)-f(h)\right\|\right\\}=0.\end{split}$

######  Corollary 5.16 (Consistency of the GMM estimator - $L$ has paths of
finite variation).

Suppose that assumptions a, b, $(d.1)$, $(f.1)$, $(f.2)$, $(f.4)$-$(f.7)$,
$(g.1)$ hold and that $L$ has paths of finite variation. Then, the GMM
estimator defined in (5.3) is weakly consistent.

###### Proof.

Let $D\in\mathbb{S}_{d}^{+}$ be a constant matrix, and consider a MUCOGARCH
process $(Y_{t})_{t\in{\mathbb{R}}^{+}}$ solving (3.3) have starting value
$D$. Then, a combination of Assumptions $(a.2)$, $(a.5)$ with the fact that
the starting value $D$ is non-random and the hypothesis imposed on the
matrices $B_{\theta},{\mathcal{B}}_{\theta},{\mathcal{C}}_{\theta}$ allow us
to apply Theorem 4.20(ii) in Stelzer (2010) to conclude that the process
$(Y_{t})_{t\in{\mathbb{R}}^{+}}$ is asymptotically second-order stationary.
Additionally, Theorem 5.7(i) ensures that the process
$(Y_{t})_{t\in{\mathbb{R}}^{+}}$ has a unique stationary distribution, is
geometrically ergodic and its stationary distribution has finite first moment,
i.e., ${\mathbb{E}}\|Y_{0}\|<\infty$. Since $Y_{t}\in{\mathbb{S}}_{d}^{+}$,
and $\text{tr}(Y^{*}Y)$ (with tr denoting the usual trace functional) defines
a scalar product on ${\mathbb{S}}_{d}$ via
$\text{tr}(Y_{t}^{*}Y_{t})=\operatorname{vec}(Y_{t}^{*})\operatorname{vec}(Y_{t})$
it follows that

$\begin{split}{\mathbb{E}}\|Y_{t}\|_{2}^{2}=&\text{tr}(Y_{t}^{*}Y_{t})=\operatorname{vec}(Y_{t}^{*})\operatorname{vec}(Y_{t})=\sum_{i,j}^{d}{\mathbb{E}}Y_{t,ij}^{2}=\sum_{i,j}^{d}{\rm
var}(Y_{t,ij})+\sum_{i,j}^{d}({\mathbb{E}}Y_{t,ij})^{2}\\\ &\text{tr}({\rm
var}(Y_{t}))+\|{\mathbb{E}}(Y_{t})\|_{2}^{2},\quad t>0.\end{split}$ (5.23)

Since both maps $t\mapsto{\mathbb{E}}\|Y_{t}\|$ and $t\mapsto{\rm var}(Y_{t})$
are continuous ((Stelzer, 2010, eqs. (4.7) and (4.16))), it follows from
(5.23) that $\limsup_{t\geq 0}{\mathbb{E}}\|Y_{t}\|^{2}<\infty$. Since Theorem
5.7(i) implies convergence of the transition probabilities in total variation,
which in turns implies weak convergence (e.g. (Klenke, 2013, Exercise
13.2.2)), we have that $Y_{t}\stackrel{{\scriptstyle d}}{{\rightarrow}}Y_{0}$
as ${t\to\infty}$, with $Y_{0}$ being the stationary version of $Y$. Hence, we
can use the continuous mapping theorem and (Billingsley, 2008, Theorem 25.11)
to conclude that ${\mathbb{E}}\|Y_{0}\|^{2}<\infty$. Finally, the result
follows by an application of Lemma 3.2, Theorem 5.2, Corollary 5.8 and Remark
A.2. ∎

Recall that for the asymptotic normality result, we need to ensure that the
stationarity distribution of the MUCOGARCH volatility process has more than
$4$ moments (cf. $(e.1)$). This is summarized in the next corollary.

######  Corollary 5.17 (Asymptotic normality of the GMM estimator).

If Assumptions a, b, $(d.1)$, $(d.2)$, $(d.4)$, f and g hold, then the GMM
estimator defined in (5.3) is asymptotically normal with covariance matrix as
in (5.7).

###### Proof.

By the same arguments of (Stelzer, 2010, Theorem 4.5) combined with (Lindner
and Maller, 2005, Proposition 4.1) and $(g.2)$, it follows that
${\mathbb{E}}\|Y_{0}\|^{p}<\infty$ for some $p>4$. The rest of the proof is
just an application of Theorem 5.5. ∎

######  Remark 5.18.

The advantage of Corollaries 5.15-5.17 is that Assumption f can be checked
numerically and Assumption g holds true if e.g., the Lévy process $L$ is a
compound Poisson process with jump distribution having a density which is
strictly positive around zero (see Theorem 5.9).

If $(g.1)$ holds, the stationary distribution of $Y$ is automatically a
maximum irreducibility measure. All maximal irreducibility measures are
equivalent and thus the support of the stationary distribution has a support
which has non-empty interior. The latter in turn implies that the variance has
to be an invertible operator (non-invertibility is equivalent to the
distribution being concentrated on a proper linear subspace) which is $(f.5)$
for $\theta_{0}$.

In the next section, we investigate the finite sample performance of the
estimators in a simulation study.

##  6 Simulation study

To assess the performance of the GMM estimator, we will focus on the MUCOGARCH
model in dimension $d=2$. We fix
$L_{t}=L^{\mathfrak{d}}_{t}+\sqrt{\sigma_{W}}W_{t}$ for $t\in{\mathbb{R}}^{+}$
where $L^{\mathfrak{d}}$ is a bivariate compound Poisson process (CPP), $W$ is
a standard bivariate Brownian motion, independent of $L^{\mathfrak{d}}$ and
$\sigma_{W}\geq 0$ is fixed. We choose $L^{\mathfrak{d}}$ as a CPP, since it
allows to simulate the MUCOGARCH volatility process $V$ exactly. Thus, we only
need to approximate the Brownian part of the (log) price process $G$ in (1.1),
which is done by an Euler scheme. Setting $L^{\mathfrak{d}}$ as a CPP is not a
very crucial restriction, since for an infinite activity Lévy process one
would need to approximate it using only finitely many jumps. For example by
using a CPP for the big jumps component of $L^{\mathfrak{d}}$ and an
appropriate Brownian motion for its small jumps component (see Cohen and
Rosinski (2007)). In applications, a CPP has also been used in combination
with the univariate COGARCH(1,1) process for modeling high frequency data (see
Müller (2010)). The jump distribution of $L^{\mathfrak{d}}$ is chosen as
$N(0,1/4I_{2})$ and the jump rate is $4$, so that ${\rm var}(L_{1})=2I_{2}$
and

${\mathbb{E}}[\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}),\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}})^{\ast}]_{1}^{\mathfrak{d}}=1/4(I_{4}+K_{2}+\operatorname{vec}(I_{2})\operatorname{vec}(I_{2})^{\ast}).$

In this case, the chosen Lévy process $L$ satisfies Assumptions a from Section
4 (with $\sigma_{L}=1$ and $\sigma_{W}>0$). Based on the identification Lemma
5.14, we assume that the model is parameterized with
$\theta=(\theta^{(1)},\dots,\theta^{(11)})$, and the matrices
$A_{\theta},B_{\theta}$ and $C_{\theta}$ are defined as:

$A_{\theta}=\begin{pmatrix}\theta^{(1)}&\theta^{(2)}\\\
\theta^{(3)}&\theta^{(4)}\\\ \end{pmatrix},\quad
B_{\theta}=\begin{pmatrix}\theta^{(5)}&\theta^{(6)}\\\
\theta^{(6)}&\theta^{(7)}\\\ \end{pmatrix}\quad\text{and}\quad
C_{\theta}=\begin{pmatrix}\theta^{(8)}&\theta^{(9)}\\\
\theta^{(9)}&\theta^{(10)}\\\ \end{pmatrix},$

with $\theta^{(2)}>0$ and $\theta^{(2)}\not=\theta^{(3)}$. Thus, Assumption
$(d.4)$ and $(f.4)$ are automatically satisfied. The data used for estimation
is a sample of the log-price process ${\bm{G}}=({\bm{G}}_{i})_{i=1}^{n}$ as
defined in (1.2) with true parameter value
$\theta_{0}\in\Theta\subset{\mathbb{R}}^{11}$ observed on a fixed grid of size
$\Delta=0.1$ (the grid size for the Euler approximation of the Gaussian part
is $0.01$).

We experiment with two different settings, namely:

######  Example 6.1.

We fix $\sigma_{W}=1$,

$A_{\theta_{0}}=\begin{pmatrix}0.85&0.10\\\ -0.10&0.75\\\ \end{pmatrix},\quad
B_{\theta_{0}}=\begin{pmatrix}-2.43&0.05\\\ 0.05&-2.42\\\
\end{pmatrix}\quad\text{and}\quad C_{\theta_{0}}=\begin{pmatrix}1&0.5\\\
0.5&1.5\\\ \end{pmatrix}.$ (6.1)

######  Example 6.2.

We fix $\sigma_{W}=0,$ $A_{\theta_{0}}$ and $C_{\theta_{0}}$ are as in Example
6.1 and

$B_{\theta_{0}}=\frac{1}{4}\begin{pmatrix}-2.43&0.05\\\ 0.05&-2.42\\\
\end{pmatrix}.$ (6.2)

For the chosen Lévy process here, Assumption $(g.1)$ is satisfied. In Example
6.1, $\theta_{0}$ is chosen in such a way that the asymptotic normality of
$\hat{\theta}_{n}$ can be verified. Then, in Example 6.2 we rescale
$B_{\theta_{0}}$ from Example 6.1 in such a way that our sufficient conditions
for weak consistency are satisfied, but our sufficient conditions for
asymptotic normality in Corollary 5.17 are not satisfied.

Due to the identifiability Lemma 5.14 we need to choose $r\geq 2$. For
comparison purposes, we perform the estimation for maximum lags
$r\in\\{2,5,10\\}$ and sample sizes $n\in\\{1\,000,10\,000,100\,000\\}$. The
computations are performed with the optim routine in combination with the
Nelder-Mead algorithm in R (R Core Team (2017)). Initial values for the
estimation were found by the DEoptim routine on a neighborhood around the true
parameter $\theta_{0}$. We only consider estimators based on the identity
matrix for the weight matrix $\Omega$ in (5.3). The results are based on $500$
independent samples of MUCOGARCH returns.

In the following we report the finite sample results of the GMM for Examples
6.1 and 6.2.

###  6.1 Simulation results for Example 6.1

We can check numerically that the matrices $A_{\theta_{0}},B_{\theta_{0}}$ and
$A_{\theta_{0}}$ are such that Assumptions b and $(f.2)$-$(f.6)$ hold.
Additionally, the eigenvalues of the matrix
$B_{\theta_{0}}+B_{\theta_{0}}^{*}+\sigma_{L}A_{\theta_{0}}^{*}A_{\theta_{0}}$
are $-4.067$ and $-4.328$, so it is negative definite and Assumptions $(f.7)$
holds. For our choice of $\theta_{0}$ we have that $B_{\theta_{0}}$ is
diagonalizable with
$B_{\theta_{0}}=S_{\theta_{0}}D_{\theta_{0}}S_{\theta_{0}}^{-1}$, where

$S_{\theta_{0}}=\begin{pmatrix}-0.671&-0.741\\\ -0.741&0.671\\\
\end{pmatrix}\quad\text{and}\quad D_{\theta_{0}}=\begin{pmatrix}-2.375&0\\\
0&-2.475\\\ \end{pmatrix}.$

In addition, for $p=4.001$,

$\int_{{\mathbb{R}}^{2}}((1+\alpha_{\theta_{0}}\|\operatorname{vec}(yy^{\ast})\|_{S_{\theta_{0}}})^{p}-1)\nu_{L}({\rm
d}y)+2p\max\\{\mathfrak{R}(\sigma(B_{\theta_{0}}))\\}=-0.024<0,$ (6.3)

so $(g.2)$ is also valid. Therefore, all assumptions for applying Corollary
5.17 can be verified, which imply assumption e, and ensure asymptotic
normality. We also note that the chosen parameters are very close to not
satisfying Assumption (6.3).

We investigate the behavior of the bias and standard deviation in Figure 1,
where we excluded those paths for which the algorithm did not converge
successfully (around $10$ percent of the paths of length $n=1\,000$ and less
than $3$ percent for larger $n$). Figure 1 show the estimated absolute values
of the bias and standard deviation for different lags $r$ and varying $n$. As
expected, they decay when $n$ increases. Additionally, the results favor the
choice of maximum lag $r=10$, which is already expected since using more lags
of the autocovariance function usually helps to give a better fit. It is also
worth noting that the estimation of the parameters in the matrix
$B_{\theta_{0}}$ is more difficult than the other parameters, specially for
$n\in\\{1\,000,10\,000\\}$.

Figures 2 and 3 assess asymptotic normality though normal QQ-plots. Based on
the previous findings we fix $r=10$, since it gave the best results. This
might have to do with the fact that using just a few lags for the
autocovariance function ($r=2$ or $r=5$) are not sufficient for a good fit.
Here we do not exclude those paths for which the algorithm did not converge
(these are denoted by large red points in the normal QQ-plots in Figures 2 and
3). These plots are clearly in line with the asymptotic normality of the
estimators. It is worth noting that the tails corresponding to the estimates
of $B_{\theta_{0}}$ deviate from the ones of a normal distribution for values
of $n\in\\{1\,000,10\,000\\}$, but they get closer to a normal distribution
for $n=100\,000$. The tails of the plots for $A_{(2,1),\hat{\theta}_{n}}$ in
Figure 2 is not close to a normal (although the plots show its convergence).
This is maybe due to identifiability condition in Lemma 5.14 which requires
$A_{(2,1),\theta}>0$ but $A_{(2,1),\theta_{0}}=0.1$ is very close to the
boundary. For $n=1\,000$, there are very large negative outliers for the
estimates of $B_{\theta_{0}}$, which affects the bias substantially.

Figure 1: Example 1: Estimated absolute bias (lhs) and standard deviation
(rhs) of $\hat{\theta}_{n,r}$. The colors green, blue and red correspond to
$r=2,5$ and $10$, respectively.

Figure 2: Example 1: Normal QQ-plots of $\hat{\theta}_{n,10}$ for $\theta_{0}$
as in (6.2). The red dots are values for which the algorithm did not converge.

Figure 3: Example 1: Normal QQ-plots of $\hat{\theta}_{n,10}$ for $\theta_{0}$
as in (6.2). The red dots are values for which the algorithm did not converge.

Figure 4: Example 2: Estimated bias and std of $\hat{\theta}_{n,r}$. The
colors green, blue and red correspond to $r=2,5$ and $10$, respectively.

Figure 5: Example 2: Normal QQ-plots of $\hat{\theta}_{n,10}$ for $\theta_{0}$
as in (6.2). The red dots are values for which the algorithm did not converge.

Figure 6: Example 2: Normal QQ-plots of $\hat{\theta}_{n,10}$ for $\theta_{0}$
as in (6.2). The red dots are values for which the algorithm did not
converged.

###  6.2 Simulation results for Example 6.2

In this section we analyze the behavior of the GMM estimator when the
consistency conditions are valid, but we cannot check the conditions for
asymptotic normality. Here, we have
$\sigma(B_{\theta_{0}}/4+B_{\theta_{0}}^{*}/4+\sigma_{L}A_{\theta_{0}}^{*}A_{\theta_{0}})=\\{-0.594,-0.619\\}\in(-\infty,0)+i{\mathbb{R}}$.
Thus, Corollary 5.16 applies and gives weak consistency of the GMM estimator.
On the other hand, for $p=4.001$ the integral in (6.3) is $14.22>0$, and thus,
we cannot apply Corollary 5.17 to ensure asymptotic normality.

The results for Example 6.2 are given in Figures 4-6. The estimation of the
entries of $B_{\theta_{0}}$ does not seem to be substantially more difficult
than the entries of $A_{\theta_{0}}$ and $C_{\theta_{0}}$, as observed in the
previous example. Also, the estimated bias and std decreases in general as $n$
grows, showing consistency of the estimators. Also, the convergence rate seems
slow and, therefore, probably smaller than $n^{1/2}$ (the asymptotic normality
rate from Theorem 5.5). The QQ-plots for the estimation of the parameters
$A_{(2,1)}$, $C_{(1,1)}$ and $C_{(2,1)}$ also show some deviation from the
normal distribution.

##  7 Real data analysis

In this Section, we fit the MUCOGARCH model to 5 minutes log-returns of stock
prices corresponding to the SAP SE and Siemens AG companies (the data was
obtained from the Refinitiv EIKON system). For both datasets we have excluded
overnight returns. The resulting bivariate dataset has a total length of 12
135 (from 30-jun-2020 to 15-dez-2020) and is shown in Figure 7.

Figure 7: Log-returns (top) and sample acf of squared log-returns (bottom) for
SAP SE (left) and Siemens AG (right) exchange rates.

Based on the sample autocorrelation function of the squared log-returns we
decided to use $15$ lags to estimate a bivariate MUCOGARCH model. Before
moving to the estimation step, we re-scaled the data multiplying it by $1000$
for numerical reasons. Since the GMM algorithm from (5.3) requires a starting
value, we first fit a one-dimensional COGARCH to each dataset with the
algorithm described in (Haug et al., 2007, Algorithm 1) and used them to
construct the corresponding 2-dimensional MUCOGARCH describing two independent
log-return price process. This will be a MUCOGARCH process whose parameters
$A,B,C$ are diagonal matrices and the driving Lévy process has independent
components (Stelzer, 2010, Example 4.2). To construct a second step GMM
estimator, we use the estimated parameters from the first step to estimate a
weight matrix to replace $\Omega$ in (5.3). Indeed, as explained in (Mátyás,
1999, Section 1.3.3), and given a consistent estimator
$\hat{\Sigma}_{\theta_{0}}$ of $\Sigma_{\theta_{0}}$, an asymptotically more
efficient estimator can be constructed as

$\hat{\theta}^{(2)}_{n}=\operatorname*{arg\,min}_{\theta\in\Theta}\Big{\\{}(\hat{k}_{n,r}-k_{\theta,r})^{T}\hat{\Sigma}_{\theta_{0}}(\hat{k}_{n,r}-k_{\theta,r})\Big{\\}}.$
(7.1)

Given a bivariate dataset of log-price returns of size $n$ following a
MUCOGARCH model, we can see from (5.6) that a natural estimator for
$\Sigma_{\theta_{0}}$ is

$\begin{split}\hat{\Sigma}_{n,M}&=\frac{1}{n}\sum_{t=1}^{n}(D_{t}-\hat{k}_{n,r})(D_{t}-\hat{k}_{n,r})^{*}+\\\
&\frac{1}{n}\sum_{t=1}^{n}\sum_{i=1}^{M}\big{\\{}(D_{t}-\hat{k}_{n,r})(D_{t+i}-\hat{k}_{n,r})^{*}+(D_{t+i}-\hat{k}_{n,r})^{*}(D_{t}-\hat{k}_{n,r})\big{\\}},\quad
M\in\mathbb{N},\end{split}$ (7.2)

where we truncated the infinity sum in $\eqref{eq:defSigma}$. The above
estimator is a symmetric matrix, regardless of the values chosen for $M$ and
$n$. As done in the proof of Lemma 5.1, sufficient conditions for it to be
consistent are assumptions a, b, c and $(d.5)$. On the other hand, we cannot
guarantee that it is going to be positive semidefinite, which is a condition
required to use it as a weight matrix in the estimator from (7.1). One way
around it is to use the simpler estimator (as in (Stelzer et al., 2015,
Section 4.1)) of $\Sigma_{\theta_{0}}$, which ignores the second summation on
the rhs of (7.2), namely,

$\hat{\Sigma}_{n,M}^{(\text{basicCOV})}=\frac{1}{n}\sum_{t=1}^{n}(D_{t}-\hat{k}_{n,r})(D_{t}-\hat{k}_{n,r})^{*},\quad
M\in\mathbb{N},$ (7.3)

The resulting estimators when applying the 2 step GMM with (7.3) for the
bivariate dataset of 5min log-returns are given by:

$A=\begin{pmatrix}0.442&0.259\\\ 0.054&0.054\\\ \end{pmatrix},\quad
B=\begin{pmatrix}-0.146&-0.014\\\ -0.014&-0.080\\\
\end{pmatrix}\quad\text{and}\quad C=\begin{pmatrix}0.257&-0.134\\\
-0.134&0.070\\\ \end{pmatrix}.$ (7.4)

The parameters of the Levy process were chosen to be the same from example 1
as they allow us to check consistency and asymptotic normality conditions.
According to Corollary 5.16, the estimated parameters given in (7.4) describe
a MUCOGARCH model for which the GMM estimators is consistent. For the
asymptotic normality, we need to verify assumption (f.8), which requires
computing $m(p,\theta_{0})$ with $\theta_{0}$ replaced by a vector formed with
the entries of the matrices in (7.4). This computation results in 2.98, and
therefore (f.8) is violated, and we cannot ensure asymptotic normality.

To assess performance of the GMM estimator used to estimate the MUCOGARCH
model for real data, we perform a simulation study using the values from (7.4)
as $\theta_{0}$. The results are reported in Table 7.1 and referred by
basicCOV. We also compare it with the first step GMM and two other 2 step GMM
estimators: one using $\hat{\Sigma}_{n,M}$ from (7.2) with $M=10$ as a weight
matrix and another which uses a weighing diagonal matrix formed with the
diagonal of $\hat{\Sigma}_{n,M}$. These are referred in Table 7.1 as fullCOV
and diagCOV, respectively. Using fullCOV did not improved the estimates when
compared with basicCOV and diagCOV, and that might be to do with the fact that
inverting the estimated covariance matrix gave several warnings during the
estimation. The estimator based on diagCOV was somehow similar to the 1 step
GMM, giving smaller bias and std for the parameters
$\theta^{(6)}-\theta^{(8)}$, which correspond to the matrix $B$. The fact that
the 2 step GMM did not improved the 1step GMM results, might be because the
sample size $n$ used to estimate the covariance matrix is too small.

Estimated Bias - GMM
---
| True Value | 1 step | 2 step (BasicCOV) |  2 step ($M=10$) | 2 step (diagCOV)
$\theta^{(1)}$ | 0.442 | -0.024 | -0.094 | -0.055 | -0.073
$\theta^{(2)}$ | 0.259 | 0.016 | -0.026 | 0.040 | -0.034
$\theta^{(3)}$ | 0.054 | 0.086 | 0.120 | 0.121 | 0.100
$\theta^{(4)}$ | 0.194 | 0.164 | 0.170 | 0.151 | 0.169
$\theta^{(5)}$ | -0.146 | 0.010 | 0.014 | -0.043 | 0.031
$\theta^{(6)}$ | -0.014 | -0.085 | -0.076 | -0.091 | -0.072
$\theta^{(7)}$ | -0.080 | -0.081 | -0.077 | -0.162 | -0.058
$\theta^{(8)}$ | 0.257 | 0.091 | 0.129 | 0.248 | 0.083
$\theta^{(9)}$ | -0.134 | 0.147 | 0.148 | 0.176 | 0.148
$\theta^{(10)}$ | 0.070 | 0.136 | 0.157 | 0.315 | 0.137
Estimated Std - GMM
| True Value | 1 step | 2 step (BasicCOV) |  2 step ($M=10$) | 2 step (diagCOV)
$\theta^{(1)}$ | 0.442 | 0.101 | 0.084 | 0.207 | 0.100
$\theta^{(2)}$ | 0.259 | 0.123 | 0.112 | 0.261 | 0.112
$\theta^{(3)}$ | 0.054 | 0.088 | 0.086 | 0.169 | 0.093
$\theta^{(4)}$ | 0.194 | 0.109 | 0.108 | 0.193 | 0.110
$\theta^{(5)}$ | -0.146 | 0.080 | 0.091 | 0.402 | 0.067
$\theta^{(6)}$ | -0.014 | 0.073 | 0.081 | 0.398 | 0.065
$\theta^{(7)}$ | -0.080 | 0.072 | 0.085 | 0.650 | 0.062
$\theta^{(8)}$ | 0.257 | 0.112 | 0.113 | 0.909 | 0.168
$\theta^{(9)}$ | -0.134 | 0.093 | 0.101 | 0.645 | 0.139
$\theta^{(10)}$ | 0.070 | 0.091 | 0.103 | 1.137 | 0.126
Table 7.1: Estimated bias and std of $\hat{\theta}_{n,r}$ (1 step GMM) and
$\hat{\theta}_{n,r}^{(2)}$, the 2 step GMM with 3 different weighting
matrices, namely: basicCOV from (7.3), fullCOV from (7.2) with $M=10$ and
diagCOV, a diagonal weighting matrix formed with the diagonal entries of
(7.3).

#### Acknowledgement

Both authors take pleasure to thank Claudia Klüppelberg for motivating this
work and for helpful comments. Thiago do Rêgo Sousa gratefully acknowledges
support from the National Council for Scientific and Technological Development
(CNPq - Brazil) and the TUM Graduate School. He also thanks the Institute of
Mathematical Finance at Ulm University for its hospitality and takes pleasure
to thank Carlos Améndola for helpful discussions.

##  Appendix A Proofs

###  A.1 Auxiliary results

Several results related to the algebra of multivariate stochastic integrals
will be used here, for which we refer to Lemma 2.1 in Behme (2012).
Furthermore, we need the following.

######  Fact A.1 ((Stelzer, 2010, Lemma 6.9) with drift).

Assume that $(X_{t})_{t\in{\mathbb{R}}^{+}}$ is an adapted cadlag
$M_{d}({\mathbb{R}})$-valued process satisfying
${\mathbb{E}}(\|X_{t}\|)<\infty$ for all $t\in{\mathbb{R}}^{+}$,
$t\mapsto{\mathbb{E}}(\|X_{t}\|)$ is locally bounded and
$(L_{t})_{t\in{\mathbb{R}}_{+}}$ is an ${\mathbb{R}}^{d}$-valued Lévy process
of finite variation with ${\mathbb{E}}(\|L_{1}\|)<\infty$. Then

${\mathbb{E}}\int_{0}^{\Delta}X_{s-}{\rm
d}L_{s}=\int_{0}^{\Delta}{\mathbb{E}}(X_{s-}){\mathbb{E}}(L_{1}){\rm d}s.$

######  Fact A.2.

Let $(A_{t})_{t\in{\mathbb{R}}^{+}}$ in $M_{d^{2}}({\mathbb{R}})$,
$(B_{t})_{t\in{\mathbb{R}}^{+}}$ in $M_{1,d^{2}}({\mathbb{R}})$ be adapted
caglad processes satisfying ${\mathbb{E}}\|A_{t}\|\|B_{t}\|<\infty$ for all
$t\in{\mathbb{R}}^{+}$, $t\mapsto{\mathbb{E}}\|A_{t}\|\|B_{t}\|$ is locally
bounded and $(L_{t})_{t\in{\mathbb{R}}_{+}}$ be an ${\mathbb{R}}^{d}$ valued
Lévy process satisfying Assumption 5.2 in Stelzer (2010). Then,

${\mathbb{E}}\int_{0}^{t}A_{s}{\rm
d}(\operatorname{vec}([L,L]_{s}))B_{s}=(\sigma_{W}+\sigma_{L})\int_{0}^{t}{\mathbb{E}}[A_{s}\operatorname{vec}(I_{d})B_{s}]{\rm
d}s.$

###### Proof.

First notice that $\operatorname{vec}([L,L]_{s})$ is an
${\mathbb{R}}^{d^{2}}$-valued Lévy process with finite variation. Then it
follows from Fact A.1 that

$\begin{split}&\operatorname{vec}\bigg{(}{\mathbb{E}}\int_{0}^{t}A_{s}{\rm
d}(\operatorname{vec}([L,L]_{s}))B_{s}\bigg{)}={\mathbb{E}}\int_{0}^{t}(B_{s}\otimes
A_{s}){\rm d}(\operatorname{vec}([L,L]_{s})\\\
&=\int_{0}^{t}{\mathbb{E}}(B_{s}\otimes
A_{s}){\mathbb{E}}(\operatorname{vec}([L,L]_{1})){\rm
d}s=(\sigma_{W}+\sigma_{L})\int_{0}^{t}{\mathbb{E}}(B_{s}\otimes
A_{s})\operatorname{vec}(I_{d}){\rm d}s\\\
&=(\sigma_{W}+\sigma_{L})\operatorname{vec}\bigg{(}\int_{0}^{t}{\mathbb{E}}(A_{s}I_{d}B_{s}){\rm
d}s\bigg{)},\end{split}$

so the result follows by an application of $\operatorname{vec}^{-1}$. ∎

Proof of Lemma 3.2: It follows from (Stelzer, 2010, Proposition 4.7) (with
$k=p$) that ${\mathbb{E}}\|Y_{t}\|^{p}<\infty$ for all $t\in{\mathbb{R}}^{+}$
and $t\mapsto{\mathbb{E}}\|Y_{t}\|^{p}$ is locally bounded. Then an
application of (Protter, 2005, Theorem 66 of Ch. 5) together with the fact
that ${\mathbb{E}}\|L_{1}\|^{2p}<\infty$ and the definition of
$(V_{t})_{t\in{\mathbb{R}}_{+}}$ in (3.2) gives for all $t>0$

$\begin{split}{\mathbb{E}}\|G_{t}\|^{2p}&={\mathbb{E}}\bigg{\|}\int_{0}^{t}V_{s-}^{1/2}{\rm
d}L_{s}\bigg{\|}^{2p}\leq c\int_{0}^{t}{\mathbb{E}}\|V_{s-}^{1/2}\|^{2p}{\rm
d}s\leq c\int_{0}^{t}{\mathbb{E}}\|C+Y_{s-}\|^{p}{\rm d}s.\end{split}$

######  Lemma A.1.

Assume that Assumptions $(a.1)$-$(a.4)$, b and $(c.2)$ hold. Then,

$\begin{split}\quad\quad&{\rm
cov}(\operatorname{vec}(Y_{\Delta}),\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast}))={\rm
cov}(\operatorname{vec}(Y_{\Delta}),\operatorname{vec}(G_{\Delta}G_{\Delta}^{\ast}))\\\
&={\rm
var}(\operatorname{vec}(V_{0}))(e^{{\mathcal{B}}^{\ast}\Delta}-I_{d^{2}})[(\sigma_{W}+\sigma_{L})({\mathcal{B}}^{\ast})^{-1}-2((A\otimes
A)^{\ast})^{-1}],\quad\Delta\geq 0.\end{split}$ (A.1)

###### Proof.

Since $(a.4)$, $(a.1)$ and $(c.1)$ hold, we can apply Lemma 3.2 with $p=2$ to
conclude that both $\|\operatorname{vec}(Y_{\Delta})\|$ and
$\|\bm{G}_{1}\bm{G}_{1}^{\ast}\|$ are square integrable random variables and
thus, the covariance at the left hand side of (A.1) is finite. Integration by
parts formula (Stelzer, 2010, p. 111) gives

$G_{\Delta}G_{\Delta}^{\ast}=\int_{0}^{\Delta}V_{s-}^{1/2}{\rm
d}L_{s}G_{s-}^{\ast}+\int_{0}^{\Delta}G_{s-}{\rm
d}L_{s}^{\ast}V_{s-}^{1/2}+\int_{0}^{\Delta}V_{s-}^{1/2}{\rm
d}[L,L^{\ast}]_{s}V_{s-}^{1/2}:=A_{\Delta}+A_{\Delta}^{\ast}+C_{\Delta}.$
(A.2)

It follows from Lemma 3.2(a) and (b) with $p=2$ together with the Cauchy-
Schwarz inequality that

$\int_{0}^{t}{\mathbb{E}}(\|V_{s-}^{1/2}\|_{2}\|G_{s-}\|_{2})^{2}{\rm
d}s\leq\int_{0}^{t}\big{(}{\mathbb{E}}\|V_{s-}\|_{2}^{2}\big{)}^{1/2}\big{(}{\mathbb{E}}\|G_{s-}\|_{2}^{4}\big{)}^{1/2}{\rm
d}s<\infty,$ (A.3)

where the finiteness is due to the fact that the integrand is locally bounded,
and thus, also bounded on $(0,t)$. Therefore $(A_{t})_{t\in{\mathbb{R}}^{+}}$
is a martingale and $A_{t}\in L^{2}$ for all $t\geq 0$. Thus, the integration
by parts formula, the formula ${\rm d}(\operatorname{vec}(A_{s}))^{\ast}={\rm
d}L_{s}^{*}(G_{s-}^{\ast}\otimes V_{s-}^{1/2})$ (Lemma 2.1(vi) in Behme
(2012)) imply

$\begin{split}\quad\,\,&{\rm
cov}(\operatorname{vec}(Y_{\Delta}),\operatorname{vec}(A_{\Delta}))\\\
&={\mathbb{E}}\big{(}\operatorname{vec}(Y_{\Delta})(\operatorname{vec}(A_{\Delta}))^{\ast}\big{)}-{\mathbb{E}}(\operatorname{vec}(Y_{\Delta})){\mathbb{E}}(\operatorname{vec}(A_{\Delta}))^{\ast}\\\
&={\mathbb{E}}\bigg{(}\int_{0}^{\Delta}\operatorname{vec}(Y_{s-}){\rm
d}(\operatorname{vec}(A_{s}))^{\ast}+\int_{0}^{\Delta}{\rm
d}\operatorname{vec}(Y_{s})(\operatorname{vec}(A_{s-}))^{\ast}+[\operatorname{vec}(Y),(\operatorname{vec}(A))^{\ast}]_{\Delta}\bigg{)}-0\\\
&=0+{\mathbb{E}}\int_{0}^{\Delta}{\rm
d}\operatorname{vec}(Y_{s})(\operatorname{vec}(A_{s-}))^{\ast}+{\mathbb{E}}([\operatorname{vec}(Y),(\operatorname{vec}(A))^{\ast}]_{\Delta}).\end{split}$
(A.4)

The first expectation in (A.4) vanishes since

$\begin{split}&\int_{0}^{\Delta}{\mathbb{E}}\|\operatorname{vec}(Y_{s-})\|^{2}\|G_{s-}\|^{2}\|V_{s-}^{1/2}\|^{2}{\rm
d}s\\\
&\leq\int_{0}^{\Delta}({\mathbb{E}}\|\operatorname{vec}(Y_{s-})\|^{4})^{1/2}({\mathbb{E}}\|G_{s-}\|^{8})^{1/4}({\mathbb{E}}\|V_{s-}^{1/2}\|^{8})^{1/4}{\rm
d}s<\infty\end{split}$

by the generalized Hölder inequality with $(1/2+1/4+1/4=1)$ (see e.g. (Kufner
et al., 1977, Theorem 2.1)), Lemma 3.2 and the fact that
$(L_{t})_{t\in{\mathbb{R}}_{+}}$ is an $L^{2}$-martingale. Let
$\tilde{C}:=(B\otimes I+I\otimes B)$ and recall from p. 84 in Stelzer (2010)
that

${\rm d}\operatorname{vec}(Y_{s})=\tilde{C}\operatorname{vec}(Y_{s-}){\rm
d}s+(A\otimes A)(V_{s-}^{1/2}\otimes V_{s-}^{1/2}){\rm
d}\operatorname{vec}([L,L]_{s}^{\mathfrak{d}}).$ (A.5)

Using (A.5), the bilinearity of the quadratic covariation process, (Behme,
2012, eq. (2.1)), Lemma 3.2, Facts A.1, $(a.3)$, (A.3) and the Itô isometry we
obtain

$\displaystyle[\operatorname{vec}(Y),(\operatorname{vec}(A))^{\ast}]_{\Delta}$
(A.6) $\displaystyle=$
$\displaystyle\bigg{[}\int_{0}^{\cdot}\tilde{C}\operatorname{vec}(Y_{s-}){\rm
d}s+\int_{0}^{\cdot}(A\otimes A)(V_{s-}^{1/2}\otimes V_{s-}^{1/2}){\rm
d}\operatorname{vec}([L,L]_{s}^{\mathfrak{d}}),\int_{0}^{\cdot}{\rm
d}L_{s}^{*}(G_{s-}^{\ast}\otimes V_{s-}^{1/2})\bigg{]}_{\Delta}$
$\displaystyle=$ $\displaystyle\int_{0}^{\Delta}(A\otimes
A)(V_{s-}^{1/2}\otimes V_{s-}^{1/2}){\rm
d}[\operatorname{vec}([L,L]^{\mathfrak{d}}),L^{\ast}]_{s}(G_{s-}^{\ast}\otimes
V_{s-}^{1/2}).$

Recall that for arbitrary matrices $M\in M_{m,n}({\mathbb{R}})$ and $N\in
M_{k,l}({\mathbb{R}})$ it holds $\|A\otimes B\|_{2}=\|A\|_{2}\|B\|_{2}$
(Bernstein, 2009, Fact 9.9.61). This together with the Hölder inequality with
(3/4 + 1/4 = 1) and Lemma 3.2 with $p=4$ gives

$\begin{split}&\int_{0}^{\Delta}{\mathbb{E}}\|V_{s-}^{1/2}\otimes
V_{s-}^{1/2}\|_{2}\|(G_{s-}^{\ast}\otimes V_{s-}^{1/2})\|_{2}{\rm
d}s=\int_{0}^{\Delta}{\mathbb{E}}\|V_{s-}^{1/2}\|_{2}^{3}\|G_{s-}^{\ast}\|_{2}{\rm
d}s\\\
&=\int_{0}^{\Delta}{\mathbb{E}}\|V_{s-}\|_{2}^{3/2}\|G_{s-}^{\ast}\|_{2}{\rm
d}s\leq\int_{0}^{\Delta}({\mathbb{E}}\|V_{s-}\|_{2}^{2})^{3/4}({\mathbb{E}}\|G_{s-}^{\ast}\|_{2}^{4})^{1/4}{\rm
d}s<\infty.\end{split}$

Thus, applying expectations at both sides of (A.6) gives

$E[\operatorname{vec}(Y),(\operatorname{vec}(A))^{\ast}]_{\Delta}=0.$ (A.7)

Let
$l_{s}:={\mathbb{E}}\operatorname{vec}(Y_{s})(\operatorname{vec}(A_{s}))^{\ast}$
and notice that it follows from Lemma 3.2 and the Cauchy-Schwarz inequality
that ${\mathbb{E}}\|l_{s}\|<\infty$ and $s\mapsto{\mathbb{E}}\|l_{s}\|$ is
locally bounded. Use (A.5), (A.7), the compensation formula, (Bernstein, 2009,
Proposition 7.1.9),
${\mathbb{E}}\operatorname{vec}(V_{s})\operatorname{vec}(A_{s-})=l_{s}$, the
Itô isometry and $(a.2)$ to get

$\begin{split}l_{\Delta}&={\mathbb{E}}\int_{0}^{\Delta}{\rm
d}\operatorname{vec}(Y_{s})(\operatorname{vec}(A_{s-}))^{\ast}\\\
&={\mathbb{E}}\int_{0}^{\Delta}\big{[}\tilde{C}\operatorname{vec}(Y_{s-}){\rm
d}s+(A\otimes A)(V_{s-}^{1/2}\otimes V_{s-}^{1/2}){\rm
d}(\operatorname{vec}([L,L]_{s}^{\mathfrak{d}})\big{]}(\operatorname{vec}(A_{s-}))^{\ast}\\\
&=\tilde{C}\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s-})(\operatorname{vec}(A_{s-}))^{\ast}{\rm
d}s\\\ &\quad\quad+\sigma_{L}\int_{0}^{\Delta}{\mathbb{E}}[(A\otimes
A)(V_{s-}^{1/2}\otimes
V_{s-}^{1/2})\operatorname{vec}(I_{d})(\operatorname{vec}(A_{s-}))^{\ast}]{\rm
d}s\\\ &=(\tilde{C}+\sigma_{L}(A\otimes A))\int_{0}^{\Delta}l_{s}{\rm
d}s.\end{split}$ (A.8)

Solving the matrix-valued integral equation in (A.8) and using that $A_{0}=0$
implies $l_{0}=0$, gives $l_{s}=0$ for all $s\geq 0$ (see Haug et al. (2007)).
Thus, it follows from (A.4)-(A.8) that

${\rm cov}(\operatorname{vec}(Y_{\Delta}),\operatorname{vec}(A_{\Delta}))=0,$
(A.9)

and, as a consequence ${\rm
cov}(\operatorname{vec}(Y_{\Delta}),\operatorname{vec}(A_{\Delta}^{\ast}))=0$.
Let ${\mathcal{V}}_{s-}:=V_{s-}^{1/2}\otimes V_{s-}^{1/2}$. Then,

$\displaystyle\operatorname{vec}(C_{\Delta})=\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]_{s})=\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}_{s})$ (A.10)
$\displaystyle\quad\quad+\sigma_{W}\int_{0}^{\Delta}(V_{s-}^{1/2}\otimes
V_{s-}^{1/2})\operatorname{vec}(I_{d}){\rm
d}s=\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}_{s})+\sigma_{W}\int_{0}^{\Delta}\operatorname{vec}(V_{s-}){\rm
d}s.$

Using the compensation formula, Fact A.1 and the stationarity of
$(V_{s})_{s\in{\mathbb{R}}_{+}}$ we get

$\begin{split}{\mathbb{E}}\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]_{s})&=(\sigma_{W}+\sigma_{L})\int_{0}^{\Delta}{\mathbb{E}}{\mathcal{V}}_{s-}\operatorname{vec}(I_{d}){\rm
d}s=\Delta(\sigma_{W}+\sigma_{L}){\mathbb{E}}\operatorname{vec}(V_{0}).\end{split}$
(A.11)

Additionally, it follows from Lemma 3.2 that
${\mathbb{E}}\|\operatorname{vec}(V_{s})\operatorname{vec}(Y_{\Delta})^{*}\|<\infty$
for all $s\geq 0$ and that
$s\mapsto{\mathbb{E}}\|\operatorname{vec}(V_{s})\operatorname{vec}(Y_{\Delta})^{*}\|$
is locally bounded. Then,

$\begin{split}\quad\quad&{\mathbb{E}}\bigg{(}\int_{0}^{\Delta}\operatorname{vec}(V_{s-}){\rm
d}s\,(\operatorname{vec}(Y_{\Delta}))^{\ast}\bigg{)}=\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(V_{s-})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s\\\
&=\Delta\operatorname{vec}(C)({\mathbb{E}}\operatorname{vec}(Y_{0}))^{\ast}+\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s.\end{split}$

Now it follows from the invertibility of $(A\otimes A)$ and from the second
equation following (3.5) in Stelzer (2010) that

$\begin{split}\quad\quad&\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}_{s})\\\ &=(A\otimes
A)^{-1}\bigg{(}\operatorname{vec}(Y_{\Delta})-\operatorname{vec}(Y_{0})-\int_{0}^{\Delta}(B\otimes
I+I\otimes B)\operatorname{vec}(Y_{s-}){\rm d}s\bigg{)}.\end{split}$ (A.12)

The representation in (A.12) gives

$\displaystyle{\mathbb{E}}\bigg{[}\bigg{(}\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}_{s})\bigg{)}(\operatorname{vec}(Y_{\Delta}))^{\ast}\bigg{]}$
(A.13) $\displaystyle=$ $\displaystyle{\mathbb{E}}\bigg{[}(A\otimes
A)^{-1}\Big{(}\operatorname{vec}(Y_{\Delta})-\operatorname{vec}(Y_{0})-\int_{0}^{\Delta}(B\otimes
I+I\otimes B)\operatorname{vec}(Y_{s-}){\rm
d}s\Big{)}(\operatorname{vec}(Y_{\Delta}))^{\ast}\bigg{]}$ $\displaystyle=$
$\displaystyle(A\otimes
A)^{-1}\bigg{[}{\mathbb{E}}\operatorname{vec}(Y_{\Delta})(\operatorname{vec}(Y_{\Delta}))^{\ast}-{\mathbb{E}}\operatorname{vec}(Y_{0})(\operatorname{vec}(Y_{\Delta}))^{\ast}$
$\displaystyle\quad-(B\otimes I+I\otimes
B)\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s-})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s\bigg{]}.$

Using the definition of $C_{\Delta}$ in (A.2), together with (A.10), (A.11)
and (A.13) gives

$\displaystyle{\rm
cov}(\operatorname{vec}(C_{\Delta}),\operatorname{vec}(Y_{\Delta}))=(A\otimes
A)^{-1}\bigg{[}{\mathbb{E}}\operatorname{vec}(Y_{\Delta})(\operatorname{vec}(Y_{\Delta}))^{\ast}-{\mathbb{E}}\operatorname{vec}(Y_{0})(\operatorname{vec}(Y_{\Delta}))^{\ast}$
$\displaystyle-(B\otimes I+I\otimes
B)\bigg{(}\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s-})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s\bigg{)}\bigg{]}$
$\displaystyle+\Delta\sigma_{W}\operatorname{vec}(C)({\mathbb{E}}\operatorname{vec}(Y_{0}))^{\ast}+\sigma_{W}\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s-})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s$
$\displaystyle-\Delta(\sigma_{W}+\sigma_{L}){\mathbb{E}}\operatorname{vec}(V_{0}){\mathbb{E}}(\operatorname{vec}(Y_{\Delta}))^{\ast}$
$\displaystyle=[\sigma_{W}I_{d^{2}}-(A\otimes A)^{-1}(B\otimes I+I\otimes
B)]\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s$ $\displaystyle+(A\otimes A)^{-1}\big{[}{\rm
var}(\operatorname{vec}(Y_{0}))-{\rm
cov}(\operatorname{vec}(Y_{0}),\operatorname{vec}(Y_{\Delta}))\big{]}-\Delta\sigma_{L}\operatorname{vec}(C){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}$
$\displaystyle-\Delta(\sigma_{W}+\sigma_{L}){\mathbb{E}}\operatorname{vec}(Y_{0}){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast},$

where the last equality follows from $V_{0}=C+Y_{0}$ and the stationarity of
$(Y_{s})_{s\in{\mathbb{R}}_{+}}$. Using (4.2) it follows first that

$\begin{split}\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(Y_{s})(\operatorname{vec}(Y_{\Delta}))^{\ast}{\rm
d}s&=\int_{0}^{\Delta}e^{{\mathcal{B}}(\Delta-s)}{\rm
var}(\operatorname{vec}(Y_{0})){\rm
d}s+\Delta{\mathbb{E}}\operatorname{vec}(Y_{0}){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}\\\
&={\mathcal{B}}^{-1}(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0}))+\Delta{\mathbb{E}}\operatorname{vec}(Y_{0}){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast},\end{split}$
(A.14)

and second that

${\rm var}(\operatorname{vec}(Y_{0}))-{\rm
cov}(\operatorname{vec}(Y_{0}),\operatorname{vec}(Y_{\Delta}))=-(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0})).$ (A.15)

Substituting $B\otimes I+I\otimes B={\mathcal{B}}-\sigma_{L}(A\otimes A)$,
using (A.14), (A.15), (4.1) and the formula for
${\mathbb{E}}\operatorname{vec}(Y_{0})$ in (4.2) gives

$\begin{split}\quad\quad&{\rm
cov}(\operatorname{vec}(C_{\Delta}),\operatorname{vec}(Y_{\Delta}))\\\
&=\big{[}\sigma_{W}I_{d^{2}}-(A\otimes
A)^{-1}({\mathcal{B}}-\sigma_{L}(A\otimes
A))\big{]}\big{[}{\mathcal{B}}^{-1}(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0}))\\\
&+\Delta{\mathbb{E}}\operatorname{vec}(Y_{0}){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}\big{]}\\\
&-(A\otimes A)^{-1}(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0}))-\Delta\sigma_{L}\operatorname{vec}(C){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}\\\
&-\Delta(\sigma_{W}+\sigma_{L}){\mathbb{E}}\operatorname{vec}(Y_{0}){\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}\\\
&=\big{[}(\sigma_{W}+\sigma_{L}){\mathcal{B}}^{-1}-2(A\otimes
A)^{-1}\big{]}(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0}))\\\ &\quad\quad-\big{[}(A\otimes
A)^{-1}{\mathcal{B}}{\mathbb{E}}\operatorname{vec}(Y_{0})+\sigma_{L}\operatorname{vec}(C)\big{]}\Delta{\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}\\\
&=\big{[}(\sigma_{W}+\sigma_{L}){\mathcal{B}}^{-1}-2(A\otimes
A)^{-1}\big{]}(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0}))\\\ &\quad\quad-\big{[}(A\otimes
A)^{-1}{\mathcal{B}}(-\sigma_{L}{\mathcal{B}}^{-1}(A\otimes
A)\operatorname{vec}(C))+\sigma_{L}\operatorname{vec}(C)\big{]}\Delta{\mathbb{E}}(\operatorname{vec}(Y_{0}))^{\ast}\\\
&=\big{[}(\sigma_{W}+\sigma_{L}){\mathcal{B}}^{-1}-2(A\otimes
A)^{-1}\big{]}(e^{{\mathcal{B}}\Delta}-I_{d^{2}}){\rm
var}(\operatorname{vec}(Y_{0})).\end{split}$ (A.16)

Finally, the result of the Lemma follows from (A.2), (A.9), (A.16) and the
fact that

${\rm
cov}(\operatorname{vec}(Y_{\Delta}),\operatorname{vec}(G_{\Delta}G_{\Delta}^{\ast}))=({\rm
cov}(\operatorname{vec}(G_{\Delta}G_{\Delta}^{\ast}),\operatorname{vec}(Y_{\Delta})))^{\ast}=({\rm
cov}(\operatorname{vec}(C_{\Delta}),\operatorname{vec}(Y_{\Delta})))^{\ast}.$

∎

###  A.2 Proof of Lemma 4.2

(i) The proof of Lemma 4.2 (i) follows directly from Lemma A.1 combined with
(5.7) in Stelzer (2010).
(ii) Denoting by $\|\cdot\|_{F}$ the Frobenius norm we have by Lemma 3.2(b)
with $p=2$

$\displaystyle{\mathbb{E}}\|\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})^{\ast}\|_{2}$
$\displaystyle=$
$\displaystyle{\mathbb{E}}\|\operatorname{vec}(\bm{G}_{1}\bm{G}_{1}^{\ast})\|_{2}^{2}={\mathbb{E}}\|\bm{G}_{1}\bm{G}_{1}^{\ast}\|^{2}_{F}$
$\displaystyle=$ $\displaystyle
tr(\bm{G}_{1}\bm{G}_{1}^{\ast}\bm{G}_{1}\bm{G}_{1}^{\ast})={\mathbb{E}}\|\bm{G}_{1}\|_{2}^{4}<\infty.$

Let $a_{s}:=\operatorname{vec}(G_{s}G_{s}^{\ast}),s\in[0,\Delta]$ and use the
integration by parts formula to write

$\begin{split}a_{\Delta}a_{\Delta}^{\ast}&=\int_{0}^{\Delta}a_{s-}{\rm
d}(a_{s}^{\ast})+\int_{0}^{\Delta}{\rm
d}a_{s}(a_{s-}^{\ast})+[a,a^{\ast}]_{\Delta}\\\
&=\bigg{(}\int_{0}^{\Delta}{\rm
d}a_{s}(a_{s-}^{\ast})\bigg{)}^{\ast}+\int_{0}^{\Delta}{\rm
d}a_{s}(a_{s-}^{\ast})+[a,a^{\ast}]_{\Delta}\end{split}$ (A.17)

hence we only need to prove that the random variables

$\int_{0}^{\Delta}{\rm
d}a_{s}(a_{s-}^{\ast})\quad\text{and}\quad[a,a^{\ast}]_{\Delta},$

have finite expectations and compute them in closed form. From (A.2), Lemma
2.1(vi) in Behme (2012) and the symmetry of $(V_{t})_{t\in{\mathbb{R}}_{+}}$
it follows that

$\displaystyle{\rm d}a_{t}$ $\displaystyle=$ $\displaystyle{\rm
d}(\operatorname{vec}(G_{t}G_{t}^{\ast}))$ (A.18) $\displaystyle=$
$\displaystyle{\rm
d}\bigg{(}\operatorname{vec}\bigg{(}\int_{0}^{t}V_{s-}^{1/2}{\rm
d}L_{s}G_{s-}^{\ast}+\int_{0}^{t}G_{s-}{\rm
d}L_{s}^{\ast}V_{s-}^{1/2}+\int_{0}^{t}V_{s-}^{1/2}{\rm
d}[L,L^{\ast}]_{s}V_{s-}^{1/2}\bigg{)}\bigg{)}$ $\displaystyle=$
$\displaystyle{\rm d}\bigg{(}\int_{0}^{t}(G_{s-}\otimes V_{s-}^{1/2}){\rm
d}L_{s}+\int_{0}^{t}(V_{s-}^{1/2}\otimes G_{s-}){\rm
d}L_{s}+\int_{0}^{t}(V_{s-}^{1/2}\otimes V_{s-}^{1/2}){\rm
d}\operatorname{vec}([L,L^{\ast}]_{s})\bigg{)}$ $\displaystyle=$
$\displaystyle(G_{t-}\otimes V_{t-}^{1/2}+V_{t-}^{1/2}\otimes G_{t-}){\rm
d}L_{t}+(V_{t-}^{1/2}\otimes V_{t-}^{1/2}){\rm
d}\operatorname{vec}([L,L^{\ast}]_{t}),\quad t\geq 0.$

By the sub-multiplicative property of $\|\cdot\|_{2}$, the generalized Hölder
inequality with $(1/4+1/4+1/2=1)$ we have

$\displaystyle\int_{0}^{\Delta}{\mathbb{E}}\|G_{s-}\otimes
V_{s-}^{1/2}\|_{2}^{2}\|a_{s-}\|_{2}^{2}{\rm
d}s=\int_{0}^{\Delta}{\mathbb{E}}\|G_{s-}\|_{2}^{2}\|V_{s-}^{1/2}\|_{2}^{2}\|\operatorname{vec}(G_{s-}G_{s-}^{\ast})\|_{2}^{2}{\rm
d}s$ $\displaystyle=$
$\displaystyle\int_{0}^{\Delta}{\mathbb{E}}\|G_{s-}\|_{2}^{2}\|V_{s-}^{1/2}\|_{2}^{2}\|G_{s-}G_{s-}^{\ast}\|_{2}^{2}{\rm
d}s\leq\int_{0}^{\Delta}({\mathbb{E}}\|G_{s-}\|_{2}^{8})^{1/4}({\mathbb{E}}\|V_{s-}\|_{2}^{4})^{1/4}({\mathbb{E}}\|G_{s-}\|_{2}^{8})^{1/2}{\rm
d}s,$

which is finite by Lemma 3.2 with $p=4$. Additionally, similar calculations
and Lemma 3.2 with $p=2$ shows that ${\mathbb{E}}(\|V_{s-}^{1/2}\otimes
V_{s-}^{1/2}\|_{2})\|a_{s-}\|_{2}\leq({\mathbb{E}}\|V_{s-}\|_{2}^{2})^{1/2}({\mathbb{E}}\|G_{s-}\|_{2}^{4})^{1/2}<\infty$
for all $s>0$ and the map $s\mapsto{\mathbb{E}}(\|V_{s-}^{1/2}\otimes
V_{s-}^{1/2}\|_{2}\|a_{s-}\|_{2})$ is locally bounded. Thus it follows from
(A.18), the Itô isometry, the fact that
$[L,L^{\ast}]_{t}=[L,L^{\ast}]_{t}^{\mathfrak{d}}+\sigma_{w}I_{d}t$ and fact
A.2 that

$\begin{split}\quad\,\,&{\mathbb{E}}\int_{0}^{\Delta}{\rm
d}a_{s}(a_{s-}^{\ast})\\\
&={\mathbb{E}}\bigg{(}\int_{0}^{\Delta}(G_{s-}\otimes
V_{s-}^{1/2}+V_{s-}^{1/2}\otimes G_{s-}){\rm
d}L_{s}a_{s-}^{\ast}+\int_{0}^{\Delta}(V_{s-}^{1/2}\otimes V_{s-}^{1/2}){\rm
d}(\operatorname{vec}([L,L^{\ast}]_{s})a_{s-}^{\ast}\bigg{)}\\\
&=(\sigma_{L}+\sigma_{W})\bigg{(}\int_{0}^{\Delta}{\mathbb{E}}\big{(}(V_{s-}^{1/2}\otimes
V_{s-}^{1/2})\operatorname{vec}(I_{d})a_{s-}^{\ast}\big{)}{\rm d}s\bigg{)}\\\
&=(\sigma_{L}+\sigma_{W})\int_{0}^{\Delta}{\mathbb{E}}(\operatorname{vec}(V_{s-})a_{s-}^{\ast}){\rm
d}s.\end{split}$ (A.20)

It follows from (5.6) in Stelzer (2010) that

$\int_{0}^{\Delta}{\mathbb{E}}a_{s-}^{\ast}{\rm
d}s=\int_{0}^{\Delta}\big{(}\operatorname{vec}((\sigma_{L}+\sigma_{W})s{\mathbb{E}}V_{0})\big{)}^{\ast}{\rm
d}s=\frac{1}{2}(\sigma_{L}+\sigma_{W})\Delta^{2}{\mathbb{E}}\operatorname{vec}(V_{0})^{\ast}.$
(A.21)

Since we assumed here that all hypothesis for using Lemma A.1 are valid, we
can use (A.1) with $\Delta=s$ to get

$\begin{split}\quad\quad&\int_{0}^{\Delta}{\rm
cov}(\operatorname{vec}(Y_{s-}),a_{s-}){\rm d}s\\\ &={\rm
var}(\operatorname{vec}(Y_{0}))\bigg{(}\int_{0}^{\Delta}(e^{{\mathcal{B}}^{\ast}s}-I_{d^{2}}){\rm
d}s\bigg{)}\big{[}(\sigma_{W}+\sigma_{L})({\mathcal{B}}^{\ast})^{-1}-2((A\otimes
A)^{\ast})^{-1}\big{]}\\\ &={\rm
var}(\operatorname{vec}(Y_{0}))\tilde{{\mathcal{B}}},\end{split}$ (A.22)

where $\tilde{{\mathcal{B}}}$ is defined in (4.7). Using (A.20), $(c.1)$
(A.21), (A.22) gives

$\begin{split}\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(V_{s-})a_{s-}^{\ast}{\rm
d}s&=\int_{0}^{\Delta}{\rm cov}(\operatorname{vec}(V_{s-}),a_{s-}){\rm
d}s+({\mathbb{E}}\operatorname{vec}(V_{s}))\int_{0}^{\Delta}{\mathbb{E}}(a_{s-}^{\ast}){\rm
d}s\\\ &=\int_{0}^{\Delta}{\rm cov}(\operatorname{vec}(Y_{s-}),a_{s-}){\rm
d}s+({\mathbb{E}}\operatorname{vec}(V_{0}))\int_{0}^{\Delta}{\mathbb{E}}(a_{s-}^{\ast}){\rm
d}s\\\
&=\frac{1}{2}(\sigma_{L}+\sigma_{W})\Delta^{2}{\mathbb{E}}\operatorname{vec}(V_{0})E\operatorname{vec}(V_{0})^{\ast}+{\rm
var}(\operatorname{vec}(Y_{0}))\tilde{{\mathcal{B}}}\\\
&=(\sigma_{L}+\sigma_{W})^{-1}D,\end{split}$ (A.23)

where $D$ is defined in (4.6). Let $f_{s}:=(G_{s-}\otimes
V_{s-}^{1/2}+V_{s-}^{1/2}\otimes G_{s-}),s\geq 0$ and recall
${\mathcal{V}}_{s-}=V_{s-}^{1/2}\otimes V_{s-}^{1/2}$. Using (A.2), Lemma
2.1(vi) in Behme (2012) and the symmetry of $V_{s-}^{1/2}$ gives

$\begin{split}\quad\quad\,\,&[a,a^{\ast}]_{\Delta}\\\
&=\bigg{[}\operatorname{vec}\bigg{(}\int_{0}^{\cdot}V_{s-}^{1/2}{\rm
d}L_{s}G_{s-}^{\ast}+\int_{0}^{\cdot}G_{s-}{\rm
d}L_{s}^{\ast}V_{s-}^{1/2}+\int_{0}^{\cdot}V_{s-}^{1/2}{\rm
d}[L,L^{\ast}]_{s}V_{s-}^{1/2}\bigg{)},\\\
&\bigg{(}\operatorname{vec}\bigg{(}\int_{0}^{\cdot}V_{s-}^{1/2}{\rm
d}L_{s}G_{s-}^{\ast}+\int_{0}^{\cdot}G_{s-}{\rm
d}L_{s}^{\ast}V_{s-}^{1/2}+\int_{0}^{\cdot}V_{s-}^{1/2}{\rm
d}[L,L^{\ast}]_{s}V_{s-}^{1/2}\bigg{)}\bigg{)}^{\ast}\bigg{]}_{\Delta}\\\
&=\bigg{[}\int_{0}^{\cdot}f_{s-}{\rm
d}L_{s}+\int_{0}^{\cdot}{\mathcal{V}}_{s-}{\rm
d}\operatorname{vec}([L,L^{\ast}]_{s}),\int_{0}^{\cdot}{\rm
d}L_{s}^{\ast}f_{s-}^{\ast}+\int_{0}^{\cdot}{\rm
d}(\operatorname{vec}([L,L^{\ast}]_{s})^{\ast}){\mathcal{V}}_{s-}\bigg{]}_{\Delta}\\\
&=\int_{0}^{\Delta}f_{s-}{\rm
d}[L,L^{\ast}]_{s}f_{s-}^{\ast}+\int_{0}^{\Delta}f_{s-}{\rm
d}[L,\operatorname{vec}([L,L^{\ast}])^{\ast}]_{s}{\mathcal{V}}_{s-}\\\
&+\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}[\operatorname{vec}([L,L^{\ast}]),L^{\ast}]_{s}f_{s-}^{\ast}+\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}[\operatorname{vec}([L,L^{\ast}]),\operatorname{vec}([L,L^{\ast}])^{\ast}]{\mathcal{V}}_{s-}\\\
&:=I_{1}+I_{2}+I_{3}+I_{4}.\end{split}$ (A.24)

By Lemma 3.2 with $p=2$ and similar calculations as in (A.2) it follows that
${\mathbb{E}}\|{\mathcal{V}}_{s-}\|\|f_{s-}\|<\infty$ for all $s>0$ and the
map $s\mapsto{\mathbb{E}}\|{\mathcal{V}}_{s-}\|\|f_{s-}\|$ is locally bounded.
Thus, it follows from $(a.3)$ that we have
${\mathbb{E}}I_{2}={\mathbb{E}}I_{3}=0$. Now, Lemma 3.2 gives
${\mathbb{E}}\|{\mathcal{V}}_{s-}\|^{2}<\infty$ for all $s>0$ and local
boundedness of the map $s\mapsto{\mathbb{E}}\|{\mathcal{V}}_{s-}\|^{2}$. Using
the second-order stationarity of $(V_{s})_{s\in{\mathbb{R}}_{+}}$ in $(c.1)$,
the compensation formula and the formulas at p. 108 in Stelzer (2010)

$\begin{split}{\mathbb{E}}I_{4}&={\mathbb{E}}\bigg{(}\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}[\operatorname{vec}([L,L^{\ast}]),\operatorname{vec}([L,L^{\ast}])^{\ast}]{\mathcal{V}}_{s-}\bigg{)}\\\
&={\mathbb{E}}\bigg{(}\int_{0}^{\Delta}{\mathcal{V}}_{s-}{\rm
d}[\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}),(\operatorname{vec}([L,L^{\ast}]^{\mathfrak{d}}))^{\ast}]^{\mathfrak{d}}{\mathcal{V}}_{s-}\bigg{)}\\\
&=\int_{0}^{\Delta}{\mathbb{E}}\big{(}{\mathcal{V}}_{s-}\rho_{L}[I_{d^{2}}+K_{d}+\operatorname{vec}(I_{d})\operatorname{vec}(I_{d})^{\ast}]{\mathcal{V}}_{s-}\big{)}{\rm
d}s\\\
&=\rho_{L}\int_{0}^{\Delta}({\bm{Q}}+K_{d}{\bm{Q}}+I_{d^{2}}){\mathbb{E}}(\operatorname{vec}(V_{s})\operatorname{vec}(V_{s})^{\ast}){\rm
d}s\\\
&=\Delta\rho_{L}({\bm{Q}}+K_{d}{\bm{Q}}+I_{d^{2}}){\mathbb{E}}\operatorname{vec}(V_{0})\operatorname{vec}(V_{0})^{\ast}.\end{split}$
(A.25)

To compute ${\mathbb{E}}I_{1}$ we will need the following matrix identity,
which is based on Fact 7.4.30 (xiv) in Bernstein (2009). Let $A\in
M_{d,1}({\mathbb{R}})$ and $B,B^{2}\in M_{d,d}({\mathbb{R}})$ be symmetric
matrices. Then,

$\begin{split}\quad\quad&(A\otimes B+B\otimes A)(A\otimes B+B\otimes
A)^{\ast}=(A\otimes B+K_{d}(A\otimes B))(A\otimes B+K_{d}(A\otimes
B))^{\ast}\\\ &=(I+K_{d})(A\otimes B)(A^{\ast}\otimes
B)(I+K_{d})=(I+K_{d}){\bm{Q}}\operatorname{vec}(AA^{\ast})\operatorname{vec}(B^{2})(I+K_{d}).\end{split}$
(A.26)

Write
$b_{s}:={\mathbb{E}}\operatorname{vec}(G_{s}G_{s}^{\ast})\operatorname{vec}(V_{s})^{\ast}$,
which is finite by Lemma 3.2 with $p=2$. Using the compensation formula,
(A.26) and the definition of $f_{s}$ gives

$\begin{split}\quad\,\,&{\mathbb{E}}\bigg{(}\int_{0}^{\Delta}f_{s-}{\rm
d}[L,L^{\ast}]_{s}f_{s-}^{\ast}\bigg{)}\\\
&=(\sigma_{L}+\sigma_{W})\int_{0}^{\Delta}{\mathbb{E}}(f_{s}f_{s}^{\ast}){\rm
d}s\\\ &=(\sigma_{L}+\sigma_{W})\int_{0}^{\Delta}{\mathbb{E}}(G_{s-}\otimes
V_{s-}^{1/2}+V_{s-}^{1/2}\otimes G_{s-})(G^{\ast}_{s-}\otimes
V_{s-}^{1/2}+V_{s-}^{1/2}\otimes G_{s-}^{\ast}){\rm d}s\\\
&=(\sigma_{L}+\sigma_{W})\int_{0}^{\Delta}(I+K_{d}){\bm{Q}}b_{s}(I+K_{d}){\rm
d}s\\\
&=(\sigma_{L}+\sigma_{W})(I+K_{d}){\bm{Q}}\bigg{(}\int_{0}^{\Delta}b_{s}{\rm
d}s\bigg{)}(I+K_{d}).\end{split}$ (A.27)

Finally, it follows from (A.23) that

$\int_{0}^{\Delta}b_{s}^{\ast}{\rm
d}s=\int_{0}^{\Delta}{\mathbb{E}}\operatorname{vec}(V_{s})a_{s-}^{\ast}{\rm
d}s=(\sigma_{L}+\sigma_{W})^{-1}D.$ (A.28)

The result now is a direct consequence of (A.17), (A.23), (A.24), (A.25),
(A.27) and (A.28).

######  Remark A.2.

An inspection of the proofs of Lemmas 4.2 and A.1 shows that the moment
assumptions $(a.6)$ and $(c.2)$ are only needed to compute expectations of
stochastic integrals with the integrator $L$. If $L$ has paths of finite
variation, these expectations can be computed by using the compensation
formulas given in Facts A.1 and A.2 without $(a.6)$ and $(c.2)$.

## References

  * Applebaum (2009) David Applebaum.  _Lévy Processes and Stochastic Calculus_.  Cambridge University Press, Cambridge, 2009.
  * Bayracı and Ünal (2014) Selçuk Bayracı and Gazanfer Ünal.  Stochastic interest rate volatility modeling with a continuous-time GARCH(1,1) model.  _Journal of Computational and Applied Mathematics_ , 259:464–473, 2014.
  * Behme (2012) Anita Behme.  Moments of MGOU processes and positive semidefinite matrix processes.  _Journal of Multivariate Analysis_ , 111:183–197, 2012\.
  * Bernstein (2009) Dennis S. Bernstein.  _Matrix Mathematics: Theory, Facts, and Formulas_.  Princeton University Press, Princeton, 2009.
  * Bibbona and Negri (2015) Enrico Bibbona and Ilia Negri.  Higher moments and prediction-based estimation for the COGARCH(1,1) model.  _Scandinavian Journal of Statistics_ , 42(4):891–910, 2015.
  * Billingsley (2008) Patrick Billingsley.  _Probability and Measure_.  Wiley, New York, 2008.
  * Bollerslev (1986) Tim Bollerslev.  Generalized autoregressive conditional heteroskedasticity.  _Journal of Econometrics_ , 31(3):307–327, 1986\.
  * Bradley (2007) Richard Bradley.  _Introduction to Strong Mixing Conditions_ , volume 1.  Kendrick Press, Heber City, 2007.
  * Cohen and Rosinski (2007) Serge Cohen and Jan Rosinski.  Gaussian approximation of multivariate lévy processes with applications to simulation of tempered stable processes.  _Bernoulli_ , 13(1):195–210, 2007.
  * Cont (2001) Rama Cont.  Empirical properties of asset returns: stylized facts and statistical issues.  _Quantitative Finance_ , 1:223–236, 2001.
  * Do Rêgo Sousa (2019) Thiago Do Rêgo Sousa.  _Simulation-based estimation of time series and stochastic volatility processes_.  PhD thesis, Department of Mathematics, Technical University of Munich, 2019.  URL https://mediatum.ub.tum.de/?id=1518760.
  * do Rêgo Sousa et al. (2019) Thiago do Rêgo Sousa, Stephan Haug, and Claudia Klüppelberg.  Indirect inference for Lévy-driven continuous-time GARCH models.  _Scandinavian Journal of Statistics_ , (46):765–801, 2019.
  * Engle (1982) Robert F. Engle.  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  _Econometrica_ , 50(4):987–1007, 1982.
  * Engle and Kroner (1995) Robert F. Engle and Kenneth F. Kroner.  Multivariate simultaneous generalized ARCH.  _Econometric Theory_ , 11(1):122–150, 1995.
  * Fasen (2010) Vicky Fasen.  Asymptotic results for sample autocovariance functions and extremes of integrated generalized ornstein–uhlenbeck processes.  _Bernoulli_ , 16(1):51–79, 2010.
  * Fasen et al. (2006) Vicky Fasen, Claudia Klüppelberg, and Alexander Lindner.  Extremal behavior of stochastic volatility models.  In _Stochastic finance_ , pages 107–155. Springer, 2006.
  * Francq and Zakoïan (2019) Christian Francq and Jean-Michel Zakoïan.  _GARCH Models: Structure, Statistical Inference and Financial Applications_.  Wiley, Hoboken, NJ, second edition, 2019.
  * Guillaume et al. (1997) Dominique M Guillaume, Michel M Dacorogna, Rakhal R Davé, Ulrich A Müller, Richard B Olsen, and Olivier V Pictet.  From the bird’s eye to the microscope: A survey of new stylized facts of the intra-daily foreign exchange markets.  _Finance and stochastics_ , 1(2):95–129, 1997\.
  * Haug et al. (2007) Stephan Haug, Claudia Klüppelberg, Alexander Lindner, and Matthias Zapp.  Method of moment estimation in the COGARCH(1,1) model.  _Econometrics Journal_ , 10(2):320–341, 2007\.
  * Horn and Johnson (1991) Roger A Horn and Charles R Johnson.  _Topics in Matrix Analysis_.  Cambridge University Press, Cambridge, 1991.
  * Ibragimov and Linnik (1971) Ildar A. Ibragimov and Yuri V. Linnik.  _Independent and Stationary Sequences of Random Variables_.  Wolters-Noordhoff, Groningen, 1971.
  * Klenke (2013) Achim Klenke.  _Probability Theory: a Comprehensive Course_.  Springer, Berlin, 2013.
  * Klüppelberg et al. (2011) Claudia Klüppelberg, Ross Maller, and Alexander Szimayer.  The COGARCH: a review, with news on option pricing and statistical inference.  In J. Blath, P. Imkeller, and S. Roelly, editors, _Survey in Stochastic Processes. Proc. of the 33rd SPA Conference in Berlin_ , pages 29–50. EMS Series of Congress Reports, EMS Publishing House, Zürich, 2011\.
  * Krengel (1985) Ulrich Krengel.  _Ergodic Theorems_.  Walter de Gruyter, Berlin, 1985.
  * Kufner et al. (1977) Alois Kufner, Oldrich John, and Svatopluk Fuc̆ik.  _Function spaces_.  Noordhoff, Leyden, 1977.
  * Lindner and Maller (2005) Alexander Lindner and Ross Maller.  Lévy integrals and the stationarity of generalised Ornstein-Uhlenbeck processes.  _Stochastic Processes and their Applications_ , 115(10):1701–1722, 2005.
  * Magnus and Neudecker (1979) Jan R Magnus and Heinz Neudecker.  The commutation matrix: some properties and applications.  _The Annals of Statistics_ , 7:381–394, 1979.
  * Maller et al. (2008) Ross A Maller, Gernot Müller, and Alex Szimayer.  GARCH modelling in continuous time for irregularly spaced time series data.  _Bernoulli_ , 14(2):519–542, 2008.
  * Mátyás (1999) László Mátyás.  _Generalized Method of Moments Estimation_.  Cambridge University Press, Cambridge, 1999.
  * Müller (2010) Gernot Müller.  MCMC estimation of the COGARCH(1,1) model.  _Journal of Financial Econometrics_ , 8(4):481–510, 2010.
  * Pigorsch and Stelzer (2009a) Christian Pigorsch and Robert Stelzer.  On the definition, stationary distribution and second order structure of positive semi-definite Ornstein-Uhlenbeck type processes.  _Bernoulli_ , 15:754–773, 2009a.
  * Pigorsch and Stelzer (2009b) Christian Pigorsch and Robert Stelzer.  A multivariate Ornstein-Uhlenbeck type stochastic volatility model.  _Technical Report: www.uni-ulm.de/en/mawi/finmath/people/ stelzer/publications.html_ , 2009b.
  * Protter (2005) Philip Protter.  _Stochastic Integration and Differential Equations_.  Springer, New York, 2 edition, 2005.
  * R Core Team (2017) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2017.
  * Sato (1999) Ken-iti Sato.  _Lévy processes and Infinitely Divisible Distributions_.  Cambridge university press, Cambridge, 1999.
  * Schlemm and Stelzer (2012) Eckhard Schlemm and Robert Stelzer.  Quasi maximum likelihood estimation for strongly mixing state space models and multivariate Levy-driven CARMA processes.  _Electronic Journal of Statistics_ , 6:2185–2234, 2012\.
  * Stelzer (2010) Robert Stelzer.  Multivariate COGARCH(1,1) processes.  _Bernoulli_ , 16(1):80–115, 2010.
  * Stelzer and Vestweber (2019) Robert Stelzer and Johanna Vestweber.  Geometric ergodicity of the multivariate COGARCH(1,1) process.  _arXiv:1701.07859_ , 2019.
  * Stelzer et al. (2015) Robert Stelzer, Thomas Tosstorff, and Marc Wittlinger.  Moment based estimation of supou processes and a related stochastic volatility model.  _Statistics & Risk Modeling_, 32(1):1–24, 2015\.
  * Wilcox (1967) Ralph M. Wilcox.  Exponential operators and parameter differentiation in quantum physics.  _Journal of Mathematical Physics_ , 8(4):962–982, 1967.