††thanks: The PhD grant of M. D. Eggen is funded by NORSAR.

# The multivariate ARMA/CARMA transformation relation

Mari Dahl Eggen<EMAIL_ADDRESS>Department of Mathematics, University of
Oslo, P.O. Box 1053 Blindern, 0316 Oslo, Norway.

###### Abstract

A transformation relation between multivariate ARMA and CARMA processes is
derived through a discretization procedure. This gives a direct relationship
between the discrete time and continuous time analogues, serving as the basis
for an estimation method for multivariate CARMA models. We will see that the
autoregressive coefficients, making up the deterministic part of a
multivariate CARMA model, are entirely given by the transformation relation.
An Euler discretization convergence rate of jump diffusions is found for the
case of small jumps of infinite variation. This substantiates applying the
transformation relation for estimation of multivariate CARMA models driven by
NIG-Lévy processes. A two-dimensional CAR model is fit to stratospheric
temperature and wind data, as an example of how to apply the transformation
relation in estimation methods.

Discretization, model estimation, MCARMA, transformation relation, VARMA

††preprint: AIP/123-QED

##  I Introduction

The widely known family of autoregressive moving average (ARMA) processes
holds properties suitable for discrete time series modelling (see, e.g.,
[11]). Another family of processes that has become important in representing
time series is the family of continuous time ARMA (CARMA) processes, the
continuous time analogue to ARMA processes. Continuous time models allow for
irregularly spaced time series and provide the opportunity to derive explicit
formulas describing events and properties relying on dynamical systems.
Univariate CARMA processes and generalizations are used to model, for example,
variables in finance and energy markets, weather variables and turbulence
(see, e.g., [43], [21], [12], [5]). For an extensive overview of developments
and applications of CARMA processes, see [10] and references therein.

Most time series arising from dynamical systems in fields like natural
sciences, finance and economics will be more accurately represented by
multidimensional models describing dependencies between two or more variables
within a system. Discrete time multivariate ARMA processes (also called vector
ARMA (VARMA) processes), see [11], are used as models for this purpose, see,
e.g. [24] and [44] for practical examples. As argued above, modelling
continuous time series is often useful, with no exception in the
multidimensional case. A particularly useful property of a multivariate
continuous time version of VARMA models in applications, would be the
possibility of deriving explicit formulas for the crosscorrelation matrix
between modelled time series. The first derivation and proper representation
of multivariate CARMA (MCARMA) processes is derived in [33]. They are shown to
be the obvious continuous analogue of VARMA processes.

As a natural extension of the work in [33], where proper analogy between VARMA
and MCARMA processes is concluded, [39] derive fundamental results for
developing estimation theory for non-Gaussian MCARMA processes based on
(discrete) equidistant observations. As stated in [14], estimating a CARMA
model consists of three tasks: 1) choosing suitable integer values $p$ and
$q$, respectively describing the autoregressive and moving average orders; 2)
estimating autoregressive and moving average coefficients; 3) suggesting an
appropriate stochastic process to drive the model. These tasks hold for MCARMA
models as well. Statistical tests for choosing $p$ and $q$ already exist for
VARMA processes when normally distributed error terms are assumed, see e.g.
[24], however, a best consensus method does not seem to exist. Results in [39]
are utilized in [40] to derive an estimation procedure for autoregressive and
moving average model coefficients in the multivariate case. Some further
efforts have been made to enrich model coefficient estimation theory for
MCARMA processes, see, e.g., [19] and [20]. Finally, as stated in [14], there
are several well-established methods to estimate an appropriate driving Lévy
process. The work in [14] focuses on a parametric estimation method of
discretely observed MCARMA processes.

Let $\bm{Y}(t)$ represent the unique causal stationary solution, see [39], of
an MCARMA process, and define its discrete time $h$-sampled analogue as
$\bm{Y}^{h}(t)\triangleq\\{\bm{Y}(nh)\\}_{n\in\mathbb{Z}}$. As shown in [20],
discretely sampled MCARMA processes admits a VARMA representation with a weak
white noise. In [13], determining the coefficients of $\bm{Y}(t)$ from those
of $\bm{Y}^{h}(t)$ and vice versa is referred to as the sampling problem and
the embedding problem respectively. In particular, [20], [16] and [42] focus
on the sampling problem, where discrete time representations of sampled MCARMA
processes are estimated. To the best of our knowledge, the embedding problem
for MCARMA processes is not yet considered from the view of discretization
methods. However, discretization of CAR processes leading to a transformation
relation between AR and CAR processes is derived in [6]. Work leading up to
such discretization transformations is found in, e.g., [35] and [41].

The aim of this study is to derive a transformation relation between VARMA and
MCARMA processes through discretization of the MCARMA state space
representation. Define a Lévy-driven MCRAMA process as the solution,
$\bm{Y}(t)$, of the state space model

$\displaystyle d\bm{X}(t)=A\bm{X}(t)dt+\beta
d\bm{L}(t),\quad\bm{Y}(t)=C\bm{X}(t),\quad t\in\mathbb{R},$

see Section II.1 for a formal definition. The multidimensional Lévy-driven
stochastic differential equation (SDE) of Ornstein–Uhlenbeck type,
representing a system of one-dimensional SDEs, is discretized using the Euler
scheme. The discretized system of SDEs is further solved in a recursive manner
to estimate the state vector process $\bm{X}(t)$. The solution admits a VARMA
representation, giving a direct transformation relation between VARMA and
MCARMA processes.

Inspired by [2] and [7], convergence rates of (multidimensional) jump
diffusions with jumps of finite variance and infinite variations are studied.
The final result gives a convergence rate of these processes’ Euler scheme,
see, e.g., [36]. In particular, it is shown that the rate of convergence is
controllable by adjusting the discretization step size, as usual, and by the
size of the approximated small jumps. The state space representation of the
Lévy-driven MCARMA process is a special case of these jumps diffusions,
meaning that the convergence results hold for such models when driven by Lévy
processes of finite variance and infinite variation.

VARMA models have been used extensively in modelling of economic and financial
variables, see, e.g., [27] and [24]. In modelling and prediction of climate
and weather variables, VARMA models have been utilized to some degree, e.g.,
[34], [32] and [15]. The continuous time counterpart would be useful in
applications as well, see [18]. The VAR model is a simplified version of the
VARMA model, with the MCAR process as its continuous counterpart. The
multivariate ARMA/CARMA transformation relation is considerably simplified in
the VAR/MCAR case, and is easy to use in model estimation of MCAR processes.
To demonstrate how the transformation relation is used as a tool in model
estimation, a case study is performed where a two-dimensional MCAR model is
fit to weather data.

The atmospheric layer closest to the surface of the Earth is called the
troposphere. Above the troposphere lays the stratosphere, reaching up to about
$50$ km above the surface of the Earth. These two atmospheric layers are said
to interact through stratosphere-troposphere coupling, as weather conditions
in the stratosphere affect weather conditions in the troposphere. As explained
in, e.g., [28] and [38], probing and proper representation of the
stratosphere, combined with a thorough understanding of stratosphere-
troposphere coupling, has the potential to enhance long-term numerical surface
weather prediction. In this regard, stochastic modelling of stratospheric
weather dynamics will enlarge the ensemble of existing methods. Inspired by
prior work, e.g., [6] and [18], an MCAR process is successfully fit to
stratospheric temperature and wind data, giving a representation of the two-
dimensional dynamical system of dependent variables.

The structure of the paper is as follows. Section II defines the VARMA and
MCARMA representations, and introduces useful notation. In Section III, the
multivariate ARMA/CARMA transformation relation is found through
discretization of the MCARMA state space representation. Further, convergence
rates substantiating the transformation relation are derived. A case study
demonstrating how to use the transformation relation in MCARMA model
estimation is performed in Section IV.

Notation

Assume that a complete filtered probability space
$(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\}_{t\geq 0},P)$ is given as usual,
and that all stochastic (vector) processes
$\bm{X}(t)\triangleq\\{\bm{X}(t)\\}_{t\geq 0}$ are defined on that space.

Block-matrices are denoted as $M\in\mathbb{R}^{m\times n}$, having associated
matrices $M_{j}$ with matrix elements $\mu_{i}^{(j)}$. Vectors are denoted as
$\bm{V}\in\mathbb{R}^{m}$, with associated sub-vectors $\bm{V}_{j}$ and
elements $V_{i}$. The $d$-dimensional identity and zero matrices are given by
$\mathbb{1}_{d},\mathbb{0}_{d}\in\mathbb{R}^{d\times d}$. Denote by
$I_{B}(\cdot)$ the indicator function of some set $B$, and let $\det(\cdot)$
denote the determinant of matrices. Further, we work under the norm
$\|\cdot\|_{2}=(E[\lvert\cdot\rvert^{2}])^{1/2}$ throughout.

For convenience, let us recall the general representation of Lévy processes
through Itô-Lévy decomposition. Let $\bm{L}(t)$ be a Lévy process with values
in $\mathbb{R}^{m}$. Then, for each $t\geq 0$, there exist constant-valued
functions $\bm{\alpha}\in\mathbb{R}^{m}$ and $\sigma\in\mathbb{R}^{m\times
n}$, referred to as drift and diffusion respectively, such that

$\displaystyle\bm{L}(t)=\bm{\alpha}t+\sigma
W(t)+\int_{\lvert\bm{z}\rvert<1}\bm{z}\tilde{N}(t,d\bm{z})+\int_{\lvert\bm{z}\rvert\geq
1}\bm{z}N(t,d\bm{z}),$ (1)

where $\tilde{N}(dt,d\bm{z})=N(dt,d\bm{z})-\nu(d\bm{z})dt$,
$\bm{z}\in\mathbb{R}^{m}$, $\nu(U)$ is a Borel measure on
$\mathbb{R}^{m}\backslash\\{0\\}$, $\bm{W}(t)$ is a standard Brownian motion
process in $\mathbb{R}^{n}$ and $N(t,U)$ is a Poisson random measure on
$\mathbb{R}^{+}\times(\mathbb{R}^{m}\backslash\\{0\\})$. Further, $\bm{W}(t)$
and $N(t,U)$ are independent, and $\nu(U)$ is a Lévy measure, meaning
$\int_{\mathbb{R}^{m}\backslash\\{0\\}}(\lvert\bm{z}\rvert^{2}\wedge
1)\nu(d\bm{z})<\infty$ is satisfied. We will further assume that
$\int_{\lvert\bm{z}\rvert\geq\epsilon}\lvert\bm{z}\rvert^{2}\nu(d\bm{z})<\infty$
for some $0<\epsilon\leq 1$, such that the Lévy process has finite second
moments. See, e.g., [1] for a thorough introduction of Lévy processes.

##  II Multivariate ARMA and CARMA models

This section introduces the modelling framework of discrete time VARMA and
continuous time MCARMA processes. With the intention of deriving a
transformation relation between VARMA and MCARMA processes in Section III, the
the MCARMA model structure is studied more explicitly, and a recursive
parameter connected to the notation of the defined MCARMA framework is
defined.

###  II.1 Model representations

A short introduction to VARMA and MCARMA processes is presented in the
following. The notation and definitions are inspired by, [24], [33] and [39].

Define the autoregressive and moving average matrix polynomials as

$\displaystyle\phi(z)=\left(I-\phi_{1}z-\phi_{2}z^{2}-\cdots-\phi_{p}z^{p}\right)\quad\text{and}\quad\theta(z)=\left(I+\theta_{1}z+\theta_{2}z^{2}+\cdots+\theta_{q}z^{q}\right),$
(2)

respectively, and a backshift operator S as $S\bm{x}(t)=\bm{x}(t-1)$. Then,
the VARMA process might be stated as

$\displaystyle\phi(S)\bm{x}(t)=\theta(S)\bm{\mathcal{E}}(t),$ (3)

where $\bm{x}(t)\in\mathbb{R}^{d}$ is a sequence of random vectors and
$\bm{\mathcal{E}}(t)\in\mathbb{R}^{m}$ is a sequence of serially uncorrelated
i.i.d. random vectors with zero mean and common covariance matrix $\Sigma$.
Stationarity and invertibility of the model is secured as long as all roots of
$\det(\phi(z))$ and $\det(\theta(z))$ are outside the unit circle,
respectively, see [31] and [23]. As seen in Eq. (2) and (3), the VARMA model,
sometimes denoted as VARMA($p,q$), is determined by an autoregressive
parameter, $p\in\mathbb{N}$, a moving average parameter, $q\in\mathbb{N}$, and
the dimensionalities $d,m\in\mathbb{N}$.

In [33], the need for MCARMA processes were addressed, and further derived,
such that the joint behaviour of $d$ different time series can be modelled
continuously. As far as we know, this is the first mentioning of MCARMA
processes in literature. We adapt the notation from [33] and [39], and define
the MCARMA process as follows.

Let $\bm{L}(t)$ be a Lévy process with values in $\mathbb{R}^{m}$, as defined
in Eq. (1). Then the $d$-dimensional Lévy-driven MCARMA process, $\bm{Y}(t)$,
is given by the state space representation

$\displaystyle d\bm{X}(t)=A\bm{X}(t)dt+\beta
d\bm{L}(t),\quad\bm{Y}(t)=C\bm{X}(t).$ (4)

Here, $\bm{X}(t)$ is required to be a unique stationary solution of the given
Ornstein–Uhlenbeck type state space model, with $A$, $\beta$ and $C$ given as

$\displaystyle
A=\begin{bmatrix}\mathbb{0}_{d}&\mathbb{1}_{d}&\mathbb{0}_{d}&\cdots&\mathbb{0}_{d}\\\
\mathbb{0}_{d}&\mathbb{0}_{d}&\mathbb{1}_{d}&\cdots&\mathbb{0}_{d}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
\mathbb{0}_{d}&\mathbb{0}_{d}&\mathbb{0}_{d}&\cdots&\mathbb{1}_{d}\\\
-A_{p}&-A_{p-1}&-A_{p-2}&\cdots&-A_{1}\\\ \end{bmatrix}\in\mathbb{R}^{pd\times
pd},\quad\beta=\left(\beta_{1}^{T}\cdots\beta_{p}^{T}\right)^{T}\in\mathbb{R}^{pd\times
m},$
$\displaystyle\beta_{p-\kappa}=-I_{\\{0,\ldots,q\\}}(\kappa)\left[\sum_{i=1}^{p-\kappa-1}A_{i}\beta_{p-\kappa-i}-B_{q-\kappa}\right]\quad\text{and}\quad
C=\left(\mathbb{1}_{d},\mathbb{0}_{d},\ldots,\mathbb{0}_{d}\right)\in\mathbb{R}^{d\times
pd},$

where $\mathbb{1}_{d},\mathbb{0}_{d}\in\mathbb{R}^{d\times d}$ represents the
$d$-dimensional identity and zero matrices respectively, and
$I_{\\{0,\ldots,q\\}}(\cdot)$ is the indicator function of the set
$\\{0,\ldots,q\\}$. Further, we have that $A_{j}\in\mathbb{R}^{d\times d}$,
and $B_{q-\kappa}\in\mathbb{R}^{d\times m}$ for $j\in\\{1,\ldots,p\\}$,
$\kappa\in\\{0,\ldots,q\\}$ and positive integers $p>q$. Notice also that
$\beta_{p-\kappa}\in\mathbb{R}^{d\times m}$ for $0\leq\kappa\leq p-1$, which
is a zero matrix if $\kappa\notin\\{0,\ldots,q\\}$. As long as the driving
process, $\bm{L}(t)$, admits finite variance, the MCARMA process is assured to
have a unique causal stationary solution when the eigenvalues of $A$ have
strictly negative real parts.

###  II.2 The MCARMA model structure

The aim of this work is to derive a transformation relation between discrete
time VARMA processes and continuous time MCARMA processes through
discretization of the MCARMA representation in Eq. (4). To prepare for this
the state space model, $\bm{X}(t)$, is stated and inspected on a more explicit
level. That is, the block-matrix $A$ is written explicitly with conveniently
defined notation, and each matrix $A_{j}$ is structured into collections. The
defined notation is used to construct a recursive parameter in Section II.3,
which is further utilized to solve the discretized system of SDEs making up
the MCARMA model. The mentioned collection-structure is made to ease
readability and understanding throughout this process.

As indicated in Section II.1, the state vector process $\bm{X}(t)$ and Lévy
vector process $\bm{L}(t)$ takes values in $\mathbb{R}^{pd}$ and
$\mathbb{R}^{m}$ respectively. Denote each element of these vector processes
as $X_{i_{X}}(t)$, $i_{X}\in\\{1,\ldots,pd\\}$, and $L_{i_{L}}(t)$,
$i_{L}\in\\{1,\ldots,m\\}$, and further denote each $d$-dimensional sub-vector
of $\bm{X}(t)$ as $\bm{X}_{l}(t)\in\mathbb{R}^{d}$, $l\in\\{1,\ldots,p\\}$.
Then the state space model in Eq. (4) might be considered as a system of $p$
equations blocks

$\displaystyle
d\bm{X}_{l}(t)=\mathbb{1}_{d}\bm{X}_{l+1}(t)dt+\beta_{p-\kappa}d\bm{L}(t),$
(5) $\displaystyle d\bm{X}_{p}(t)=[-A_{p}\cdots-
A_{1}]\bm{X}(t)dt+\beta_{p}d\bm{L}(t),$ (6)

with $(l,\kappa)\in\\{(1,p-1),(2,p-2),\ldots,(p-1,1)\\}$. Note that the
equation block in Eq. (6) corresponds to $l=p$ and $\kappa=0$. The index-
dependence between the block number index, $l$, and the moving average index,
$\kappa$, is summarised in Table 1, together with the corresponding index
number of each one-dimensional SDE in the system.

Further, group the $p$ equation blocks into three disjoint collections. That
is, define the solution-collection $\mathcal{C}^{S}$ as holding SDEs in block
$1$ ($l=1$), the recursive-collection $\mathcal{C}^{R}$ as holding SDEs in
blocks $2$ to $p-1$ ($l=2,\ldots,p-1$) and the AR-collection
$\mathcal{C}^{AR}$ as holding SDEs in block $p$. The collection-structure is
presented in Table 1. A $d$-dimensional dynamical system known to follow an
MCARMA process, $\bm{Y}(t)$, is given by the vector process $\bm{X}_{1}(t)$
(see the model setup in Eq. (4) to understand why). This is why
$\mathcal{C}^{S}$ is called the solution-collection. The time lag dependence
structure of the dynamical system is entirely described through the SDEs in
Eq. (6), giving the name to the AR-collection, $\mathcal{C}^{AR}$. Finally,
the system of one-dimensional SDEs making up the MCARMA model has to be solved
recursively to obtain $\bm{X}_{1}(t)$, with the SDEs in $\mathcal{C}^{AR}$ as
a starting point. The SDEs in the recursive-collection, $\mathcal{C}^{R}$, are
used for this purpose. Notice that an increasing number of lags, $p$, in the
MCARMA model requires a larger collection $\mathcal{C}^{R}$ to make the system
of equations solvable.

Table 1: Dependence structure between collections $\mathcal{C}^{S}$, $\mathcal{C}^{R}$ and $\mathcal{C}^{AR}$, block number index, $l$, moving average index, $\kappa$, and SDE index numbers.  Collection |  Block number index, $l$ |  Moving average index, $\kappa$ | SDE index number
---|---|---|---
$\mathcal{C}^{S}$ | 1 | p-1 | $\\{1,\ldots,d\\}$
$\mathcal{C}^{R}$ | 2 | p-2 | $\\{d+1,\ldots,pd-d\\}$
⋮ | ⋮
p-2 | 2
p-1 | 1
$\mathcal{C}^{AR}$ | p | 0 | $\\{pd-d+1,\ldots,pd\\}$

A general solution, $\bm{X}_{1}(t)$, of the state space model representing the
MCARMA process will be derived through the Euler scheme of $d\bm{X}(t)$. Note
that the derived solution will be an approximation with accuracy depending on
the incremental value defining the discretization scheme, see Section III. To
prepare for element-wise discretization, define each matrix $A_{j}$,
$j\in\\{1,\ldots,p\\}$, and $\beta_{p-\kappa}$, $\kappa\in\\{0,\ldots,p-1\\}$,
in the respective block-matrices $A$ and $\beta$, as

$\displaystyle
A_{j}=\begin{bmatrix}\alpha_{1}^{(j)}&\alpha_{2}^{(j)}&\cdots&\alpha_{d}^{(j)}\\\
\alpha_{d+1}^{(j)}&\cdots&\cdots&\alpha_{2d}^{(j)}\\\
\vdots&\vdots&\vdots&\vdots\\\
\alpha_{(d-1)d+1}^{(j)}&\cdots&\cdots&\alpha_{dd}^{(j)}\end{bmatrix}\quad\text{and}\quad\beta_{p-\kappa}=\begin{bmatrix}\beta_{1}^{(\kappa)}&\beta_{2}^{(\kappa)}&\cdots&\beta_{m}^{(\kappa)}\\\
\beta_{m+1}^{(\kappa)}&\cdots&\cdots&\beta_{2m}^{(\kappa)}\\\
\vdots&\vdots&\vdots&\vdots\\\
\beta_{(d-1)m+1}^{(\kappa)}&\cdots&\cdots&\beta_{dm}^{(\kappa)}\end{bmatrix}.$
(7)

This gives elements

$\displaystyle-\sum_{l=1}^{p}\sum_{s=1}^{d}\alpha_{(k-1)d+s}^{(p-l+1)}X_{(l-1)d+s}(t)dt\quad\text{and}\quad\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(\kappa)}dL_{r}(t),$
(8)

of the respective matrix-vector products $[-A_{p}\cdots-A_{1}]\bm{X}(t)dt$ and
$\beta_{p-\kappa}d\bm{L}(t)$, for all $\kappa$. Here, $\kappa$ depends on the
equation block, $l$, under consideration and $k\in\\{1,\ldots,d\\}$ is set
based on the dimension of interest. For example, to state the SDE of dimension
$2$ in equation block $p$, set $\kappa=0$ (confer with Table 1) and $k=2$.
Also, remember that the element-wise subscripts of $\bm{X}(t)$ and $\bm{L}(t)$
are defined as $i_{X}\in\\{1,\ldots,pd\\}$ and $i_{L}\in\\{1,\ldots,m\\}$
respectively, corresponding to the respective elements in Eq.(8).

###  II.3 The recursive parameter

The solution, $\bm{Y}(t)$, of the $d$-dimensional MCARMA($p,q$) process is
given by the state space model in Eq. (4). In this section, a recursive
parameter is defined to solve a discretized version of the state space model
recursively. We will see that the recursive parameter can be used to express
all the one-dimensional SDEs in Eq. (5).

It is important to realize that the $d$-dimensions of an MCARMA model each
represents a dynamical variable, each dependent on $p$ previous time steps,
also across dimensions, with a moving average of degree $q$. The $p$ times
lagged $d$ dimensions are represented through $d$ SDEs in each of the $p$
equation blocks introduced in Section II.2 (see Eq. (5) and (6)). That is, the
collections $\mathcal{C}^{S}$, $\mathcal{C}^{R}$ and $\mathcal{C}^{AR}$ holds
$d$, $(p-2)d$ and $d$ SDEs, respectively (see Table 1). Remember that the
$(p-2)d$ SDEs in $\mathcal{C}^{R}$ are distributed into $p-2$ equation blocks
(that is block $2$ to block $p-1$), where block $p-1$ represents the first
recursive step, block $p-2$ the second recursive step, and so on. To keep
track of all the one-dimensional SDEs in a recursive procedure from block $p$
to block $1$, a recursive parameter is introduced as

$\displaystyle Q_{i}^{(l)}\triangleq(Q_{i}^{(l)}\mid k)=(l-i)d+k,$ (9)

where $k\in\\{1,\ldots,d\\}$ and $l\in\\{1,\ldots,p\\}$, such that $1\leq
i\leq l$ for every $l$, and $i>1$ when $l=p$. Note that $d$ is the fixed
dimension parameter of the MCARMA process, $i$ and $l$ are dynamical indexes
holding track of equation blocks $1$ to $p-1$ during the recursive procedure,
and that $k$ is set depending on the dimension of interest.

To understand how the recursive parameter is utilized, consider the following.
Let the recursive parameter represent element index, $i_{X}$, of $\bm{X}(t)$,
meaning that elements of $\bm{X}(t)$ are given as $X_{Q_{i}^{(l)}}(t)$. By the
definition in Eq. (9), the recursive parameter takes values in
$\\{1,\ldots,pd-d\\}$, corresponding to SDEs in $\mathcal{C}^{S}$ and
$\mathcal{C}^{R}$ (see Table 1). Further, the recursive parameter holds two
properties that are utilized in the recursive procedure:

  1. Property 1.

SDEs in an equation block with lower block number index, $l$, might be
expressed in more ways using $X_{Q_{i}^{(l)}}(t)$ than SDEs in an equation
block with higher block number index. This property is useful because the
recursive procedure iterates through SDEs in an equation block with lower
block number index more times. That is, in the first recursive step SDEs from
block $p$ are substituted into SDEs in block $p-1$, in the second recursive
step these SDEs are further substituted into SDEs in block $p-2$. This
procedure continues until the last recursive step where the result is
substituted into SDEs in block $1$. Substituted SDEs from block $p-1$ into
SDEs in block $p-2$ (in the second recursive step), which now also contain
substituted SDEs from block $p$, will be substituted into one less block. A
similar argument holds for all further recursive steps starting in blocks
$p-2$ to $2$. See an illustrative example below of this recursive parameter
property for a case with $p=3$, where SDEs in equation block $2$ can be
expressed using $X_{Q_{i}^{(l)}}(t)$ in two different ways, and SDEs in
equation block $1$ can be expressed using $X_{Q_{i}^{(l)}}(t)$ in three
different ways.

  2. Property 2.

By definition (Eq. (9)), we have that
$Q_{1}^{(l)}=Q_{1}^{(l)}-d+d=Q_{2}^{(l)}+d$. Similarly,
$Q_{2}^{(l)}=Q_{3}^{(l)}+d$, and so on, until $Q_{p-1}^{(l)}=Q_{p}^{(p)}+d$.
This iterative property is used to conveniently express the SDEs in
collections $\mathcal{C}^{S}$ and $\mathcal{C}^{R}$ (see Proposition 1), and
further to recursively solve the discretized $d$-dimensional MCARMA process by
the recursive procedure from equation block corresponding to
$X_{Q_{1}^{(l)}}$, for all $l$, to the final solution block corresponding to
$X_{Q_{p}^{(p)}}(t)=X_{k}(t)\in\bm{X}_{1}(t)$.

The following proposition states how the recursive parameter might be used to
represent the one-dimensional SDEs in collections $\mathcal{C}^{S}$ and
$\mathcal{C}^{R}$.

######  Proposition 1.

Let the recursive parameter be as given in Eq. (9), where
$k\in\\{1,\ldots,d\\}$ and $l\in\\{1,\ldots,p\\}$. For $1\leq l\leq p$ and
$1\leq i\leq l$ ($i>1$ when $l=p$), corresponding to SDE numbers in
collections $\mathcal{C}^{S}$ and $\mathcal{C}^{R}$ (see Table 1), we have
that

$\displaystyle
X_{Q_{i}^{(l)}+d}(t)dt=dX_{Q_{i}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}dL_{r}(t).$
(10)

###### Proof.

See Appendix A. ∎

######  Example.

The MCARMA process dynamics with given parameters $p=4$, $d=2$ and $m=d$ is
given by

$\displaystyle\begin{cases}d\bm{X}_{1}(t)=\bm{X}_{2}(t)dt+\beta_{1}d\bm{L}(t)\\\
d\bm{X}_{2}(t)=\bm{X}_{3}(t)dt+\beta_{2}d\bm{L}(t)\\\
d\bm{X}_{3}(t)=\bm{X}_{4}(t)dt+\beta_{3}d\bm{L}(t)\\\
d\bm{X}_{4}(t)=[-A_{4}\cdots-
A_{1}]\bm{X}(t)dt+\beta_{4}d\bm{L}(t),\end{cases}$ (11)

where
$\bm{X}_{1}(t)=[X_{1}(t),X_{2}(t)],\ldots,\bm{X}_{4}(t)=[X_{7}(t),X_{8}(t)]$,
$\bm{X}(t)=[\bm{X}_{1}(t),\ldots,\bm{X}_{4}(t)]$,
$\bm{L}(t)=[L_{1}(t),L_{2}(t)]$, and $A_{j}$ and $\beta_{p-\kappa}$ are given
in Eq. (7). Notice that Eq. (11) is a system of eight SDEs in four equation
blocks. That is, one equation block in $\mathcal{C}^{S}$, two in
$\mathcal{C}^{R}$ and one in $\mathcal{C}^{AR}$ (see Table 1). Proposition 1
states that SDEs in collections $\mathcal{C}^{S}$ and $\mathcal{C}^{R}$
(corresponding to Eq. (5)) can be written in terms of the recursive parameter.
The following points illustrate how property 1 of the recursive parameter is
utilized to recursively solve the discretized $d$-dimensional MCARMA process.
Concentrating on dimension number $2$, corresponding to considering a
recursive parameter with $k=2$ (see Eq. (9)) we see that $dX_{2}(t)$ might be
expressed by the recursive parameter in three different ways:

  1. 1.

With equation block $3$ as a starting point in the backwards recursive
procedure, corresponding to equation number $Q_{1}^{(3)}$, $dX_{2}(t)$ is
obtained when we hit $Q_{3}^{(3)}$, after two recursive steps;

  2. 2.

With equation block $2$ as starting point in the backwards recursive
procedure, corresponding to equation number $Q_{1}^{(2)}$, $dX_{2}(t)$ is
obtained when we hit $Q_{2}^{(2)}$, after one recursive step;

  3. 3.

With equation block $1$ as starting point in the backwards recursive
procedure, $dX_{2}(t)$ is already obtained by $Q_{1}^{(1)}$.

Further, $dX_{4}(t)$ might be expressed by the recursive parameter in two
different ways:

  1. 1.

With equation block $3$ as a starting point in the backwards recursive
procedure, corresponding to equation number $Q_{1}^{(3)}$, $dX_{2}(t)$ is
obtained when we hit $Q_{2}^{(3)}$, after one recursive step;

  2. 2.

With equation block $2$ as starting point in the backwards recursive
procedure, $dX_{4}(t)$ is already obtained by $Q_{1}^{(2)}$.

Notice that the SDEs in equation block $2$ might be expressed with the
recursive parameter when starting the recursive procedure in block $3$ or $2$,
however, not when starting it in equation block $1$. SDEs in equation block
$1$ might be expressed with the recursive parameter regardless of the starting
point of the recursive procedure. This is because the recursive parameter is
defined to follow a backwards recursive pattern. All of the starting points,
blocks $p$ to $2$, of the backwards recursive procedure of a general MCARMA
process ends in block $1$ (SDEs in $\mathcal{C}^{S}$) with
$X_{Q^{(l)}_{l}}(t)$, where $l$ equals the block number index of the recursive
procedure starting point.

##  III The transformation relation

The Lévy-driven MCARMA process dynamics is discretized using an Euler scheme.
Through the discretized version a transformation relation between the
continuous time MCARMA process and the discrete time VARMA process is derived.
The convergence rate of the discretized MCARMA process dynamics is finally
assessed for a driving Lévy process with finite variance and infinite
variations, which is relevant for the case study in Section IV.

###  III.1 From the MCARMA Euler scheme to the transformation relation

In this section, an approximated solution, $\bm{x}_{1}(t)$, of the
$d$-dimensional MCARMA($p,q$) process with an $m$-dimensional driving Lévy
process is found through discretization of the state space model dynamics,
$d\bm{X}(t)$, in Eq. (4). A backwards recursive procedure is performed on the
discretized process to find the solution, which finally takes the form of a
VARMA process with coefficients given by the MCARMA coefficients. This is what
we refer to as the multivariate ARMA/CARMA transformation relation, see
Theorem 1.

Consider the equidistant time discretization

$\displaystyle t^{h}=\\{0=t_{0}<t_{1}<\dots<t_{N}=T\\},$ (12)

of a given time interval $[0,T]$ with fixed incremental value $(0,1]\ni
h=t_{i+1}-t_{i}$, for $i\in\\{0,\ldots N-1\\}$. That is, the discretized time
interval might be written as $[0,h,2h,\ldots,Nh]$, and a time step from an
arbitrary point in time, $t$, is represented by $t+h$. Assume that the law of
stationary and independent increments $\Delta L(t)\triangleq L(t+h)-L(t)$ is
known. Then, using the Euler scheme (see [29] and [37]), the SDEs in
$\mathcal{C}^{S}\cup\mathcal{C}^{R}$ (see Proposition 1, Eq. (10)) might be
written as the piecewise constant process

$\displaystyle
x_{Q_{i}^{(l)}+d}(t)=\frac{1}{h}\left(x_{Q_{i}^{(l)}}(t+h)-x_{Q_{i}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}\Delta
L_{r}(t)\right).$ (13)

Through a backwards recursive procedure of discretized SDEs on this form, an
approximated solution of the MCARMA process (Eq. (4)) is found as
$\bm{X}_{1}(t)\simeq\bm{x}_{1}(t)=[x_{(Q_{l}^{(l)}\mid
k=1)}(t),\ldots,x_{(Q_{l}^{(l)}\mid k=d)}(t)]=[x_{1}(t),\ldots,x_{d}(t)]$. The
backwards recursive procedure has $p-1$ starting points, namely equation
blocks $p$ to $2$ (corresponding to SDEs in
$\mathcal{C}^{R}\cup\mathcal{C}^{AR}$), where the SDE index number in each of
these blocks are given by $Q_{1}^{(l)}$, for $l=2,\ldots,p$. The final goal is
to express all SDEs represented by $Q_{1}^{(l)}$, in terms of $Q_{l}^{(l)}$.
See the example in Section II.3 to understand how the recursive procedure is
iterated. As seen in Eq. (5) and (6), the structures of the SDEs in
$\mathcal{C}^{R}$ and $\mathcal{C}^{AR}$ are different. In the following
lemma, SDEs represented by $Q_{1}^{(l)}$ is written in terms of $Q_{l}^{(l)}$,
for the the collection $\mathcal{C}^{R}$. This will further be modified in
terms of the SDEs in $\mathcal{C}^{AR}$ to find the final approximated
solution.

######  Lemma 1.

The approximated solution of SDEs in equation block $l\in\\{1,\ldots,p-1\\}$,
for $i=1$, corresponding to the starting point of each recursive step in the
derivation of the transformation relation in Theorem 1, is given by

$\displaystyle x_{Q_{1}^{(l)}+d}(t)=$
$\displaystyle\frac{1}{h^{l}}\sum_{n=0}^{l}(-1)^{n}b_{n}^{l}x_{k}(t+(l-n)h)$
$\displaystyle-\sum_{w=0}^{l-1}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+w)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h),$

where $p$ is the number of lags, $d$ is the total number of dimensions,
$k\in\\{1,\ldots,d\\}$ is the dimension of interest, $m$ is the number of
independent driving Lévy processes, $b_{n}^{i}$ is defined recursively from
Eq. (47) (Appendix B), and $h$ is the Euler discretization step size.

###### Proof.

See Appendix B. ∎

Lemma 1 gives an implicit formula for the approximated solution of SDEs in
$\mathcal{C}^{R}$. However, the autoregressive behaviour of the MCARMA process
is described by SDEs in $\mathcal{C}^{AR}$. This information has to be added
to the formula in Lemma 1, to find the final approximated solution of the
MCARMA process. Theorem 1 states a formula for the $p$ times lagged variable
of dimension $k\in\\{1,\ldots,d\\}$.

######  Theorem 1.

The multivariate ARMA/CARMA transformation relation is given by

$\displaystyle x_{k}(t+ph)=$ $\displaystyle
x_{k}(t+(p-1)h)-h^{p}\sum_{s=1}^{d}\alpha_{(k-1)d+s}^{(p)}x_{s}(t)+h^{p-1}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(0)}\Delta
L_{r}(t)$
$\displaystyle-\sum_{l=2}^{p}\sum_{s=1}^{d}\alpha_{(k-1)d+s}^{(p-l+1)}\Bigg{(}h^{p-l+1}\sum_{n=0}^{l-1}(-1)^{n}b_{n}^{l-1}x_{s}(t+(l-1-n)h)$
$\displaystyle\qquad\qquad\qquad\qquad-\sum_{w=0}^{l-2}h^{p-w-1}\sum_{r=1}^{m}\beta_{(s-1)m+r}^{(p-l+w+1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h)\Bigg{)}$
$\displaystyle-\sum_{n=1}^{p-1}(-1)^{n}b_{n}^{p-1}\left(x_{k}(t+(p-n)h)-x_{k}(t+(p-1-n)h)\right)$
$\displaystyle+\sum_{w=0}^{p-2}h^{p-w-2}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(w+1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\left(\Delta
L_{r}(t+(w-v+1)h)-\Delta L_{r}(t+(w-v)h)\right),$

where $p$ is the number of lags, $d$ is the total number of dimensions,
$k\in\\{1,\ldots,d\\}$ is the dimension of interest, $m$ is the number of
independent driving Lévy processes, $b_{n}^{i}$ is defined recursively from
Eq. (47) (Appendix B), and $h$ is the Euler discretization step size.

###### Proof.

See Appendix C. ∎

Notice that the multivariate ARMA/CARMA transformation relation expresses the
$p$ times lagged variable of dimension $k\in\\{1,\ldots,d\\}$ as a linear
combination of $0$ to $p-1$ lagged variables of all model dimensions and
driving Lévy processes. That is, the approximated solution of the MCARMA
process is represented by a VARMA process.

###  III.2 An analysis of convergence rates

A multivariate ARMA/CARMA transformation relation was derived in the previous
section through an Euler discretization of the Lévy-driven MCARMA process. For
any practical application, it is important that the Euler scheme converges. An
Euler scheme convergence rate for jump diffusions with jumps of finite
variance and infinite variations is derived in this section. In Section IV.1
we will see that the result holds for an NIG-Lévy-driven MCARMA process. Note
that vector processes and deterministic vector functions are written without
bold font in this section.

Consider the jump diffusion

$\displaystyle Z(t)=$ $\displaystyle
Z(0)+\int_{0}^{t}\tilde{a}(s,Z(s))ds+\int_{0}^{t}b(s,Z(s))dW(s)+\int_{0}^{t}\int_{\mathbb{R}^{m}\backslash\\{0\\}}\gamma(s-,Z(s-),z)\tilde{N}(ds,dz),$
(14)

where $\tilde{a}(t,x)=\left(a(t,x)+\int_{\lvert
z\rvert\geq\epsilon}\gamma(t,x,z)\nu(dz)\right)$. We assume that the usual
integrability conditions of coefficients are satisfied, and that $W(t)$ and
$N(t,U)$ are independent stochastic processes. Further, we assume that
$\gamma(t,x,z)=g(z)\eta(t,x)$, where $t\to\eta(t,x)$ is càdlàg, and that the
finite variance condition

$\displaystyle
G^{2}(\infty)=\int_{\mathbb{R}^{m}\backslash\\{0\\}}g^{2}(z)\nu(dz)<\infty,$

is satisfied.

Convergence rates of jump diffusions depend on the behaviour of the Lévy
measure at origin, see, e.g., [30]. A complete discussion on the Euler scheme
(among others) and convergence rates for Lévy processes with Lévy measure of
finite total mass, $\nu(U)<\infty$, is found in [36]. Note that this
corresponds to Lévy processes where the jump part is a compound Poisson
process. For Lévy processes with Lévy measure of infinite total mass,
$\nu(U)=\infty$, the question of Euler scheme convergence rates is more
intricate. In [2], the cases $\int_{\lvert z\rvert\leq 1}\lvert
z\rvert\nu(dz)<\infty$ and $\int_{\lvert z\rvert\leq 1}\lvert
z\rvert\nu(dz)=\infty$ are discussed. In the following, inspired by work in
[36], [2] and [7], an Euler scheme convergence rate for jump diffusions of
infinite variations, as given in Eq. (14), is derived for the case
$\nu(U)=\infty$.

In cases where the jump part of $Z(t)$ is of finite variation, convergence
rates might be obtained by replacing the small jumps by their expected value,
or by simply removing the small jumps. In the infinite variation case,
removing the small jumps would not be appropriate, as the small jumps dominate
in that case. Define the function

$\displaystyle G^{2}(\epsilon)\triangleq\int_{\lvert
z\rvert<\epsilon}g^{2}(z)\nu(dz),$ (15)

with natural property $\lim_{\epsilon\to 0}G^{2}(\epsilon)=0$. We propose, as
in [7], to approximate the small jumps part of Eq. (14) as

$\displaystyle\int_{0}^{t}\int_{\lvert
z\rvert<\epsilon}g(z)\eta(s-,Z(s-))\tilde{N}(ds,dz)\simeq\int_{0}^{t}G(\epsilon)\eta(s,Z_{\epsilon}(s))dB(s)\quad\leftrightarrow\quad
S(t)\simeq S_{\epsilon}(t),$ (16)

where $B(t)$ is a Brownian motion process independent of $W(t)$ and $N(t,U)$.
That is, $Z(t)$ in Eq. (14) is approximated by the process

$\displaystyle\begin{split}Z_{\epsilon}(t)=&Z(0)+\int_{0}^{t}\tilde{a}(s,Z_{\epsilon}(s))ds+\int_{0}^{t}b(s,Z_{\epsilon}(s))dW(s)+\int_{0}^{t}G(\epsilon)\eta(s,Z_{\epsilon}(s))dB(s)\\\
&\quad+\int_{0}^{t}\int_{\lvert
z\rvert\geq\epsilon}g(z)\eta(s-,Z_{\epsilon}(s-))\tilde{N}(ds,dz).\end{split}$
(17)

Notice that the jump part of $Z_{\epsilon}(t)$ is represented by a general
(compensated) compound Poisson process. Before stating the convergence rate
for the approximation $Z_{\epsilon}(t)$, we state an intermediate result. The
following assumption, inspired by [36], holds throughout this work.

######  Assumption 1.

Assume that the jump diffusions $Z(t)$ and $Z_{\epsilon}(t)$, in Eq. (14) and
(17) respectively, satisfy

$\displaystyle E(\lvert Z(0)\rvert^{2})=E(\lvert
Z_{\epsilon}(0)\rvert^{2})<\infty\quad\text{and}\quad E(\lvert
Z_{\epsilon}(0)-Z_{\epsilon}^{h}(0)\rvert^{2})\leq K_{0}h,$

where $Z_{\epsilon}^{h}(0)$ represents the initial condition of the Euler
discretization of $Z_{\epsilon}(t)$. Further, for $t\in[0,T]$ and
$x,y\in\mathbb{R}^{m}$, the drift and diffusion coefficients, as well as
$\eta(t,x)$ in the separable jump coefficient $\gamma(t,x,z)=g(z)\eta(t,x)$,
satisfy the Lipschitz and linear growth conditions

$\displaystyle\lvert f(t,x)-f(t,y)\rvert\leq K_{1}\lvert
x-y\rvert\quad\text{and}\quad\lvert f(t,x)\rvert^{2}\leq K_{2}(1+\lvert
x\rvert^{2}).$ (18)

Define $K=\max(K_{1},K_{2})$, which is used in the following proofs without
notice.

######  Lemma 2.

Let $S(t)$ and $S_{\epsilon}(t)$ be given as in Eq. (16), and denote the
process $S_{\epsilon}(t)$ by $S_{\epsilon,Z(t)}(t)$, when it is state-
dependent on $Z(t)$ rather than $Z_{\epsilon}(t)$ (see Eq. (14) and (17)). The
convergence rate of $S_{\epsilon,Z(t)}(t)$ is given by

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
S(t)-S_{\epsilon,Z(t)}(t)\rvert\|_{2}\leq C^{1/2}G(\epsilon),$

where $C=8KT(1+C_{1}e^{C_{1}})$, $C_{1}=2K(2T+2+G^{2}(\infty))/T$ and
$G(\epsilon)$ given in Eq. (15).

###### Proof.

Define $f\triangleq f(t,Z(t))$. First, we show a boundedness result of $Z(t)$.
By the triangle inequality, the identity $\lvert f+g\rvert^{p}\leq
2^{p-1}(\lvert f\rvert^{p}+\lvert g\rvert^{p})$, Cauchy-Schwarz inequality,
Itô isometry, monotone convergence and linear growth, we have

$\displaystyle\begin{split}\|\sup_{0\leq t\leq T}\lvert
Z(t)\rvert\|_{2}^{2}&=E\Bigg{[}\lvert\sup_{0\leq t\leq
T}\lvert\int_{0}^{t}\tilde{a}ds+\int_{0}^{t}bdW(s)+\int_{0}^{t}\int_{\mathbb{R}^{m}\backslash\\{0\\}}g(z)\eta\tilde{N}(ds,dz)\rvert\rvert^{2}\Bigg{]}\\\
&\leq E\Bigg{[}\lvert\sup_{0\leq t\leq
T}\left(\lvert\int_{0}^{t}\tilde{a}ds\rvert+\lvert\int_{0}^{t}bdW(s)\rvert+\lvert\int_{0}^{t}\int_{\mathbb{R}^{m}\backslash\\{0\\}}g(z)\eta\tilde{N}(ds,dz)\rvert\right)\rvert^{2}\Bigg{]}\\\
&\leq E\Bigg{[}\sup_{0\leq t\leq
T}\left(4\lvert\int_{0}^{t}\tilde{a}ds\rvert^{2}+4\lvert\int_{0}^{t}bdW(s)\rvert^{2}+2\lvert\int_{0}^{t}\int_{\mathbb{R}^{m}\backslash\\{0\\}}g(z)\eta\tilde{N}(ds,dz)\rvert^{2}\right)\Bigg{]}\\\
&\leq
E\Bigg{[}4T\int_{0}^{T}\tilde{a}^{2}ds+4\int_{0}^{T}b^{2}ds+2G^{2}(\infty)\int_{0}^{T}\eta^{2}ds\Bigg{]}\\\
&\leq 2(2T+2+G^{2}(\infty))\int_{0}^{T}E\left[K(1+\lvert
Z(s)\rvert^{2})\right]ds\\\ &\leq
2K(2T+2+G^{2}(\infty))\left(T+\int_{0}^{T}E\left[\sup_{0\leq t\leq s}\lvert
Z(t)\rvert^{2}\right]ds\right)\\\
&=C_{1}+\frac{C_{1}}{T}\int_{0}^{T}\|\sup_{0\leq t\leq s}\lvert
Z(t)\rvert\|_{2}^{2}ds,\end{split}$ (19)

where $C_{1}=2K(2T+2+G^{2}(\infty))/T$. Define $F(u)\triangleq\|\sup_{0\leq
t\leq u}\lvert Z(t)\rvert\|_{2}^{2}$. Then, by Grönwall’s inequality, we find
that

$\displaystyle F(T)\leq
C_{1}+\frac{C_{1}}{T}\int_{0}^{T}F(s)ds\quad\leftrightarrow\quad F(T)\leq
C_{1}e^{C_{1}}.$

Now, continue to prove the stated convergence result. By Doob’s maximal
inequality, the independence property of processes, the expectation rule of
Itô integrals, Itô isometry, monotone convergence and linear growth, we have
that

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
S(t)-S_{\epsilon,Z(t)}(t)\rvert\|_{2}^{2}$ $\displaystyle\leq
4E\left[\left(\int_{0}^{T}\int_{\lvert
z\rvert<\epsilon}g(z)\eta\tilde{N}(ds,dz)-\int_{0}^{T}G(\epsilon)\eta
dB(s)\right)^{2}\right]$ $\displaystyle=4E\left[\left(\int_{0}^{T}\int_{\lvert
z\rvert<\epsilon}g(z)\eta\tilde{N}(ds,dz)\right)^{2}+\left(\int_{0}^{T}G(\epsilon)\eta
dB(s)\right)^{2}\right]$
$\displaystyle=4E\left[\int_{0}^{T}G^{2}(\epsilon)\eta^{2}ds+\int_{0}^{T}G^{2}(\epsilon)\eta^{2}ds\right]$
$\displaystyle=8G^{2}(\epsilon)\int_{0}^{T}E\left[\eta^{2}\right]ds\leq
8G^{2}(\epsilon)\int_{0}^{T}E\left[K(1+\lvert Z(s)\rvert^{2})\right]ds$
$\displaystyle\leq 8G^{2}(\epsilon)K\left(T+\int_{0}^{T}E\left[\sup_{0\leq
t\leq s}\lvert Z(t)\rvert^{2}\right]ds\right)$
$\displaystyle=8G^{2}(\epsilon)K\left(T+\int_{0}^{T}\|\sup_{0\leq t\leq
s}\lvert Z(t)\rvert\|_{2}^{2}ds\right).$

Since $\|\sup_{0\leq t\leq s}\lvert Z(t)\rvert\|_{2}^{2}\leq\|\sup_{0\leq
t\leq T}\lvert Z(t)\rvert\|_{2}^{2}\leq C_{1}e^{C_{1}}$, we have that

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
S(t)-S_{\epsilon,Z(t)}(t)\rvert\|_{2}^{2}$ $\displaystyle\leq
8G^{2}(\epsilon)K\left(T+\int_{0}^{T}C_{1}e^{C_{1}}ds\right),$

and the proof is complete. ∎

Note that $G(\epsilon)$ is finite by definition, and that the bound in Lemma 2
converges to zero as $\epsilon\to 0$.

The convergence rate of the approximated jump diffusion $Z_{\epsilon}(t)$ (Eq.
(17)) is derived in the following proposition.

######  Proposition 2.

Let $Z(t)$ and $Z_{\epsilon}(t)$ be as given in Eq. (14) and (17)
respectively. The convergence rate of $Z_{\epsilon}(t)$ is given by

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
Z(t)-Z_{\epsilon}(t)\rvert\|_{2}\leq C_{\epsilon}^{1/2}G(\epsilon),$

where $C_{\epsilon}=8Ce^{C_{2}T}$,
$C_{2}=4K(T+1+G^{2}(\infty)+G^{2}(\epsilon))$, and $C$ and $G(\epsilon)$ is
given in Lemma 2 and Eq. (15) respectively.

###### Proof.

Define $f\triangleq f(t,Z(t))$ and $f_{\epsilon}\triangleq
f(t,Z_{\epsilon}(t))$. As in the proof of Lemma 2 (Eq. (19)), use the triangle
inequality, the identity $\lvert f+g\rvert^{p}\leq 2^{p-1}(\lvert
f\rvert^{p}+\lvert g\rvert^{p})$, Cauchy-Schwarz inequality, Itô isometry and
monotone convergence to obtain

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
Z(t)-Z_{\epsilon}(t)\rvert\|_{2}^{2}$
$\displaystyle=E\Bigg{[}\lvert\sup_{0\leq t\leq
T}\lvert\int_{0}^{t}(\tilde{a}-\tilde{a}_{\epsilon})ds+\int_{0}^{t}(b-b_{\epsilon})dW(s)+\int_{0}^{t}\int_{\lvert
z\rvert\geq\epsilon}g(z)(\eta-\eta_{\epsilon})\tilde{N}(ds,dz)$
$\displaystyle\quad\quad\quad+\int_{0}^{t}G(\epsilon)(\eta-\eta_{\epsilon})dB(s)+\int_{0}^{t}\int_{\lvert
z\rvert<\epsilon}g(z)\eta\tilde{N}(ds,dz)-\int_{0}^{t}G(\epsilon)\eta
dB(s)\rvert\rvert^{2}\Bigg{]}$ $\displaystyle\leq
4T\int_{0}^{T}E\left[\lvert\tilde{a}-\tilde{a}_{\epsilon}\rvert^{2}\right]ds+4\int_{0}^{T}E\left[\lvert
b-b_{\epsilon}\rvert^{2}\right]ds+4(G^{2}(\infty)-G^{2}(\epsilon))\int_{0}^{T}E\left[\lvert\eta-\eta_{\epsilon}\rvert^{2}\right]ds$
$\displaystyle\quad+8G^{2}(\epsilon)\int_{0}^{T}E\left[\lvert\eta-\eta_{\epsilon}\rvert^{2}\right]ds+8E\left[\sup_{0\leq
t\leq T}\lvert S(t)-S_{\epsilon,Z(t)}(t)\rvert^{2}\right].$

By Lipschitz continuity and Lemma 2 we further find

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
Z(t)-Z_{\epsilon}(t)\rvert\|_{2}^{2}$ $\displaystyle\leq
8CG^{2}(\epsilon)+4K(T+1+G^{2}(\infty)+G^{2}(\epsilon))\int_{0}^{T}E\left[\lvert
Z(s)-Z_{\epsilon}(s)\rvert^{2}\right]ds$
$\displaystyle=8CG^{2}(\epsilon)+C_{2}\int_{0}^{T}\|\sup_{0\leq t\leq s}\lvert
Z(t)-Z_{\epsilon}(t)\rvert\|_{2}^{2}ds,$

where $C_{2}=4K(T+1+G^{2}(\infty)+G^{2}(\epsilon))$. Define
$F(u)\triangleq\|\sup_{0\leq t\leq u}\lvert Z(t)\rvert\|_{2}^{2}$. Then, by
Grönwall’s inequality, we have

$\displaystyle F(T)\leq
8CG^{2}(\epsilon)+C_{2}\int_{0}^{T}F(s)ds\quad\leftrightarrow\quad F(T)\leq
8CG^{2}(\epsilon)e^{C_{2}T}.$

This concludes the proof. ∎

As for the convergence rate in Lemma 2, the convergence rate of
$Z_{\epsilon}(t)$ converges to zero as $\epsilon\to 0$.

Recall that the aim of this section is to derive an Euler scheme convergence
rate for the jump diffusion $Z(t)$ in Eq. (14) with infinite variations and
$\nu(U)=\infty$. The convergence rate of the Euler discretization of the
approximated process $Z_{\epsilon}(t)$ is therefore derived next. As we will
see in Section IV, this convergence result provides a foundation to estimate
MCARMA models driven by NIG-Lévy processes.

######  Proposition 3.

Let $Z(t)$ and $Z_{\epsilon}(t)$ be as given in Eq. (14) and (17)
respectively, and let
$Z_{\epsilon}^{h}(t)\triangleq\\{Z_{\epsilon}^{h}(t)\\}_{t\in[0,T]}$ denote
the Euler scheme of $Z_{\epsilon}(t)$ with equidistant time discretization
step size $h\in(0,1]$. The convergence rate of $Z_{\epsilon}^{h}(t)$ is given
by

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
Z(t)-Z_{\epsilon}^{h}(t)\rvert\|_{2}\leq\sqrt{2}\left(C_{\epsilon}G^{2}(\epsilon)+K_{3}h\right)^{1/2},$

where $K_{3}$ is a finite positive constant, and $C_{\epsilon}$ and
$G(\epsilon)$ is given in Proposition 2 and Eq. (15) respectively.

###### Proof.

We assume that $\tilde{a}(t,x)$ (and therefore also
$\tilde{a}_{\epsilon}(t,x)$), see Eq. (14), satisfy the Lipschitz and linear
growth conditions in Eq. (18). Then, since all conditions in Assumption 1 are
satisfied, Corollary 6.4.3 in [36] (see also [22]) gives

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
Z_{\epsilon}(t)-Z_{\epsilon}^{h}(t)\rvert\|_{2}^{2}\leq K_{3}h,$ (20)

where $K_{3}$ is a finite positive constant independent of $h$. Note that
Corollary 6.4.3 in [36] holds when the jump part of $Z_{\epsilon}(t)$ is a
compound Poisson process, which is the case for $Z(t)$ and $Z_{\epsilon}(t)$
when the drift is $\tilde{a}(t,x)$ and $\tilde{a}_{\epsilon}(t,x)$
respectively. Using the identity $\lvert f+g\rvert^{p}\leq 2^{p-1}(\lvert
f\rvert^{p}+\lvert g\rvert^{p})$, as in the proof of Lemma 2, we find

$\displaystyle\|\sup_{0\leq t\leq T}\lvert
Z(t)-Z_{\epsilon}^{h}(t)\rvert\|_{2}^{2}$ $\displaystyle=\|\sup_{0\leq t\leq
T}\lvert
Z(t)-Z_{\epsilon}(t)+Z_{\epsilon}(t)-Z_{\epsilon}^{h}(t)\rvert\|_{2}^{2}$
$\displaystyle\leq 2\|\sup_{0\leq t\leq T}\lvert
Z(t)-Z_{\epsilon}(t)\rvert\|_{2}^{2}+2\|\sup_{0\leq t\leq T}\lvert
Z_{\epsilon}(t)-Z_{\epsilon}^{h}(t)\rvert\|_{2}^{2}.$

By Proposition 2 and Eq. (20) the proof is complete. ∎

Note that the Euler scheme $Z_{\epsilon}^{h}(t)$ converges to $Z(t)$ as
$\epsilon\to 0$ and $h\to 0$.

It is straight forward to show that the results in Lemma 2, Proposition 2 and
Proposition 3 hold for MCARMA processes. A special case is considered in
Section IV.1.

##  IV The NIG-Lévy-driven MCAR process: A case study

The remaining of this work focus on a case study for an NIG-Lévy-driven MCAR
process, as an example of how to apply the multivariate ARMA/CARMA
transformation relation in model estimation. This MCAR modelling framework is
fit to a two-dimensional system of stratospheric temperature and wind
variables. More specifically, the dynamical system of dependent variables that
is to be considered is one of stratospheric temperature and U wind. As opposed
to temperature, wind constitutes a direction in addition to its magnitude. The
wind direction is given by two wind components, that is U wind representing
the west to east flow, and V wind representing the south to north flow. As
explained in, e.g., [28] and [26], it is particularly interesting to study the
U wind component of the Northern hemisphere, as extreme events such as sudden
stratospheric warmings has the potential to influence surface weather to a
great extent.

The initial stratospheric temperature and U wind data are retrieved as
European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim
atmospheric reanalysis model products ([9] and [17]), see Appendix S1. The
data analysed in this work are reprocessed as daily circumpolar mean
stratospheric temperature and U wind at $60^{\circ}$N and $10$ hPa altitude,
from 1 January 1979 to 31 December 2018, see Table 2. See [18] for more
information about data preparation.

In Section IV.1, we define the model representation and show that the
convergence results in Section III.2 hold, meaning that the multivariate
ARMA/CARMA transformation relation holds for this special case. Further, in
Section IV.2, some statistical considerations of the NIG-Lévy process are
assessed to properly express the error (moving average) coefficients of the
NIG-Lévy-driven MCAR process. Finally, in Section IV.3, an extended MCAR
process able to describe additive seasonality and heteroscedasticity in
dynamical systems is defined. Then, based on the multivariate ARMA/CARMA
transformation relation, explicit formulas for estimated autoregressive and
error coefficients are derived for a given set of model parameters. The
results are used to fit a two-dimensional NIG-Lévy-driven MCAR model to
stratospheric temperature and U wind data.

###  IV.1 Some considerations on the NIG-Lévy-driven MCAR(MA) process

In this section we explicitly define the MCAR process as a state space
representation, and show that the convergence results in Section III.2 hold
for an NIG-Lévy-driven MCAR(MA) process. This means that the multivariate
ARMA/CARMA transformation relation holds for the NIG-Lévy-driven MCAR model.

By the state space representation of Lévy-driven MCARMA processes in Section
II.1, the corresponding MCAR process is given by

$\displaystyle\begin{cases}\bm{Y}(t)=\bm{X}_{1}(t)\\\
d\bm{X}_{l}(t)=\mathbb{1}_{d}\bm{X}_{l+1}(t)dt\\\
d\bm{X}_{p}(t)=[-A_{p}\cdots-A_{1}]\bm{X}(t)dt+B_{0}d\bm{L}(t),\end{cases}$
(21)

see Eq. (4), (5) and (6). The multivariate ARMA/CARMA transformation relation
in Theorem 1 is considerably simplified for MCAR processes, as the shifted
Lévy-terms disappear. For any practical application, we assume the simpler
MCAR($p$) process to be a useful approximation of the MCARMA($p,q$) process.
That is, as stated in [24], for $p$ large enough, the VAR($p$) model
approaches a VARMA($p,q$) model. A rigorous proof of this statement for
MCAR/MCARMA processes is a topic for further research. The argument for
assuming this to hold for MCAR/MCARMA processes as well, is that the direct
transformation relation between VARMA/MCARMA processes (see Theorem 1) is
derived by introducing no structural changes going from the discretized MCARMA
process to the VARMA process. In particular, we will see in Section IV.3 that
an MCAR model explains well the two dimensional dynamical system of
stratospheric temperature and U wind.

We will continue to specify the MCAR process in Eq. (21) further. That is, we
will consider the framework when driven by an NIG-Lévy process. As shown in
[3], the generating triplet of an NIG-Lévy process is $(\bm{\alpha},0,\nu)$.
As the property $\int_{\lvert\bm{z}\rvert\leq
1}\lvert\bm{z}\rvert\nu(dz)=\infty$ holds for this process the small jumps
dominate, and the results in Section III.2 apply. By properties of the Lévy
measure, $\nu((-\epsilon,\epsilon)^{c})<\infty$ for all $0<\epsilon\leq 1$,
the NIG-Lévy process might be represented as (see Eq. (1))

$\displaystyle\bm{L}(t)=\hat{\bm{\alpha}}t+\int_{\lvert\bm{z}\rvert<\epsilon}\bm{z}\tilde{N}(t,d\bm{z})+\int_{\lvert\bm{z}\rvert\geq\epsilon}\bm{z}N(t,d\bm{z}),$
(22)

where
$\hat{\bm{\alpha}}=\left(\bm{\alpha}-\int_{\epsilon\leq\lvert\bm{z}\rvert<1}\bm{z}\nu(d\bm{z})\right)$.
Further, the MCARMA state space model in Eq. (4) written terms of the NIG-Lévy
process is given by

$\displaystyle
d\bm{X}(t)=\left(A\bm{X}(t)+\beta\hat{\bm{\alpha}}\right)dt+\beta\int_{\lvert\bm{z}\rvert<\epsilon}\bm{z}\tilde{N}(dt,d\bm{z})+\beta\int_{\lvert\bm{z}\rvert\geq\epsilon}\bm{z}N(dt,d\bm{z}),$
(23)

which is referred to as the NIG-Lévy-driven MCARMA process. By Proposition 2
and Proposition 3, it is straight forward to see that $\bm{X}_{\epsilon}(t)$
(corresponding to $Z_{\epsilon}(t)$ in Eq. (17)) converges to $\bm{X}(t)$ when
the small jumps approach zero, and that the corresponding Euler scheme
$\bm{X}_{\epsilon}^{h}(t)$ converges to $\bm{X}(t)$ as $h\to 0$, when all
relevant coefficients satisfy Assumption 1. That is, the transformation
relation in Theorem 1 is valid for the NIG-Lévy-driven MCARMA process, as
consequently for the NIG-Lévy-driven MCAR process.

###  IV.2 Distributional properties of the NIG-Lévy process

The distributional properties of linear combinations of NIG-Lévy processes are
considered in this section. The results will be used in Section IV.3 to
estimate an NIG-Lévy-driven MCAR model describing stratospheric temperature
and U wind dynamics. Note that an extended analysis including lags of the
driving NIG-Lévy process would have to be performed in order to consider
estimation of a corresponding MCARMA model.

By definition of the MCARMA process in Eq. (4), the moving average matrices
$\beta_{p-1},\ldots,\beta_{1}$ of the MCAR process are given by
$\mathbb{0}_{m}$. That is, the transformation relation in Theorem 1 gives a
component-wise stochastic part

$\displaystyle\Delta\mathcal{E}_{k}(t)\triangleq
stoch(x_{k}(t+p))=\sum_{r=1}^{m}\beta_{r,k}\Delta L_{r}(t),$ (24)

where $\beta_{i_{r},k}\triangleq\beta_{(k-1)m+i_{r}}^{(0)}$, and the
discretization step size is assumed to be $h=1$ (day). Notice that this
indicates that the driving process of the multivariate dynamical system (the
MCAR process) is given by $d\bm{\mathcal{E}}(t)=\beta d\bm{L}(t)$
($\beta\triangleq\beta_{p}$), with $\bm{\mathcal{E}}(t)$ being a
$d$-dimensional Lévy process.

Let the $1$-day increment of each component of the multidimensional Lévy
process, $\bm{L}(t)$, be distributed as univariate NIG random variables,
meaning

$\displaystyle\Delta
L_{i_{r}}(t)\overset{d}{\simeq}L_{i_{r}}(1)\overset{d}{\simeq}NIG(a_{i_{r}},b_{i_{r}},\delta_{i_{r}},\mu_{i_{r}}),$
(25)

where $a_{i_{r}}$, $b_{i_{r}}$, $\delta_{i_{r}}$ and $\mu_{i_{r}}$ is the tail
heaviness, asymmetry parameter, scale parameter and distributional location
respectively. The NIG distribution is closed under convolution and affine
transformations in the following sense (see, e.g., [4]): if $\Delta
L_{i_{r}}$, for $1\leq{i_{r}}\leq m$, are independent random variables that
are NIG distributed such that $a_{i_{r}}=a$ and $b_{i_{r}}=b$ for all
${i_{r}}$, then

$\displaystyle\sum_{r=1}^{m}\Delta
L_{r}(t)\overset{d}{\simeq}NIG\left(a,b,\sum_{r=1}^{m}\delta_{r},\sum_{r=1}^{m}\mu_{r}\right)$

and

$\displaystyle\beta_{i_{r},k}\Delta
L_{i_{r}}(t)\overset{d}{\simeq}NIG\left(\frac{a_{i_{r}}}{\lvert\beta_{{i_{r}},k}\rvert},\frac{b_{i_{r}}}{\beta_{{i_{r}},k}},\lvert\beta_{{i_{r}},k}\rvert\delta_{{i_{r}}},\beta_{{i_{r}},k}\mu_{{i_{r}}}\right).$

By these properties, the component-wise distribution generated by the Lévy
process, $\bm{\mathcal{E}}(t)$, is given by

$\displaystyle\mathcal{E}_{k}(t)\overset{d}{\simeq}NIG\left(a,b,\sum_{r=1}^{m}\lvert\beta_{r,k}\rvert\delta_{r},\sum_{r=1}^{m}\beta_{r,k}\mu_{r}\right),$
(26)

with restrictions
$a=a_{1}/\lvert\beta_{1,k}\rvert=\ldots=a_{m}/\lvert\beta_{m,k}\rvert$ and
$b=b_{1}/\beta_{1,k}=\ldots=b_{m}/\beta_{m,k}$, for all $1\leq k\leq d$. Note
that each independent NIG-Lévy process component, $dL_{i_{r}}(t)$, might
generate random variables from NIG distributions with distinct parameters
$a_{i}$ and $b_{i}$, where the components of $\beta$ have to be restricted for
Eq. (26) to hold.

###  IV.3 A two-dimensional MCAR($4$) case study

Estimation of a two-dimensional MCAR($4$) model driven by a two-dimensional
NIG-Lévy process is performed in this section. First, we will modify the MCAR
process, defined in Eq. (21), such that data with inherent seasonal behaviour
can be represented by the model. Then the transformation relation in Theorem 1
is used to explicitly state the autoregressive MCAR model coefficients as a
function of autoregressive VAR model coefficients, and restrictions are set on
the NIG distributions representing model residuals, such that formulas for the
MCAR model coefficients of the stochastic part can be derived. Finally,
empirical estimation of the modified MCAR model is exemplified using daily
circumpolar mean stratospheric temperature and U wind data, see Table 2 for
data specifications.

In [18], daily circumpolar mean stratospheric temperature at $60^{\circ}$ N
and $10$ hPa altitude were shown to follow a CAR($4$) process on the form

$\displaystyle\begin{cases}S(t)=\Lambda(t)+X_{1}(t)\\\
d\bm{X}(t)=A_{p}\bm{X}(t)dt+\bm{e}_{p}\sigma(t-)dL(t),\end{cases}$ (27)

with statistical significance. Here, $\Lambda(t)\in\mathbb{R}$ is a
deterministic, bounded and continuously differentiable seasonality function,
$\sigma(t)$ is a deterministic and càdlàg volatility function, and the
stochastic process $\bm{X}(t)$ in $\mathbb{R}^{p}$ represents deseasonalized
stratospheric temperature and its autoregressive nature. As indicated,
$\bm{X}(t)$ is given by a multidimensional non-Gaussian OU process with time
dependent volatility, where $A_{p}\in\mathbb{R}^{p}$, $\bm{e}_{p}$ is the unit
vector in $\mathbb{R}^{p}$, and the driving Lévy process takes values in
$\mathbb{R}$. This model is also used for modelling of temperature and wind in
the troposphere. For more details about this model, see, e.g., [18], [6], [8].

To represent weather variables independently with an one-dimensional model is
a significant simplification, as weather is a complex and non-linear system of
dependent variables. The MCAR process allows the improvement of representing
several weather variables as a linear cross-correlated system in an
autoregressive manner. Using a similar methodology as in [18], we see that
daily zonal mean stratospheric U wind, at $10$ hPa altitude, follows a
CAR($4$) process with statistical significance, just as stratospheric
temperature. This motivates the case study of fitting a two-dimensional
MCAR($4$) process to stratospheric temperature and U wind. Note that the
exercise of concluding an optimal value of the lag parameter, $p$, for the
MCAR process is not part of this work. The choice of using $p=4$ is made
because the marginals are well modelled in that case.

Inspired by the CAR process in Eq. (27), we redefine the MCAR model in Eq.
(21) as

$\displaystyle\begin{cases}\bm{Y}(t)=\bm{\Lambda}(t)+\bm{X}_{1}(t)\\\
d\bm{X}_{l}(t)=\mathbb{1}_{d}\bm{X}_{l+1}(t)dt\\\
d\bm{X}_{p}(t)=[-A_{p}\cdots-A_{1}]\bm{X}(t)dt+\bm{\sigma}(t-)\circ
B_{0}d\bm{L}(t),\end{cases}$ (28)

where $\bm{\Lambda}(t)\in{\mathbb{R}^{d}}$ is a deterministic and continuous
seasonality function, $\bm{\sigma}(t)\in{\mathbb{R}^{d}}$ is a deterministic
and càdlàg volatility function, and $\circ$ represents the Hadamard product.
This allows modelling data with seasonal and heteroscedastic behaviour. Note
that this extension of the MCAR model framework still might be represented as
in Eq. (23), meaning that the transformation relation in Theorem 1 is valid as
long as the relevant model coefficients satisfy Assumption 1. Considerations
of possible theoretical challenges regarding the spectral method of deriving
the MCARMA model, see [33], is beyond the scope of this paper. In what
follows, the two-dimensional NIG-Lévy-driven MCAR($4$) process in Eq. (28) is
fitted to stratospheric temperature and U wind data. Specifications of the
data are given in Table 2.

####  IV.3.1 The model coefficients

In this section, we show how autoregressive MCAR model coefficients are found
using the multivariate ARMA/CARMA transformation relation in Theorem 1.
Further, the error coefficients are computed through statistical properties of
the NIG distribution. This example is easily extended to other dimension and
lag parameters. The method of deriving error coefficients only holds for the
case of two driving NIG-Lévy processes.

The MCAR model parameters have to be given initially, and we choose parameters
as follows. We are considering a dynamical system of two variables, and so the
number of model dimensions is $d=2$. Further, it is reasonable to assume that
each dimension in the MCAR model is source of a random process, giving rise to
the idiosyncratic model error of each dimension separately. Therefore, we set
$m=d=2$. As seen in [24], it is more intricate to find an optimal lag
parameter for VARMA models than for the one-dimensional analogue, especially
when the driving process is not normal. As this case study is meant as a
demonstration of how to apply the transformations in Theorem 1 rather than a
detailed assessment of geophysical parameters, we are not concluding an
optimal value for $p$ in this work. As argued in the introduction of Section
IV.3, we set $p=4$.

Insert the parameters $p=4$, $d=2$ and $m=d$ into the multivariate ARMA/CARMA
transformation relation formula in Theorem 1 (remember that we assume $h=1$
and $q=0$), and rearrange to find the compressed expression

$\displaystyle\bm{x}(t+4)=\Phi_{1}\bm{x}(t+3)+\Phi_{2}\bm{x}(t+2)+\Phi_{3}\bm{x}(t+1)+\Phi_{4}\bm{x}(t)+\bm{\Phi}(\bm{x})+\beta\Delta\bm{L}(t),$
(29)

where $\bm{x}(\cdot)=[x_{1}(\cdot),x_{2}(\cdot)]$, $\Delta\bm{L}(t)=[\Delta
L_{1}(t),\Delta L_{2}(t)]$, and $\Phi_{i},\beta\in\mathbb{R}^{2\times 2}$ and
$\bm{\Phi}(\bm{x})\in\mathbb{R}^{2}$ are given by

$\displaystyle\Phi_{1}=\begin{bmatrix}-\alpha_{1,1}^{(1)}\hskip
5.69054pt&\hskip 5.69054pt-\alpha_{1,2}^{(1)}\\\ -\alpha_{2,1}^{(1)}\hskip
5.69054pt&\hskip
5.69054pt-\alpha_{2,2}^{(1)}\end{bmatrix},\Phi_{2}=\begin{bmatrix}-\alpha^{(2)}_{1,1}+3\alpha_{1,1}^{(1)}\hskip
5.69054pt&\hskip 5.69054pt-\alpha^{(2)}_{1,2}+3\alpha_{1,2}^{(1)}\\\
-\alpha^{(2)}_{2,1}+3\alpha_{2,1}^{(1)}\hskip 5.69054pt&\hskip
5.69054pt-\alpha^{(2)}_{2,2}+3\alpha_{2,2}^{(1)}\end{bmatrix},$
$\displaystyle\Phi_{3}=\begin{bmatrix}-\alpha^{(3)}_{1,1}+2\alpha_{1,1}^{(2)}-3\alpha_{1,1}^{(1)}\hskip
5.69054pt&\hskip
5.69054pt-\alpha^{(3)}_{1,2}+2\alpha_{1,2}^{(2)}-3\alpha_{1,2}^{(1)}\\\
-\alpha^{(3)}_{2,1}+2\alpha_{2,1}^{(2)}-3\alpha_{2,1}^{(1)}\hskip
5.69054pt&\hskip
5.69054pt-\alpha^{(3)}_{2,2}+2\alpha_{2,2}^{(2)}-3\alpha_{2,2}^{(1)}\end{bmatrix},$
$\displaystyle\Phi_{4}=\begin{bmatrix}-\alpha^{(4)}_{1,1}+\alpha_{1,1}^{(3)}-\alpha_{1,1}^{(2)}+\alpha^{(1)}_{1,1}\hskip
5.69054pt&\hskip
5.69054pt-\alpha^{(4)}_{1,2}+\alpha_{1,2}^{(3)}-\alpha_{1,2}^{(2)}+\alpha^{(1)}_{1,2}\\\
-\alpha^{(4)}_{2,1}+\alpha_{2,1}^{(3)}-\alpha_{2,1}^{(2)}+\alpha^{(1)}_{2,1}\hskip
5.69054pt&\hskip
5.69054pt-\alpha^{(4)}_{2,2}+\alpha_{2,2}^{(3)}-\alpha_{2,2}^{(2)}+\alpha^{(1)}_{2,2}\end{bmatrix},$
$\displaystyle\bm{\Phi}(\bm{x})=\begin{bmatrix}\sum_{i=0}^{3}k(i)x_{1}(t+i)\\\
\sum_{i=0}^{3}k(i)x_{2}(t+i)\end{bmatrix}\quad\text{and}\quad\beta=\begin{bmatrix}\beta_{1,1}&\beta_{1,2}\\\
\beta_{2,1}&\beta_{2,2}\end{bmatrix},$

where we defined
$\alpha_{k,i_{s}}^{(j)}\triangleq\alpha_{(k-1)d+i_{s}}^{(p-l+1)}$ and
$\beta_{k,i_{r}}\triangleq\beta_{(k-1)m+i_{r}}^{(0)}$, and
$[k(0),k(1),k(2),k(3)]=[4,-6,4,-1]$.

Compare the formula in Eq. (29) with a two-dimensional VAR($4$) model,
generally defined in Eq. (2). That is,

$\displaystyle\bm{x}(t+4)=\phi_{1}\bm{x}(t+3)+\phi_{2}\bm{x}(t+2)+\phi_{3}\bm{x}(t+1)+\phi_{4}\bm{x}(t)+\mathbb{1}_{2}\bm{\mathcal{E}}(t),$
(30)

where $\bm{\mathcal{E}}(t)\in\mathbb{R}^{d}$ represents a sequence of serially
uncorrelated i.i.d. random vectors and

$\displaystyle\phi_{j}=\begin{bmatrix}\phi^{(j)}_{11}&\phi^{(j)}_{12}\\\
\phi^{(j)}_{21}&\phi^{(j)}_{22}\end{bmatrix}.$

With the discretized MCAR process and the VAR process expressed as in Eq. (29)
and (30) respectively, it is straight forward to compute the MCAR model
coefficients, $\alpha_{k,i_{s}}^{(j)}$, as a function of the VAR model
coefficients, $\phi^{j}_{\cdot\cdot}$. That is, still with attention on the
autoregressive model coefficients, we have for example that

$\displaystyle\begin{bmatrix}\alpha_{1,1}^{(1)}&\alpha_{1,2}^{(1)}\\\
\alpha_{2,1}^{(1)}&\alpha_{2,2}^{(1)}\end{bmatrix}=\begin{bmatrix}-\phi_{11}^{(1)}-1&-\phi_{12}^{(1)}\\\
-\phi_{21}^{(1)}&-\phi_{22}^{(1)}-1\end{bmatrix}.$ (31)

Continuing with similar approach for $\Phi_{2}$, $\Phi_{3}$, and then
$\Phi_{4}$, with some substitutions, we obtain all $16$ autoregressive MCAR
model coefficients as a function of the fitted autoregressive VAR model
coefficients. The results are given in Appendix D.

When fitting a VAR model to data, the distributions of components of the
random vector $\bm{\mathcal{E}}(t)$ are easily obtained by fitting an NIG
distribution to the model residuals. This is done component-wise by the same
methodology as in [18]. Compare the VAR and MCAR representations in Eq. (29)
and (30), respectively, to see that the error coefficients of the MCAR model
are given by

$\displaystyle\mathcal{E}_{k}(t)=\beta_{k,1}\Delta L_{1}(t)+\beta_{k,2}\Delta
L_{2}(t),$ (32)

where we use the notation
$\bm{\mathcal{E}}(t)=[\mathcal{E}_{1}(t),\mathcal{E}_{2}(t)]$, and where
$\Delta L_{1}(t)$ and $\Delta L_{2}(t)$ represents the idiosyncratic error of
stratospheric temperature and U wind respectively. Remember that an MCAR model
with $m=2$ is considered, and that $\mathcal{E}_{k}(t)$ and $\Delta
L_{i_{r}}(t)$ ($k,i_{r}\in\\{1,2\\}$) are assumed to be NIG distributed random
variables. Note that $\beta$ is simply the identity matrix $\mathbb{1}_{2}$ if
the residual datasets represented by $\mathcal{E}(t)$ are independent, as the
fitted distributions are given as $\hat{\mathcal{E}}_{1}(t)\simeq\Delta
L_{1}(t)$ and $\hat{\mathcal{E}}_{2}(t)\simeq\Delta L_{2}(t)$ in that case. If
the residual datasets are dependent, the question is more intricate.

As we will see, to obtain a solution of the linear system in Eq. (32) (if it
exists) for $k\in\\{1,2\\}$, the idiosyncratic error distributions generated
by $\Delta L_{1}(t)$ and $\Delta L_{2}(t)$ have to be restricted. We assume
that the idiosyncratic error distributions are given as in Eq. (25), and that
the residual datasets for stratospheric temperature and U wind are NIG-
distributed as
$\mathcal{E}_{k}(t)\overset{d}{\simeq}NIG(a_{k}^{\mathcal{E}},b^{\mathcal{E}}_{k},\delta^{\mathcal{E}}_{k},\mu^{\mathcal{E}}_{k})$.
The distribution parameters of $\mathcal{E}_{k}(t)$ are naturally restricted
by the explicitly given parameters in Eq. (26).

Denote the two-dimensional NIG-Lévy process driving the two-dimensional
MCAR($4$) model as $\bm{L}(t)=[L_{1}(t),L_{2}(t)]$. Further, assume that the
elements of $\bm{L}(t)$ are independent random variables with
$E[\bm{L}(t)]=[0,0]$ and $\text{Var}(\bm{L}(t))=[1,1]$. These assumptions are
reasonable because we are considering the MCAR model in Eq. (28), meaning that
the model is shifted and scaled by the seasonally varying functions
$\bm{\Lambda}(t)$ and $\bm{\sigma}(t)$, respectively. The goal from here is to
estimate the parameters in $\beta$ from NIG-distributions fitted to the
residual datasets, which is represented by $\bm{\mathcal{E}}(t)$. The above
assumptions give $\text{Cov}(\bm{L}(t))=\mathbb{1}_{2}$, meaning that the
covariance of $\bm{\mathcal{E}}(t)$ is given as

$\displaystyle\Sigma\triangleq\text{Cov}(\bm{\mathcal{E}}(t))=\text{Cov}(\beta\Delta\bm{L}(t))=\beta\text{Cov}(\Delta\bm{L}(t))\beta^{\prime}=\beta\beta^{\prime}=\begin{bmatrix}\beta_{1,1}^{2}+\beta_{1,2}^{2}&\beta_{1,1}\beta_{2,1}+\beta_{1,2}\beta_{2,2}\\\
&\beta_{2,1}^{2}+\beta_{2,2}^{2}\end{bmatrix}.$ (33)

This formula gives a direct relation between $\beta$ and the empirically
computed covariance matrix, $\hat{\Sigma}$, of model residuals. As the
covariance matrix is symmetric (see Eq. (33)), an additional restriction is
needed to derive explicit formulas for the components of $\beta$. The
additional restriction is set on the idiosyncratic error distributions
generated by $\bm{L}(t)$. The goal is to set a restriction that is simple, but
still leave the model as flexible as possible. Consider the following
assessment of possible restrictions, and how they affect the NIG distribution
parameters of $\bm{\mathcal{E}}(t)$:

  1. 1.

Assuming $a_{i_{r}}=C_{a}$, $i_{r}\in\\{1,2\\}$, for a constant $C_{a}\geq 0$
gives $\beta_{k,i_{r}}=C_{a}/\hat{a}_{k}^{\mathcal{E}}$ for all $i_{r}$, where
$\hat{a}_{k}^{\mathcal{E}}$ is the estimated tail heaviness for the
distribution of $\mathcal{E}_{k}(t)$. This is not a beneficial restriction as
it forces $\beta_{k,1}=\beta_{k,2}$ for each dimension $k$.

  2. 2.

A similar argument as in point 1. holds for the idiosyncratic error
distribution asymmetry parameters $b_{i_{r}}$.

  3. 3.

The assumption $E[\bm{L}(t)]=[0,0]$ gives
$\mu_{i_{r}}=-\delta_{i_{r}}b_{i_{r}}/\sqrt{a_{i_{r}}^{2}-b_{i_{r}}^{2}}$,
$i_{r}\in\\{1,2\\}$, by definition. Therefore, the restriction $\mu_{i_{r}}=0$
forces $b_{i_{r}}$ to be zero and vice versa, since $\delta_{i_{r}}>0$ by
definition. This leads to strict conditions on the idiosyncratic error
distributions.

The restrictions listed above increase the risk of non-existing statistical
significant NIG distributions for the elements of $\bm{\mathcal{E}}(t)$. Based
on this, there are only two reasonable (and simple enough) choices of
restrictions to the idiosyncratic error distributions. Either choose
$\mu_{i_{r}}=C_{\mu}$ for a constant $\lvert C_{\mu}\rvert>0$ and
$i_{r}\in\\{1,2\\}$, or $\delta_{i_{r}}=C_{\delta}$ for a constant
$C_{\delta}>0$ and $i_{r}\in\\{1,2\\}$. That is, either restrict the scale
parameters or the location parameters of the idiosyncratic error distributions
to be equal.

In this study, we continue with the restriction $\delta_{i_{r}}=C_{\delta}$
for a constant $C_{\delta}>0$ and $i_{r}\in\\{1,2\\}$. Thus, by Eq. (26) we
have that

$\displaystyle\begin{cases}\delta_{1}^{\mathcal{E}}=\lvert\beta_{1,1}\rvert
C_{\delta}+\lvert\beta_{1,2}\rvert C_{\delta}\\\
\delta_{2}^{\mathcal{E}}=\lvert\beta_{2,1}\rvert
C_{\delta}+\lvert\beta_{2,2}\rvert C_{\delta}\end{cases}.$ (34)

To find the two-dimensional MCAR($4$) model coefficients in $\beta$, the
system of equations in Eq. (33) and (34) has to be solved. Notice that the
system of equations are independent of the lag parameter, $p$.

Solving the system of equations in Eq. (33) and (34) gives

$\displaystyle\pm\sqrt{\Sigma_{11}-\beta_{1,2}^{2}}\beta_{2,1}\pm\sqrt{\Sigma_{22}-\beta_{2,1}^{2}}\beta_{1,2}=\Sigma_{12},$
(35)

where

$\displaystyle\beta_{1,2}=\frac{1}{2C_{\delta}}\left(\pm\delta_{1}^{\mathcal{E}}\pm\sqrt{2\Sigma_{11}C_{\delta}^{2}-(\delta_{1}^{\mathcal{E}})^{2}}\right)\quad\text{and}\quad\beta_{2,1}\in\frac{1}{2C_{\delta}}\left(\pm\delta_{2}^{\mathcal{E}}\pm\sqrt{2\Sigma_{22}C_{\delta}^{2}-(\delta_{2}^{\mathcal{E}})^{2}}\right),$

for all four possible combinations of signs in each case. For convenience of
the reader, the derivations are included in Appendix E.

To conclude, all components of the coefficient matrix $\beta$ are dependent on
the idiosyncratic error distribution scale parameter, $C_{\delta}$. Any real-
valued solution of components of $\beta$ are allowed as long as Eq. (35) is
satisfied, including the restrictions

$\displaystyle 1)\;C_{\delta}>0;\quad
2)\;C_{\delta}\geq\delta_{1}^{\mathcal{E}}/\sqrt{2\Sigma_{11}},C_{\delta}\geq\delta_{2}^{\mathcal{E}}/\sqrt{2\Sigma_{22}};\quad
3)\;\Sigma_{11}\geq\beta_{2,1}^{2},\Sigma_{22}\geq\beta_{1,2}^{2}.$ (36)

####  IV.3.2 Fit model to stratospheric temperature and U wind data

Empirical estimation of coefficients of a two-dimensional MCAR($4$) model
driven by a two-dimensional NIG-Lévy process is performed using results from
Section IV.3.1. The model is fit to daily zonal mean stratospheric temperature
and U wind data. See Table 2 for data specifications, and note that the
datasets are prepared as presented in [18].

Table 2: Specifications of stratospheric temperature and U wind datasets. The specifications for the two datasets are similar, except units. Initial datasets are retrieved as ECMWF ERA-Interim reanalysis model products.  Date | Grid | Pressure level | Time | Area | Unit
---|---|---|---|---|---
1 January 1979 to 31 December 2018 |  $0.5^{\circ}$ |  $10$ hPa |  00:00, 06:00, 12:00, 18:00 |  $60^{\circ}$N and $[-180^{\circ}\text{E},180^{\circ}\text{E})$ |  Temp.: Kelvin, U wind: m/s

As argued in the introduction of Section IV.3, it is reasonable to assume that
the dynamical system of stratospheric temperature and U wind follow a VAR/MCAR
model. In the following, we will work with the two-dimensional state space
model in Eq. (28). That is, the dynamical system is given by
$\bm{Y}(t)=[Y_{1}(t),Y_{2}(t)]$, where $Y_{1}(t)$ and $Y_{2}(t)$ represents
temperature and U wind respectively. As explained in Section IV.3.1, the
parameter values of the MCAR model are initially given as $d=2$, $m=d=2$ and
$p=4$ (remember that $q=0$ in a MCAR model). We use a similar methodology as
in [18] to fit a VAR model, equivalent to the modified MCAR model in Eq. (28),
to temperature and U wind datasets. That is, with all following steps
performed component-wise for $Y_{1}(t)$ and $Y_{2}(t)$, except when fitting
the VAR model, the methodology goes as

  1. 1.

Fit a continuous seasonality function
$\bm{\Lambda}(t)=[\Lambda_{1}(t),\Lambda_{2}(t)]$ to datasets representing
$Y_{1}(t)$ and $Y_{2}(t)$. Deseasonalized datasets, represented by $X_{1}(t)$
and $X_{1}(t)$, are obtained by subtracting the fitted functions,
$\Lambda_{1}(t)$ and $\Lambda_{2}(t)$, from the datasets.

  2. 2.

Fit a VAR($p$) model to the (two-dimensional) deseasonalized dataset using a
standard statistical programming package. Subtract the fitted VAR($p$) model
from the deseasonalized datasets (component-wise) to obtain datasets of model
residuals.

  3. 3.

Compute empirically the expected values of squared residuals each day over the
year (assumed to be $365$ days as each February 29 is removed for convenience)
to construct the bivariate volatility function,
$\bm{\sigma}(t)=[\sigma_{1}(t),\sigma_{2}(t)]$.

  4. 4.

Scale the residual datasets with the constructed volatility function
components to obtain datasets of $\bm{\sigma}(t)$-scaled model residuals, to
which NIG distributions are fitted. These distributions are said to be
generated by the random vector $\bm{\mathcal{E}}(t)$.

  5. 5.

The MCAR model coefficients are computed from formulas derived in Section
IV.3.1.

As in [18], the seasonality function components of $\bm{\Lambda}(t)$ are
assumed to be given by

$\displaystyle\Lambda_{k}(t)=c_{0}^{(k)}+c_{1}^{(k)}t+\sum_{j=1}^{10}\left(c_{2j}^{(k)}\cos(j\pi
t/365)+c_{2j+1}^{(k)}\sin(j\pi t/365)\right),$ (37)

where $k=1,2$ gives the dimension and
$c_{0}^{(k)},c_{1}^{(k)},c_{2}^{(k)},\ldots,c_{21}^{(k)}$ are constants. Note
that $c_{0}^{(k)}$ captures the initial mean value, $c_{1}^{(k)}$ captures the
slope of the long-term liner changes, and
$c_{0}^{(k)},c_{1}^{(k)},c_{2}^{(k)},\ldots,c_{21}^{(k)}$ describe the yearly
cycle as weights in the truncated Fourier series. The ten last years of the
stratospheric temperature and U wind datasets are shown with a fitted
seasonality function in Figure 1, and the corresponding parameter values of
Eq. (37) are listed in Table 3.

(a) Stratospheric temperature

(b) Stratospheric U wind

Figure 1: Daily circumpolar mean stratospheric temperature and U wind from 1
January 2009 to 31 December 2018 with fitted seasonality function
$\Lambda_{1}(t)$ and $\Lambda_{2}(t)$, respectively. Table 3: Coefficients of
seasonality functions $\Lambda_{1}(t)$ and $\Lambda_{2}(t)$ fitted to
stratospheric temperature and U wind.  Seasonality function parameters:
Temperature
---
$c_{0}$ | $c_{2}$ | $c_{4}$ | $c_{6}$ | $c_{8}$ | $c_{10}$ | $c_{12}$ | $c_{14}$ | $c_{16}$ | $c_{18}$ | $c_{20}$
$226.15$ | $-0.049$ | $-12.094$ | $0.23$ | $1.88$ | $0.33$ | $0.16$ | $0.13$ | $-0.094$ | $-0.15$ | $-0.014$
$c_{1}$ | $c_{3}$ | $c_{5}$ | $c_{7}$ | $c_{9}$ | $c_{11}$ | $c_{13}$ | $c_{15}$ | $c_{17}$ | $c_{19}$ | $c_{21}$
$-0.000072$ | $-0.11$ | $1.63$ | $-0.23$ | $2.81$ | $-0.040$ | $1.54$ | $0.14$ | $0.45$ | $0.049$ | $0.11$
Seasonality function parameters: U wind
$c_{0}$ | $c_{2}$ | $c_{4}$ | $c_{6}$ | $c_{8}$ | $c_{10}$ | $c_{12}$ | $c_{14}$ | $c_{16}$ | $c_{18}$ | $c_{20}$
$11.18$ | $0.19$ | $22.58$ | $-0.40$ | $0.88$ | $-0.29$ | $1.12$ | $0.59$ | $1.06$ | $0.66$ | $0.95$
$c_{1}$ | $c_{3}$ | $c_{5}$ | $c_{7}$ | $c_{9}$ | $c_{11}$ | $c_{13}$ | $c_{15}$ | $c_{17}$ | $c_{19}$ | $c_{21}$
$-0.00011$ | $0.17$ | $-4.29$ | $0.015$ | $-0.33$ | $-0.69$ | $0.18$ | $-0.57$ | $-0.37$ | $0.19$ | $0.31$

Studying crosscorrelations of deseasonalized stratospheric temperature and U
wind confirms dependence in daily lagged values between the two variates, see
Figure 2. This supports the assumption that an MCAR process might represent
the two-dimensional system well. A VAR model is fit to the deseasonalized
dataset, corresponding to $[Y_{1}(t)-\Lambda_{1}(t),Y_{2}(t)-\Lambda_{2}(t)]$
(see Eq. (28)), using the add-on SSMMATLAB in MATLAB, see [25] and [24]. The
fitted VAR coefficients are given in Table 4 with corresponding t-values. Most
of the significance levels are acceptable. Note that the constants, $\bm{c}$,
of the VAR model are taken as zeroes, as deseasonalized data are considered.

Figure 2: Autocorrelation and crosscorrelation of deseasonalized stratospheric
temperature and U wind. Table 4: Coefficients of VAR model fitted to the two-
dimensional dataset of deseasonalized stratospheric temperature and U wind,
with corresponding t-values.  VAR model coefficients and corresponding
t-values
---
$\bm{c}$ | $\hat{\phi}_{1}$ | $\hat{\phi}_{2}$ | $\hat{\phi}_{3}$ | $\hat{\phi}_{4}$
$\begin{bmatrix}0.0003\\\ -0.0014\end{bmatrix}$ | $\begin{bmatrix}1.53&0.0008\\\ -0.10&1.73\end{bmatrix}$ | $\begin{bmatrix}-0.73&0.024\\\ -0.097&-1.03\end{bmatrix}$ | $\begin{bmatrix}0.27&0.0029\\\ 0.14&0.32\end{bmatrix}$ | $\begin{bmatrix}-0.10&0.023\\\ 0.0026&-0.041\end{bmatrix}$
$t_{\bm{c}}$ | $t_{\phi_{1}}$ | $t_{\phi_{2}}$ | $t_{\phi_{3}}$ | $t_{\phi_{4}}$
$\begin{bmatrix}0.058\\\ -0.12\end{bmatrix}$ | $\begin{bmatrix}186.43&0.25\\\ -4.94&208.67\end{bmatrix}$ | $\begin{bmatrix}-48.83&-3.78\\\ -2.56&-63.25\end{bmatrix}$ | $\begin{bmatrix}17.88&0.45\\\ 3.79&19.37\end{bmatrix}$ | $\begin{bmatrix}-12.18&6.91\\\ 0.13&-4.92\end{bmatrix}$

The yearly varying multivariate volatility function $\bm{\sigma}(t)$ is
estimated from the model residuals. As in [18], each volatility function,
$\sigma_{k}(t)$ for $k=1,2$, is constructed from three truncated Fourier
series on the form

$\displaystyle
w_{f_{i},k}^{(n_{i})}(t)=d_{0,k}^{(i)}+\sum_{j=1}^{n_{i}}\left(d_{2j-1,k}^{(i)}\cos(f_{i}j\pi
t/365)+d_{2j,k}^{(i)}\sin(f_{i}j\pi t/365)\right),$

where $k=1,2$ gives the dimension, $i\in\\{1,2,3\\}$ represents the series
number, $f_{i}$ adjusts the series frequency, $n_{i}$ gives the number of
terms in each series, and $d_{0,k}^{(i)},d_{2j-1,k}^{(i)},d_{2j,k}^{(i)}$ are
constants. Each of the three fitted series, for each $k$, represents seasons
winter/spring, summer, and autumn/winter, defined by the intervals [1
January,30 April], [1 May,31 October] and [1 November,31 December],
respectively. See [18] for an explanation of why. The final volatility
function $\sigma_{k}(t)$ is constructed by connecting each of the three fitted
functions $w_{f_{i},k}^{(n_{i})}(t)$ using sigmoid functions. That is, let the
sigmoid function, $\omega(x)$, and the connective function, $\xi(x)$, be given
by

$\displaystyle\omega(x)=\frac{1}{1+\exp\left(-(\frac{x-a}{b})\right)}\quad\text{and}\quad\xi(x)=\big{(}1-\omega(x)\big{)}g_{1}(x)+\omega(x)g_{2}(x),$

where $a$ and $b$ are shift and scaling constants respectively, and $g_{1}(x)$
and $g_{2}(x)$ are two functions that are to be connected. The numerical
results are as follows. For stratospheric temperature, the functions
$w_{0.44,1}^{(2)}(t)$ and $w_{0.30,2}^{(2)}(t)$ are connected with
$w_{2.0,1}^{(2)}(t)$ and $w_{0.50,2}^{(3)}(t)$ respectively, using parameters
$a=120$ and $b=2$. Further, for stratospheric U wind, the functions
$w_{2.0,1}^{(2)}(t)$ and $w_{0.50,2}^{(3)}(t)$ are connected with
$w_{0.44,1}^{(2)}(t)$ and $w_{0.05,2}^{(4)}(t)$ respectively, using parameters
$a=303$ and $b=5$. These connections make up two smooth functions
$\sigma_{1}(t)$ and $\sigma_{2}(t)$, that are illustrated in Figure 3 as the
variance function $V_{k}(t)=\sigma^{2}_{k}(t)$, together with empirically
estimated daily variance values over the year. The corresponding coefficients
of the truncated Fourier series are given in Table 5 and Table 6 for
stratospheric temperature and U wind respectively.

(a) Stratospheric temperature

(b) Stratospheric U wind

Figure 3: Estimated expected squared residuals ($\simeq$ variance) each day of
the year and a fitted variance function. Table 5: Coefficients of the
volatility function $\sigma_{1}(t)$ fitted to VAR model residuals of
stratospheric temperature.  Winter/spring: $w_{0.44,1}^{(2)}$
---
$d_{0,1}^{(1)}$ | $d_{1,1}^{(1)}$ | $d_{2,1}^{(1)}$ | $d_{3,1}^{(1)}$ | $d_{4,1}^{(1)}$
$-595.152$ | $749.200$ | $289.857$ | $-153.179$ | $-140.785$
Summer: $w_{2.0,1}^{(2)}$
$d_{0,1}^{(2)}$ | $d_{1,1}^{(2)}$ | $d_{2,1}^{(2)}$ | $d_{3,1}^{(2)}$ | $d_{4,1}^{(2)}$
$0.089$ | $0.103$ | $0.022$ | $0.033$ | $0.014$
Autumn/winter: $w_{0.44,1}^{(2)}$
$d_{0,1}^{(3)}$ | $d_{1,1}^{(3)}$ | $d_{2,1}^{(3)}$ | $d_{3,1}^{(3)}$ | $d_{4,1}^{(3)}$
$13.000$ | $-228.138$ | $75.488$ | $80.783$ | $87.432$
Table 6: Coefficients of the volatility function $\sigma_{2}(t)$ fitted to VAR
model residuals of stratospheric U wind.  Winter/spring: $w_{0.30,2}^{(2)}$
---
$d_{0,2}^{(1)}$ | $d_{1,2}^{(1)}$ | $d_{2,2}^{(1)}$ | $d_{3,2}^{(1)}$ | $d_{4,2}^{(1)}$ | $-$ | $-$ | $-$ | $-$
$-32.767$ | $-333.062$ | $1960.769$ | $370.204$ | $-945.071$ | $-$ | $-$ | $-$ | $-$
Summer: $w_{0.50,2}^{(3)}$
$d_{0,2}^{(2)}$ | $d_{1,2}^{(2)}$ | $d_{2,2}^{(2)}$ | $d_{3,2}^{(2)}$ | $d_{4,2}^{(2)}$ | $d_{5,2}^{(2)}$ | $d_{6,2}^{(2)}$ | $-$ | $-$
$22.726$ | $-127.933$ | $89.798$ | $126.536$ | $2.515$ | $-24.633$ | $-22.122$ | $-$ | $-$
Autumn/winter: $w_{0.05,2}^{(4)}$
$d_{0,2}^{(3)}$ | $d_{1,2}^{(3)}$ | $d_{2,2}^{(3)}$ | $d_{3,2}^{(3)}$ | $d_{4,2}^{(3)}$ | $d_{5,2}^{(3)}$ | $d_{6,2}^{(3)}$ | $d_{7,2}^{(3)}$ | $d_{8,2}^{(3)}$
$113.995$ | $93.302$ | $20.305$ | $33.465$ | $22.639$ | $-58.418$ | $-7.422$ | $-170.865$ | $-83.501$

The $\bm{\sigma}(t)$-scaled VAR model residual datasets are shown to be NIG-
distributed with statistical significance. The NIG parameters for
stratospheric temperature and U wind are listed in Table 7 with corresponding
KS test statistics and p-values. The autocorrelations and crosscorrelations of
$\bm{\sigma}(t)$-scaled VAR model residuals are shown in Figure 4. It is
apparent that the two-dimensional VAR($4$) model explains most of the time
lagged dependence. As discussed in [18], the remaining memory effects indicate
that a stochastic volatility model would provide a more accurate model.

Table 7: Fitted NIG parameters for $\bm{\sigma}(t)$-scaled VAR model residuals
and corresponding KS test results.  Stratospheric temperature
---
$a_{1}^{\mathcal{E}}$ | $b_{1}^{\mathcal{E}}$ | $\delta_{1}^{\mathcal{E}}$ | $\mu_{1}^{\mathcal{E}}$ | statistic | p-value
$2.93$ | $0.398$ | $1.70$ | $-0.234$ | $0.007$ | $0.47$
Stratospheric U wind
$a_{2}^{\mathcal{E}}$ | $b_{2}^{\mathcal{E}}$ | $\delta_{2}^{\mathcal{E}}$ | $\mu_{2}^{\mathcal{E}}$ | statistic | p-value
$3.13$ | $-0.0781$ | $1.77$ | $0.0471$ | $0.005$ | $0.90$
Figure 4: Autocorrelation and crosscorrelation of squared
$\bm{\sigma}(t)$-scaled model residuals of stratospheric temperature and U
wind.

The autoregressive MCAR model coefficients are found from the VARMA/MCARMA
transformation relation in Theorem 1, that is explicitly given for our model
framework in Section IV.3.1. The empirical results are shown in Table 8, where
the autoregressive coefficients are displayed as the matrix components of the
MCAR model block matrix $A$, see Eq. (4) and (7).

Table 8: Estimated autoregressive coefficients of the MCAR model fit to
stratospheric temperature and U wind.  Autoregressive MCAR model coefficients
---
$\hat{A}_{4}$ | $\hat{A}_{3}$ | $\hat{A}_{2}$ | $\hat{A}_{1}$
$\begin{bmatrix}0.030&-0.0019\\\ 0.054&0.029\end{bmatrix}$ | $\begin{bmatrix}-4.41&0.043\\\ 0.36&-4.43\end{bmatrix}$ | $\begin{bmatrix}-2.87&0.022\\\ 0.41&-3.15\end{bmatrix}$ | $\begin{bmatrix}-2.53&-0.00083\\\ 0.10&-2.72\end{bmatrix}$

As mentioned in Section II.1, the MCAR process has a unique causal stationary
solution when the eigenvalues of the block matrix $A$ have negative real
parts. As reported in Table 9, we see that the eigenvalues
$\lambda^{(A)}_{3,4}$ fail having this property. Following [23], we know that
the VAR process has a unique causal stationary solution if and only if the
modulus of each root of the polynomial $\det(\phi(z))$ is greater than one.
Note that
$\phi(z)=\left(I-\phi_{1}z-\phi_{2}z^{2}-\phi_{3}z^{3}-\phi_{4}z^{4}\right)$
in this case, see Eq. (2). Further, [23] states that when the VAR process is
written in a so-called Akaike’s state space form, one can show that

$\displaystyle\det(\phi(\lambda))=\det(\mathbb{1}_{d}-F\lambda),\quad\text{where}\quad
F=\begin{bmatrix}\mathbb{0}_{d}&\mathbb{1}_{d}&\mathbb{0}_{d}&\mathbb{0}_{d}\\\
\mathbb{0}_{d}&\mathbb{0}_{d}&\mathbb{1}_{d}&\mathbb{0}_{d}\\\
\mathbb{0}_{d}&\mathbb{0}_{d}&\mathbb{0}_{d}&\mathbb{1}_{d}\\\
\phi_{4}&\phi_{3}&\phi_{2}&\phi_{1}\\\ \end{bmatrix}.$

This means that the above condition on existence of a proper VAR process
solution is equivalent to the modulus of each eigenvalue of $F$ being less
than one. The eigenvalues of $F$ are listed in Table 9 with their respective
modulus. All eight eigenvalues have modulus less than one, and thus the fitted
VAR process has a unique causal stationary solution.

Table 9: Eigenvalues of the MCAR state space representation matrix, $A$, and
of the VAR Akaike’s state space form matrix, $F$.  Eigenvalues of $A$
---
$\lambda_{1}^{(A)}$ | $\lambda_{2}^{(A)}$ | $\lambda_{3,4}^{(A)}$ | $\lambda_{5,6}^{(A)}$ | $\lambda_{7,8}^{(A)}$
$-2.21$ |  $-2.15$ |  $0.0067\pm 0.0021i$ |  $-0.20\pm 1.42i$ |  $-0.25\pm 1.40i$
Eigenvalues of $F$
$\lambda_{1}^{(F)}$ | $\lambda_{2}^{(F)}$ | $\lambda_{3,4}^{(F)}$ | $\lambda_{5,6}^{(F)}$ | $\lambda_{7,8}^{(F)}$
$0.18$ | $0.68$ | $-0.024\pm 0.41i$ | $0.28\pm 0.38i$ | $0.94\pm 0.035i$
Modulus of eigenvalues of $F$
$\lvert\lambda_{1}^{(F)}\rvert$ | $\lvert\lambda_{2}^{(F)}\rvert$ | $\lvert\lambda_{3,4}^{(F)}\rvert$ | $\lvert\lambda_{5,6}^{(F)}\rvert$ | $\lvert\lambda_{7,8}^{(F)}\rvert$
$0.18$ | $0.68$ | $0.41$ | $0.47$ | $0.94$

The above test for existence of stationary solutions show that there exist a
solution for the the VAR process, however not necessarily for the
corresponding MCAR process. This fact leads us to suspect that the real part
of $\lambda^{(A)}_{3,4}$ is positive due to approximation errors. That is, the
transformation relation transforming autoregressive VAR model coefficients to
corresponding MCAR model coefficients is based on discretization, introducing
an Euler discretization error to the autoregressive (deterministic) part. Also
standard programming packages used to compute eigenvalues of large matrices
introduce approximation errors. For example, the perturbation $\rho=-0.03$ of
the blocks $A_{j}$ $(j=1,2,3,4)$ in $A$ results in eigenvalues
$\lambda^{(A)}_{3,4}=-0.00032\pm 0.0062i$, where the rest of the eigenvalues
have negative real part as well. It is also likely that additional fine-tuning
of the seasonality functions, $\bm{\Lambda}(t)$, and yearly volatility
functions, $\bm{\sigma}(t)$, give an MCAR model satisfying the existence
condition, regardless of approximation errors.

Finally, the model error coefficients, $\beta$, are calculated based on the
derivations in Section IV.3.1. As shown in Eq. (35), the coefficients are
dependent on the scale parameters, $\delta_{1}^{\mathcal{E}}$ and
$\delta_{2}^{\mathcal{E}}$, of the NIG distributions $\mathcal{E}_{1}(t)$ and
$\mathcal{E}_{2}(t)$, as well as the crosscorrelation matrix, $\Sigma$, of
$\bm{\mathcal{E}}(t)$. The formula in Eq. (35) is derived by requiring that
the scale parameters of the idiosyncratic error distributions, $\Delta
L_{1}(t)$ and $\Delta L_{2}(t)$ are equal to a constant $C_{\delta}$. As
$\delta_{1}^{\mathcal{E}}$, $\delta_{2}^{\mathcal{E}}$ and $\Sigma$ are known
parameters (see the procedure of fitting the autoregressive model
coefficients), the task of solving the system in Eq. (35) reduces to finding a
constant $C_{\delta}$ satisfying all given restrictions. An analytical
solution is not found, and a standard numerical solver has to be used to find
an optimal value of $C_{\delta}$.

The empirical covariance matrix of $\bm{\mathcal{E}}(t)$ is given by

$\displaystyle\hat{\Sigma}=\begin{bmatrix}1.018&-0.02238\\\
-0.02238&1.006\end{bmatrix}.$

By definition, the distributions of $\mathcal{E}_{1}(t)$ and
$\mathcal{E}_{2}(t)$ are nearly independent, meaning that $\beta$ is
approximately equal to the identity matrix (see Eq. (32)). To confirm this
heuristic result, and to give an illustrative example of the derivations in
Section IV.3.1, the system in Eq. (35) is solved, and the results presented.
An appropriate solution is found for the case when

$\displaystyle\sqrt{\Sigma_{11}-\beta_{1,2}^{2}}\beta_{2,1}+\sqrt{\Sigma_{22}-\beta_{2,1}^{2}}\beta_{1,2}=\Sigma_{12},$
(38)

and

$\displaystyle\beta_{1,2}=\frac{1}{2C_{\delta}}\left(\sqrt{2\Sigma_{11}C_{\delta}^{2}-(\delta_{1}^{\mathcal{E}})^{2}}-\delta_{1}^{\mathcal{E}}\right)\quad\text{and}\quad\beta_{2,1}\in\frac{1}{2C_{\delta}}\left(\sqrt{2\Sigma_{22}C_{\delta}^{2}-(\delta_{2}^{\mathcal{E}})^{2}}-\delta_{2}^{\mathcal{E}}\right).$

The estimated constant $\hat{C}_{\delta}$ with corresponding model error
coefficients, $\hat{\beta}$, and final restriction values (see Eq. (36)) are
given in Table 10. Also the relative error rates between
$\delta_{1,2}^{\mathcal{E}}$ and $\hat{\delta}_{1,2}^{\mathcal{E}}$ are given,
see Eq. (34). The empirically computed results confirm that the linear
relationship between the distributions of $\bm{\sigma}(t)$-scaled model
residuals of stratospheric temperature and U wind are approximately
independent.

Table 10: Fitted idiosyncratic error distribution scale parameter, $\hat{C}_{\delta}$, with computed model error coefficients, $\hat{\beta}$, corresponding restriction values, and an error result.  $\hat{C}_{\delta}$ | $\hat{\beta}$ | $\delta_{1}^{\mathcal{E}}/\sqrt{2\Sigma_{11}},\delta_{2}^{\mathcal{E}}/\sqrt{2\Sigma_{22}}$ | $\hat{\beta}_{1,2}^{2},\hat{\beta}_{2,1}^{2}$ | $\lvert\hat{\delta}_{1}^{\mathcal{E}}-\delta_{1}^{\mathcal{E}}\rvert/\delta_{1}^{\mathcal{E}},\lvert\hat{\delta}_{2}^{\mathcal{E}}-\delta_{2}^{\mathcal{E}}\rvert/\delta_{2}^{\mathcal{E}}$
---|---|---|---|---
$1.709$ | $\begin{bmatrix}1.009&0.01174\\\ -0.03310&1.002\end{bmatrix}$ | $1.194,1.247$ | $0.0001382,0.001096$ | $0.02358,1.255\cdot 10^{-16}$

##  V Conclusions

A transformation relation giving a direct relationship between discrete time
VARMA processes and continuous time MCARMA processes is derived. The
transformation relation is a potential method in future model estimation of
MCARMA processes. Convergence results are given as a tool to substantiate
validity of the transformation. A demonstration of applying the transformation
relation in model estimation of an MCAR model is given through a case study.

This work leans on derivations in [33], concluding that VARMA processes are
the proper discrete time analogue of MCARMA processes. Future work could
investigate if the VARMA process remains such proper analogue to the MCARMA
process when the framework is extended to include additive seasonality and
heteroscedasticity. Further, to develop methods for estimating model error
coefficients for an MCAR model with more than two driving Lévy processes is an
important topic. In this work, the transformation relation is demonstrated by
fitting an MCAR model to data. A next step would be to investigate the
possibility of estimating driving processes from MCARMA models, which includes
estimating lagged versions of Lévy processes.

##  Appendix A Proof Proposition 1

In this proof, we will see that $Q_{i}^{(l)}$ corresponds to an SDE number
(see Table 1). The statement
$Q_{i}^{(l)}\in\\{\mathcal{C}^{S},\mathcal{C}^{R},\mathcal{C}^{AR}\\}$ refer
to that the SDE with SDE number $Q_{i}^{(l)}$ is in one of the collections
$\\{\mathcal{C}^{S},\mathcal{C}^{R},\mathcal{C}^{AR}\\}$. See Eq. (9) for
definition of the recursive parameter $Q_{i}^{(l)}$. The result is derived in
four stages:

  1. 1.

The case $l=p$: Consider first $Q_{1}^{(p)}$ (i=1). We see that

$\displaystyle
Q_{1}^{(p)}=Q_{2}^{(p)}+d\in\mathcal{C}^{AR}\quad\leftrightarrow\quad
Q_{2}^{(p)}\in\mathcal{C}^{R}.$

This means that, by Eq. (5), Eq. (8) and Table 1, we might write

$\displaystyle
X_{Q_{2}^{(p)}+d}(t)dt=dX_{Q_{2}^{(p)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(1)}dL_{r}(t).$

Similarly,

$\displaystyle
Q_{2}^{(p)}=Q_{3}^{(p)}+d\in\mathcal{C}^{R}\quad\leftrightarrow\quad
Q_{3}^{(p)}\in\mathcal{C}^{R},$

such that, by Eq. (5), Eq. (8) and Table 1, we might write

$\displaystyle
X_{Q_{3}^{(p)}+d}(t)dt=dX_{Q_{3}^{(p)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(2)}dL_{r}(t).$

Continuing in a similar fashion for SDE numbers in the collection
$\mathcal{C}^{R}$, we finally hit

$\displaystyle
X_{Q_{p-1}^{(p)}+d}(t)dt=dX_{Q_{p-1}^{(p)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-2)}dL_{r}(t),$

and since

$\displaystyle
Q_{p-1}^{(p)}=Q_{p}^{(p)}+d\in\mathcal{C}^{R}\quad\leftrightarrow\quad
Q_{p}^{(p)}\in\mathcal{C}^{S},$

this finally gives

$\displaystyle
X_{Q_{p}^{(p)}+d}(t)dt=dX_{Q_{p}^{(p)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-1)}dL_{r}(t),$

by Eq. (5), Eq. (8) and Table 1. This proves that

$\displaystyle
X_{Q_{i}^{(p)}+d}(t)dt=dX_{Q_{i}^{(p)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-p+i-1)}dL_{r}(t),$

holds for $2\leq i\leq p$.

  2. 2.

The case $l=p-1$: Consider $Q_{1}^{(p-1)}$ and $Q_{2}^{(p-1)}$ ($i=1$ and
$i=2$ respectively). We see that

$\displaystyle Q_{1}^{(p-1)}\in\mathcal{C}^{R}\quad\leftrightarrow\quad
Q_{1}^{(p-1)}+d\in\mathcal{C}^{AR},$

and that

$\displaystyle
Q_{1}^{(p-1)}=Q_{2}^{(p-1)}+d\in\mathcal{C}^{R}\quad\leftrightarrow\quad
Q_{2}^{(p-1)}\in\mathcal{C}^{R}.$

By Eq. (5), Eq. (8) and Table 1 we thus might write

$\displaystyle
X_{Q_{1}^{(p-1)}+d}(t)dt=dX_{Q_{1}^{(p-1)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(1)}dL_{r}(t).$

and

$\displaystyle
X_{Q_{2}^{(p-1)}+d}(t)dt=dX_{Q_{2}^{(p-1)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(2)}dL_{r}(t).$

Continuing in a similar fashion for SDE numbers in the collection
$\mathcal{C}^{R}$, we finally hit

$\displaystyle
X_{Q_{p-2}^{(p-1)}+d}(t)dt=dX_{Q_{p-2}^{(p-1)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-2)}dL_{r}(t),$

and since

$\displaystyle
Q_{p-2}^{(p-1)}=Q_{p-1}^{(p-1)}+d\in\mathcal{C}^{R}\quad\leftrightarrow\quad
Q_{p-1}^{(p-1)}\in\mathcal{C}^{S},$

this finally gives

$\displaystyle
X_{Q_{p-1}^{(p-1)}+d}(t)dt=dX_{Q_{p-1}^{(p-1)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-1)}dL_{r}(t),$

by Eq. (5), Eq. (8) and Table 1. This proves that

$\displaystyle
X_{Q_{i}^{(p-1)}+d}(t)dt=dX_{Q_{i}^{(p-1)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-(p-1)+i-1)}dL_{r}(t),$

holds for $1\leq i\leq p-1$.

  3. 3.

The cases $l\in\\{2,\ldots,p-2\\}$: Continuing in a similar fashion as in
points 1. and 2. for $Q_{i}^{(l)}$, $l\in\\{2,\ldots,p-2\\}$, with
$Q_{1}^{(l)}$ as the starting point for each $l$, we realize that

$\displaystyle
X_{Q_{i}^{(l)}+d}(t)dt=dX_{Q_{i}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}dL_{r}(t),$
(39)

holds for for all $l\in\\{2,\ldots,p\\}$ and $1\leq i\leq l$ ($i>1$ when
$l=p$).

  4. 4.

The case $l=1$: The only valid case for $l=1$ is when $i=1$. By Eq. (5), Eq.
(8) and Table 1 this trivially gives

$\displaystyle
X_{Q_{1}^{(1)}+d}(t)dt=dX_{Q_{1}^{(1)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-1)}dL_{r}(t),$

which satisfy Eq. (39) as well.

This concludes the proof.

##  Appendix B Proof Lemma 1

The definition of the recursive parameter, $Q_{i}^{(l)}$, is used throughout
this proof, see Eq. 9. The proof is performed in four stages:

  1. 1.

The backwards recursive procedure from $x_{Q_{1}^{(l)}}$ (SDEs in
$\mathcal{C}^{R}$) to the solution SDEs, $x_{Q_{l}^{(l)}}(t)=x_{k}(t)$, in
$\mathcal{C}^{S}$ goes as follows. We know from Eq. (13) and Proposition 1
that

$\displaystyle
x_{Q_{i}^{(l)}+d}(t)=\frac{1}{h}\left(x_{Q_{i}^{(l)}}(t+h)-x_{Q_{i}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}\Delta
L_{r}(t)\right),$ (40)

for $1\leq l\leq p$ and $1\leq i\leq l$ ($i>1$ when $l=p$). Thus, for the case
$1\leq l<p$, the backwards recursive procedure starting point ($i=1$) is given
by

$\displaystyle
x_{Q_{1}^{(l)}+d}(t)=\frac{1}{h}\left(x_{Q_{1}^{(l)}}(t+h)-x_{Q_{1}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l)}\Delta
L_{r}(t)\right).$ (41)

By definition of $Q_{i}^{(l)}$ we have that
$x_{Q_{1}^{(l)}}(t)=x_{Q_{2}^{(l)}+d}(t)$, and by Eq. (40) we might write

$\displaystyle
x_{Q_{2}^{(l)}+d}(t)=\frac{1}{h}\left(x_{Q_{2}^{(l)}}(t+h)-x_{Q_{2}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+1)}\Delta
L_{r}(t)\right).$ (42)

That is, combining Eq. (41) and (42), we find

$\displaystyle\begin{split}x_{Q_{1}^{(l)}+d}=&\frac{1}{h}\Bigg{\\{}\frac{1}{h}\left(x_{Q_{2}^{(l)}}(t+2h)-x_{Q_{2}^{(l)}}(t+h)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+1)}\Delta
L_{r}(t+h)\right)\\\
&\qquad-\frac{1}{h}\left(x_{Q_{2}^{(l)}}(t+h)-x_{Q_{2}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+1)}\Delta
L_{r}(t)\right)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l)}\Delta
L_{r}(t)\Bigg{\\}}\\\
=&\frac{1}{h^{2}}\Bigg{(}x_{Q_{2}^{(l)}}(t+2h)-2x_{Q_{2}^{(l)}}(t+h)+x_{Q_{2}^{(l)}}(t)\Bigg{)}\\\
&-\frac{1}{h^{2}}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+1)}\left(\Delta
L_{r}(t+h)-\Delta
L_{r}(t)\right)-\frac{1}{h}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l)}\Delta
L_{r}(t).\end{split}$ (43)

Further, by definition of $Q_{i}^{(l)}$ we have that
$x_{Q_{2}^{(l)}}(t)=x_{Q_{3}^{(l)}+d}(t)$, and by Eq. (40) we might write

$\displaystyle
x_{Q_{3}^{(l)}+d}(t)=\frac{1}{h}\left(x_{Q_{3}^{(l)}}(t+h)-x_{Q_{3}^{(l)}}(t)-\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+2)}\Delta
L_{r}(t)\right).$ (44)

That is, $x_{Q_{1}^{(l)}+d}(t)$ in Eq. (43) might be rewritten as a function
of $x_{Q_{3}^{(l)}}(t+3h),\ldots,x_{Q_{3}^{(l)}}(t),$ $\Delta
L_{r}(t+2h),\Delta L_{r}(t+h),\Delta L_{r}(t)$ using Eq. (44). Continuing in a
similar fashion for $x_{Q_{3}^{(l)}}=x_{Q_{4}^{(l)}+d}$, and so on, to
recursively rewrite $x_{Q_{1}^{(1)}+d}$, we finally reach a right hand side of
Eq. (43) which is given by
$x_{Q_{l}^{(l)}}(t+lh),\ldots,x_{Q_{l}^{(l)}}(t),\Delta
L_{r}(t+(l-1)h),\ldots,\Delta L_{r}(t)$.

  2. 2.

The backwards recursive procedure in point 1. is a analogous to the
$1$-dimensional AR/CAR transformation relation formula, see [6]. From [6]
(Lemma 10.2) we know that when a discrete variable $x_{\tilde{q}+1}(t)$
($\tilde{q}\in\mathbb{N}$) satisfies

$\displaystyle x_{\tilde{q}+1}(t)=x_{\tilde{q}}(t+1)-x_{\tilde{q}}(t),\quad
1\leq\tilde{q}\leq p-1,$ (45)

it holds that

$\displaystyle
x_{\tilde{q}+1}(t)=\sum_{n=0}^{\tilde{q}}(-1)^{n}b_{n}^{\tilde{q}}x_{1}(t+\tilde{q}-n),$
(46)

where the coefficients $b_{n}^{\tilde{q}}$ ($n\in\mathbb{N}$) are defined
recursively as

$\displaystyle b_{n}^{\tilde{q}}=b_{n-1}^{\tilde{q}-1},\quad 1\leq n\leq
p-1,\tilde{q}\geq 2,$ (47)

and $b_{0}^{\tilde{q}}=b_{\tilde{q}}^{\tilde{q}}=1$ for $0\leq\tilde{q}\leq
p$. As shown in point 1., our goal is to move recursively backwards from
$x_{Q_{1}^{(l)}}(t)$ to $x_{Q_{l}^{(l)}}(t)=x_{k}(t)$
($k\in\\{1,\ldots,d\\}$), just as the formula in Eq. (46) achieves for the
$1$-dimensional AR/CAR case (however from $x_{\tilde{q}}$ to $x_{1}$ in that
specific case). In our multivariate case, when $1\leq l<p$, the discrete
variable $x_{Q_{i}^{(l)}+d}(t)$ satisfies Eq. (40) for all $1\leq i\leq l$.
Notice that the deterministic part of the right hand side in Eq. (40) is
equivalent to the relationship of the 1-dimensional case in Eq. (45). By this
observation, we see that the formula in Eq. (46) might be modified to hold for
the deterministic part of the $d$-dimensional VARMA/MCARMA case. This is the
goal of point 3.

  3. 3.

The formula in Eq. (46) will be modified according to the deterministic part
of the right hand side, $1/h(x_{Q_{i}^{l}}(t+h)-x_{Q_{i}^{(l)}}(t))$, of Eq.
(40). Consider the four arguments that follows.

    1. (a)

As mentioned in point 2., the final goal of the current multivariate case is
to reach $x_{Q_{l}^{(l)}}(t)=x_{k}(t)$ ($k\in\\{1,\ldots,d\\}$), meaning that
the right hand side of Eq. (46) should be stated as

$\displaystyle\sum_{n=0}^{\tilde{q}}(-1)^{n}b_{n}^{\tilde{q}}x_{k}(t+\tilde{q}-n).$
(48)

    2. (b)

In the $1$-dimensional case, the backwards recursive procedure moves from
variable $x_{\tilde{q}}(t)$ to $x_{1}(t)$ by inserting Eq. (45) repeatedly (as
shown in point 1. for the current multivariate case). That is, the variable
index number is reduced by $\tilde{q}-1$, and the formula requires
$\tilde{q}+1$ terms ($n$ running from $0$ to $\tilde{q}$) in order to do that
procedure. In the $d$-dimensional case, the backwards recursive procedure
moves from variable $x_{Q_{1}^{(l)}}$ to $x_{Q_{l}^{(l)}}$ in the same way as
for the $1$-dimensional case, with a variable index number increased by $l-1$
in the current case (notice that the indexes altered in the 1- and
d-dimensional recursive procedures are defined in opposite order). That is,
the formula in Eq. (48) requires $l+1$ terms ($n$ running from $0$ to $l$),
and the updated formula becomes

$\displaystyle\sum_{n=0}^{l}(-1)^{n}b_{n}^{l}x_{k}(t+l-n).$ (49)

    3. (c)

In the current multivariate case, an Euler scheme step size $h$ is considered
rather than step size equal to $1$, as in the 1-dimensional case. As a
consequence, for each step in the backwards recursive procedure (see point
1.), the substituted expression adds an extra factor $1/h$ to the formula in
Eq. (49). The starting point of the recursive procedure (see Eq. (41)) gives
the first factor $1/h$, and each of the $l-1$ recursive steps adds one
additional factor. That is, the right hand side of Eq. (49) has to be modified
as

$\displaystyle\frac{1}{h^{l}}\sum_{n=0}^{l}(-1)^{n}b_{n}^{l}x_{k}(t+l-n).$
(50)

    4. (d)

An additional consequence of setting the Euler scheme step size to $h$ is
that, when the Euler scheme has moved $r$ steps forward, the length of each
step is equal to $rh$. That is, the following modification of the right hand
side of Eq. (46) has to be done

$\displaystyle\frac{1}{h^{l}}\sum_{n=0}^{l}(-1)^{n}b_{n}^{l}x_{k}(t+(l-n)h).$
(51)

This concludes the deterministic part of the formula we want to prove.

  4. 4.

The result of the backwards recursive procedure for the stochastic part of the
right hand side of Eq. (40) has to be added to the formula in Eq. (51). The
illustration of the recursive procedure in point 1. shows that, for each
recursive step, incremental Lévy-terms are added. Continuing the explicit
calculations in point 1. by one recursive step adds the stochastic term

$\displaystyle-\frac{1}{h^{3}}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-j+2)}\left(\Delta
L_{r}(t+2h)-2\Delta L_{r}(t+h)+\Delta L_{r}(t)\right),$

to the recursive procedure of $x_{Q_{1}^{(l)}+d}(t)$. It is evident that the
recursively added incremental Lévy-terms follow a similar pattern as the
recursively added deterministic terms, however with a delay of one step. Also
note that each of the added Lévy-terms remains in the final formula, rather
than being updated at each recursive step as for the deterministic part. By
these remarks, the formula for the added incremental Lévy-terms is found
through the following arguments:

    1. (a)

As the recursively added incremental Lévy-terms follow a similar pattern as
the recursively added deterministic terms, we use the formula in Eq. (51) as a
starting point. As seen in point 1., the variable in the stochastic case is
$\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}\Delta L_{r}(t)$ rather than
$x_{k}(t)$. That is, the first step towards the recursive formula of the
stochastic part on the right hand side of Eq. (41) is

$\displaystyle\frac{1}{h^{l}}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}\sum_{v=0}^{l}(-1)^{v}b_{v}^{l}\Delta
L_{r}(t+(l-v)h).$ (52)

    2. (b)

The incremental Lévy-terms are not regressed on in the recursive procedure. As
seen in point 1., each added Lévy-term originate from the recursive steps on
$1/h(x_{Q_{i}^{(l)}}(t+h)-x_{Q_{i}^{(l)}}(t))$ ($1\leq i<l$). This is why each
recursively added stochastic term remains in the final expression. That is,
each of the $l-1$ added incremental Lévy-terms during the recursive procedure
are given by the recursive formula in Eq. (52), and the stochastic recursive
procedure formula has to be modified as

$\displaystyle\sum_{w=0}^{l-1}\frac{1}{h^{l}}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h).$ (53)

    3. (c)

In point 1. we see that each of the added incremental Lévy-terms brings one
extra factor $1/h$ for each recursive step. The index $w$ in Eq. (53) might
account for this by modifying the formula as

$\displaystyle\sum_{w=0}^{l-1}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta_{(k-1)m+r}^{(p-l+i-1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h).$ (54)

    4. (d)

Finally, from the calculations in point 1. we see that the top index of
$\beta$ increases from $p-l$ by one for each added term. The index $w$ in Eq.
(54) might be used to account for this as well, and we modify the formula as

$\displaystyle\sum_{w=0}^{l-1}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta_{r+(k-1)m}^{(p-l+w)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h).$ (55)

Add the deterministic and stochastic recursive formulas in Eq. (51) and (55),
respectively, according to the formula in Eq. (41) with regards to signs, to
find the backwards recursive formula of $x_{Q_{1}^{(l)}+d}(t)$ as stated in
the lemma, and the proof is complete.

##  Appendix C Proof Theorem 1

Combine Eq. (8) and the recursive parameter in Eq. (9) to rewrite the SDEs in
$\mathcal{C}^{AR}$ (Eq. (6)) as the piecewise constant process

$\displaystyle
x_{Q_{1}^{(p)}}(t+h)-x_{Q_{1}^{(p)}}(t)=-h\sum_{l=1}^{p}\sum_{s=1}^{d}\alpha^{(p-l+1)}x_{Q_{1}^{(l)}}(t)+\sum_{r=1}^{m}\beta^{(0)}\Delta
L_{r}(t),$ (56)

using the Euler scheme (see [29] and [37]). Here, we have defined

$\displaystyle\alpha^{(p-l+1)}\triangleq\alpha_{(k-1)d+s}^{(p-l+1)}\quad\text{and}\quad\beta^{(\kappa)}\triangleq\beta^{(\kappa)}_{(k-1)m+r}$

for notational convenience. Note that we sum
$Q_{1}^{(l)}\triangleq(Q_{1}^{(l)}\mid s)$ over $s$ as well. Extract the terms
from the matrix-vector product $[-A_{p}\cdots-A_{1}]\bm{X}(t)dt$ corresponding
to equation block $1$ (variables $x_{Q_{1}^{(1)}}=x_{k}$), such that we can
write $Q_{1}^{(l)}$ as $Q_{1}^{(l-1)}+d$. That is, rewrite Eq. (56) as

$\displaystyle x_{Q_{1}^{(p-1)}+d}(t+h)-x_{Q_{1}^{(p-1)}+d}(t)=$
$\displaystyle-h\sum_{s=1}^{d}\alpha^{(p)}x_{s}(t)+\sum_{r=1}^{m}\beta^{(0)}\Delta
L_{r}(t)-h\sum_{l=2}^{p}\sum_{s=1}^{d}\alpha^{(p-l+1)}x_{Q_{1}^{(l-1)}+d}(t).$
(57)

By Lemma 1, the left hand side of Eq. (57) can be written as

$\displaystyle\frac{1}{h^{p-1}}\sum_{n=0}^{p-1}(-1)^{n}b_{n}^{p-1}\left(x_{k}(t+(p-n)h)-x_{k}(t+(p-1-n)h)\right)$
$\displaystyle-\sum_{w=0}^{p-2}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta^{(w+1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\left(\Delta
L_{r}(t+(w-v+1)h)-\Delta L_{r}(t+(w-v)h)\right),$

and the right hand side as

$\displaystyle-h\sum_{s=1}^{d}\alpha^{(p)}x_{s}(t)+\sum_{r=1}^{m}\beta^{(0)}\Delta
L_{r}(t)$
$\displaystyle-h\sum_{l=2}^{p}\sum_{s=1}^{d}\alpha^{(p-l+1)}\Bigg{(}\frac{1}{h^{l-1}}\sum_{n=0}^{l-1}(-1)^{n}b_{n}^{l-1}x_{s}(t+(l-1-n)h)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\sum_{w=0}^{l-2}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta^{(p-l+w+1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h)\Bigg{)}.$

Remember that $k=s$ in the subscript of $\beta^{(p-l+w+1)}$. That is, Eq. (57)
can be rearranged as

$\displaystyle\frac{1}{h^{p-1}}\left(x_{k}(t+ph)-x_{k}(t+(p-1)h)\right)=-h\sum_{s=1}^{d}\alpha^{(p)}x_{s}(t)+\sum_{r=1}^{m}\beta^{(0)}\Delta
L_{r}(t)$
$\displaystyle-h\sum_{l=2}^{p}\sum_{s=1}^{d}\alpha^{(p-l+1)}\Bigg{(}\frac{1}{h^{l-1}}\sum_{n=0}^{l-1}(-1)^{n}b_{n}^{l-1}x_{s}(t+(l-1-n)h)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\sum_{w=0}^{l-2}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta^{(p-l+w+1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\Delta
L_{r}(t+(w-v)h)\Bigg{)}$
$\displaystyle-\frac{1}{h^{p-1}}\sum_{n=1}^{p-1}(-1)^{n}b_{n}^{p-1}\left(x_{k}(t+(p-n)h)-x_{k}(t+(p-1-n)h)\right)$
$\displaystyle+\sum_{w=0}^{p-2}\frac{1}{h^{w+1}}\sum_{r=1}^{m}\beta^{(w+1)}\sum_{v=0}^{w}(-1)^{v}b_{v}^{w}\left(\Delta
L_{r}(t+(w-v+1)h)-\Delta L_{r}(t+(w-v)h)\right).$

where the highest lagged variable, meaning $x_{k}(t+ph)$, is extracted, and
$b_{0}^{p-1}=1$ is inserted (see Eq. (46)). Solve for $x_{k}(t+ph)$ to
complete the proof.

##  Appendix D Computations to find MCAR model coefficients of deterministic
part

Compare Eq. (29) and (30) to see that

$\displaystyle\begin{cases}-\alpha_{1,1}^{(1)}-1=\phi_{11}^{(1)}\\\
-\alpha_{1,2}^{(1)}=\phi_{12}^{(1)}\end{cases}\quad\Rightarrow\quad\begin{cases}\alpha_{1,1}^{(1)}=-\phi_{11}^{(1)}-1\\\
\alpha_{1,2}^{(1)}=-\phi_{12}^{(1)}\end{cases}$
$\displaystyle\begin{cases}-\alpha_{1,1}^{(2)}+3\alpha_{1,1}^{(1)}+4=\phi_{11}^{(2)}\\\
-\alpha_{1,2}^{(2)}+3\alpha_{1,2}^{(1)}=\phi_{12}^{(2)}\end{cases}\quad\Rightarrow\quad\begin{cases}\alpha_{1,1}^{(2)}=-3\phi_{11}^{(1)}-\phi_{11}^{(2)}+1\\\
\alpha_{1,2}^{(2)}=-3\phi_{12}^{(1)}-\phi_{12}^{(2)}\end{cases}$
$\displaystyle\begin{cases}-\alpha_{1,1}^{(3)}+2\alpha_{1,1}^{(2)}-3\alpha_{1,1}^{(1)}-6=\phi_{11}^{(3)}\\\
-\alpha_{1,2}^{(3)}+2\alpha_{1,2}^{(2)}-3\alpha_{1,2}^{(1)}=\phi_{12}^{(3)}\\\
\end{cases}\quad\Rightarrow\quad\begin{cases}\alpha_{1,1}^{(3)}=-3\phi_{11}^{(1)}-2\phi_{11}^{(2)}-\phi_{11}^{(3)}-1\\\
\alpha_{1,2}^{(3)}=-3\phi_{12}^{(1)}-2\phi_{12}^{(2)}-\phi_{11}^{(3)}\\\
\end{cases}$
$\displaystyle\begin{cases}-\alpha_{1,1}^{(4)}+\alpha_{1,1}^{(3)}-\alpha_{1,1}^{(2)}+\alpha_{1,1}^{(1)}+4=\phi_{11}^{(4)}\\\
-\alpha_{1,2}^{(4)}+\alpha_{1,2}^{(3)}-\alpha_{1,2}^{(2)}+\alpha_{1,2}^{(1)}=\phi_{12}^{(4)}\end{cases}\quad\Rightarrow\quad\begin{cases}\alpha_{1,1}^{(4)}=-\phi_{11}^{(1)}-\phi_{11}^{(2)}-\phi_{11}^{(3)}-\phi_{11}^{(4)}+1\\\
\alpha_{1,2}^{(4)}=-\phi_{12}^{(1)}-\phi_{12}^{(2)}-\phi_{12}^{(3)}-\phi_{12}^{(4)}.\end{cases}$

In the above solutions, change subscripts $1,1$ and $11$ with $2,2$ and $22$,
respectively, to find the solutions of
$\alpha_{2,2}^{(1)},\ldots,\alpha_{2,2}^{(4)}$. Further, change subscripts
$1,2$ and $12$ with $2,1$ and $21$, respectively, to find the solutions of
$\alpha_{2,1}^{(1)},\ldots,\alpha_{2,1}^{(4)}$. Now, by the definition
$\alpha_{k,i_{s}}^{(j)}\triangleq\alpha_{(k-1)d+i_{s}}^{(p-l+1)}$, we find the
matrices $A_{p},\ldots,A_{1}$ (see Eq. (7)).

##  Appendix E Computations to find MCAR model coefficients of stochastic part

The system of equations in Eq. (33) and (34) is given by

$\displaystyle\begin{cases}\delta_{1}^{\mathcal{E}}=\lvert\beta_{1,1}\rvert
C_{\delta}+\lvert\beta_{1,2}\rvert C_{\delta}\\\
\delta_{2}^{\mathcal{E}}=\lvert\beta_{2,1}\rvert
C_{\delta}+\lvert\beta_{2,2}\rvert
C_{\delta}\end{cases}\quad\text{and}\quad\begin{bmatrix}\Sigma_{11}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\end{bmatrix}=\begin{bmatrix}\beta_{1,1}^{2}+\beta_{1,2}^{2}&\beta_{1,1}\beta_{2,1}+\beta_{1,2}\beta_{2,2}\\\
&\beta_{2,1}^{2}+\beta_{2,2}^{2}\end{bmatrix}.$ (58)

Solve the expressions for $\delta_{1}^{\mathcal{E}}$,
$\delta_{2}^{\mathcal{E}}$, $\Sigma_{11}$ and $\Sigma_{22}$ in Eq. (58) with
respect to $\beta_{1,1}$ and $\beta_{2,2}$, and equate them to find

$\displaystyle\begin{cases}\pm(\delta_{1}^{\mathcal{E}}/C_{\delta}-\lvert\beta_{1,2}\rvert)=\pm\sqrt{\Sigma_{11}-\beta_{1,2}^{2}}\\\
\pm(\delta_{2}^{\mathcal{E}}/C_{\delta}-\lvert\beta_{2,1}\rvert)=\pm\sqrt{\Sigma_{22}-\beta_{2,1}^{2}}.\end{cases}$
(59)

Solve the system in Eq. (59) with respect to $\beta_{1,2}$ and $\beta_{2,1}$.
Both equations are solved in a similar way, and we therefore continue the
calculations with

$\displaystyle\\{\beta,\sigma,\delta,C\\}\in\\{(\beta_{1,2},\beta_{2,1}),(\Sigma_{11},\Sigma_{22}),(\delta_{1}^{\mathcal{E}},\delta_{2}^{\mathcal{E}}),(C_{\delta},C_{\delta})\\}.$

That is,

$\displaystyle\begin{split}(\delta/C-\lvert\beta\rvert)^{2}&=\sigma-\beta^{2}\\\
\rightarrow\quad\beta^{2}-\lvert\beta\rvert\delta/C+\frac{1}{2}(\delta^{2}/C^{2}-\sigma)&=0.\end{split}$
(60)

For $\beta$ positive and $\beta$ negative, Eq. (60) gives the following
respective solutions by the quadratic equation

$\displaystyle\beta=\frac{1}{2C}\left(\delta\pm\sqrt{2\sigma
C^{2}-\delta^{2}}\right)\quad\text{and}\quad\beta=-\frac{1}{2C}\left(\delta\mp\sqrt{2\sigma
C^{2}-\delta^{2}}\right).$

Finally, combine the expression of $\Sigma_{12}$ in Eq. (58) with the right
hand side of Eq. (59) to find

$\displaystyle\pm\sqrt{\Sigma_{11}-\beta_{1,2}^{2}}\beta_{2,1}\pm\sqrt{\Sigma_{22}-\beta_{2,1}^{2}}\beta_{1,2}=\Sigma_{12}.$

Supporting information. Additional information for this article is available
online. Raw data used in this work are available from link in Appendix S1.

##  Appendix S1 Data availability

https://apps.ecmwf.int/datasets/

References

## References

  * Applebaum [2004] Applebaum, D.,  _Lévy processes and stochastic calculus_ , Vol. 93 (Cambridge University Press, 2004).
  * Asmussen and Rosiński [2001] Asmussen, S.and Rosiński, J., “Approximations of small jumps of lévy processes with a view towards simulation,” Journal of Applied Probability 38, 482––493 (2001).
  * Barndorff-Nielsen [1997] Barndorff-Nielsen, O. E., “Normal inverse Gaussian distributions and stochastic volatility modelling,” Scandinavian Journal of Statistics 24, 1–13 (1997).
  * Barndorff-Nielsen, Mikosch, and Resnick [2001] Barndorff-Nielsen, O. E., Mikosch, T., and Resnick, S. I.,  _Lévy Processes: Theory and Applications_ (Birkhäuser Boston: Imprint: Birkhäuser, 2001).
  * Barndorff-Nielsen and Schmiegel [2007] Barndorff-Nielsen, O. E.and Schmiegel, J., “Ambit processes; with applications to turbulence and tumour growth,” in  _Stochastic Analysis and Applications_ , edited by F. E. Benth, G. Di Nunno, T. Lindstrøm, B. Øksendal, and T. Zhang (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 93–124.
  * Benth, Šaltytė Benth, and Koekebakker [2008] Benth, F. E., Šaltytė Benth, J., and Koekebakker, S.,  _Stochastic modelling of electricity and related markets_ ((Vol. 11, Advanced series on statistical science & applied probability). Singapore: World Scientific Publishing Pte., 2008).
  * Benth, Di Nunno, and Khedher [2011] Benth, F. E., Di Nunno, G., and Khedher, A., “Robustness of option prices and their deltas in markets modelled by jump-diffusions,” Communications on Stochastic Analysis 5 (2011), 10.31390/cosa.5.2.03.
  * Benth and Taib [2013] Benth, F. E.and Taib, C. M. I. C., “On the speed towards the mean for continuous time autoregressive moving average processes with applications to energy markets,” Energy economics 40, 259–268 (2013).
  * Berrisford  _et al._ [2011] Berrisford, P., Dee, D. P., Poli, P., Brugge, R., Fielding, M., Fuentes, M., Kållberg, P. W., Kobayashi, S., Uppala, S., and Simmons, A., “The ERA-Interim archive Version 2.0,” ERA Report Series , 1–23 (2011).
  * Brockwell [2014] Brockwell, P. J., “Recent results in the theory and applications of CARMA processes,” Annals of the Institute of Statistical Mathematics 66, 647–685 (2014).
  * Brockwell and Davis [1991] Brockwell, P. J.and Davis, R. A.,  _Time Series: Theory and Methods_ (Springer New York, 1991).
  * Brockwell, Ferrazzano, and Klüppelberg [2013] Brockwell, P. J., Ferrazzano, V., and Klüppelberg, C., “High-frequency sampling and kernel estimation for continuous-time moving average processes,” Journal of Time Series Analysis 34, 385–404 (2013).
  * Brockwell and Lindner [2019] Brockwell, P. J.and Lindner, A., “Sampling, embedding and inference for carma processes,” Journal of Time Series Analysis 40, 163–181 (2019).
  * Brockwell and Schlemm [2013] Brockwell, P. J.and Schlemm, E., “Parametric estimation of the driving lévy process of multivariate carma processes from discrete observations,” Journal of Multivariate Analysis 115, 217–251 (2013).
  * Broszkiewicz-Suwaj and Wyłomańska [2021] Broszkiewicz-Suwaj, E.and Wyłomańska, A., “Application of non-gaussian multidimensional autoregressive model for climate data prediction,” International Journal of Advances in Engineering Sciences and Applied Mathematics 13, 236–247 (2021).
  * Chambers and Thornton [2012] Chambers, M. J.and Thornton, M. A., “Discrete time representation of continuous time arma processes,” Econometric Theory 28, 219–238 (2012).
  * Dee  _et al._ [2011] Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge‐Sanz, B. M., Morcrette, J. J., Park, B. K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J. N., and Vitart, F., “The ERA‐Interim reanalysis: configuration and performance of the data assimilation system,” Quarterly Journal of the Royal Meteorological Society 137, 553–597 (2011).
  * Eggen  _et al._ [2022] Eggen, M. D., Dahl, K. R., Näsholm, S. P., and Mæland, S., “Stochastic modeling of stratospheric temperature,” Mathematical Geosciences (2022), 10.1007/s11004-021-09990-6.
  * Fasen-Hartmann and Mayer [2021] Fasen-Hartmann, V.and Mayer, C., “Whittle estimation for continuous-time stationary state space models with finite second moments,” Annals of the Institute of Statistical Mathematics (2021), 10.1007/s10463-021-00802-6.
  * Fasen-Hartmann and Scholz [2021] Fasen-Hartmann, V.and Scholz, M., “Factorization and discrete-time representation of multivariate carma processes,”  (2021).
  * García, Klüppelberg, and Müller [2011] García, I., Klüppelberg, C., and Müller, G., “Estimation of stable carma models with an application to electricity spot prices,” Statistical Modelling 11, 447–470 (2011).
  * Gardoń [2004] Gardoń, A., “The order of approximations for solutions of itô-type stochastic differential equations with jumps,” Stochastic Analysis and Applications 22, 679–699 (2004).
  * Gómez [2016] Gómez, V.,  _Multivariate Time Series With Linear State Space Structure_ (Springer International Publishing, 2016).
  * Gómez [2019] Gómez, V.,  _Linear Time Series with MATLAB and OCTAVE_ (Springer International Publishing: Imprint: Springer, 2019).
  * Gómez [2020] Gómez, V., “Ssmmatlab,” https://github.com/vgomezenriquez/ssmmatlab (2020).
  * Hitchcock and Simpson [2014] Hitchcock, P.and Simpson, I. R., “The downward influence of stratospheric sudden warmings,” Journal of the Atmospheric Sciences 71, 3856 – 3876 (2014).
  * James, Koreisha, and Partch [1985] James, C., Koreisha, S., and Partch, M., “A varma analysis of the causal relations among stock returns, real output, and nominal interest rates,” The Journal of Finance (New York) 40, 1375–1384 (1985).
  * Karpechko, Tummon, and Secretariat [2016] Karpechko, A., Tummon, F., and Secretariat, W., “Climate predictability in the stratosphere,” Bulletin of World Meteorological Organization (WMO) 65 (2016).
  * Kloeden and Platen [1992] Kloeden, P. E.and Platen, E.,  _Numerical solution of stochastic differential equations_ , Vol. 23 (Springer, 1992).
  * Kühn and Schilling [2019] Kühn, F.and Schilling, R. L., “Strong convergence of the euler–maruyama approximation for a class of lévy-driven sdes,” Stochastic Processes and Their Applications 129, 2654–2680 (2019).
  * Levendis [2018] Levendis, J. D.,  _Time Series Econometrics: Learning Through Replication_ (Springer International Publishing: Imprint: Springer, 2018).
  * Love, Matthews, and Janacek [2008] Love, B. S., Matthews, A. J., and Janacek, G. J., “Real-time extraction of the madden–julian oscillation using empirical mode decomposition and statistical forecasting with a varma model,” Journal of Climate 21, 5318–5335 (2008).
  * Marquardt and Stelzer [2007] Marquardt, T.and Stelzer, R., “Multivariate carma processes,” Stochastic Processes and Their Applications 117, 96–120 (2007).
  * Parlange and Katz [2000] Parlange, M. B.and Katz, R. W., “An extended version of the richardson model for simulating daily weather variables,” Journal of Applied Meteorology (1988) 39, 610–622 (2000).
  * Pham and Le Breton [1991] Pham, D. T.and Le Breton, A., “Levinson-durbin-type algorithms for continuous-time autoregressive models and applications,” Mathematics of Control, Signals, and Systems 4, 69–79 (1991).
  * Platen and Bruti-Liberati [2010] Platen, E.and Bruti-Liberati, N.,  _Numerical Solution of Stochastic Differential Equations with Jumps in Finance_ , Vol. 64 (Springer, 2010).
  * Protter and Talay [1997] Protter, P.and Talay, D., “The Euler scheme for Lévy driven stochastic differential equations,” The Annals of Probability 25, 393 – 423 (1997).
  * Scaife  _et al._ [2022] Scaife, A. A., Baldwin, M. P., Butler, A. H., Charlton-Perez, A. J., Domeisen, D. I. V., Garfinkel, C. I., Hardiman, S. C., Haynes, P., Karpechko, A. Y., Lim, E.-P., Noguchi, S., Perlwitz, J., Polvani, L., Richter, J. H., Scinocca, J., Sigmond, M., Shepherd, T. G., Son, S.-W., and Thompson, D. W. J., “Long-range prediction and the stratosphere,” Atmospheric chemistry and physics 22, 2601–2623 (2022).
  * Schlemm and Stelzer [2012a] Schlemm, E.and Stelzer, R., “Multivariate carma processes, continuous-time state space models and complete regularity of the innovations of the sampled processes,” Bernoulli: Official Journal of the Bernoulli Society for Mathematical Statistics and Probability 18, 46–63 (2012a).
  * Schlemm and Stelzer [2012b] Schlemm, E.and Stelzer, R., “Quasi maximum likelihood estimation for strongly mixing state space models and multivariate lévy-driven carma processes,” Electronic Journal of Statistics 6, 2185–2234 (2012b).
  * Söderström  _et al._ [1997] Söderström, T., Fan, H., Carlsson, B., and Mossberg, M., “Some approaches on how to use the delta operator when identifying continuous-time processes,” in  _Proceedings of the 36th IEEE Conference on Decision and Control_ (1997) pp. 890–895.
  * Thornton and Chambers [2017] Thornton, M. A.and Chambers, M. J., “Continuous time arma processes: Discrete time representation and likelihood evaluation,” Journal of Economic Dynamics and Control 79, 48–65 (2017).
  * Todorov and Tauchen [2006] Todorov, V.and Tauchen, G., “Simulation methods for lévy-driven continuous-time autoregressive moving average (carma) stochastic volatility models,” Journal of Business & Economic Statistics 24, 455–469 (2006).
  * Wei [2019] Wei, W. W. S.,  _Multivariate time series analysis and applications_ (Wiley, 2019).