Autoregressive moving average model for matrix time series

Abstract

In the paper, the autoregressive moving average model for matrix time series (MARMA) is investigated. The properties of the MARMA model are investigated by using the conditional least square estimation, the conditional maximum likelihood estimation, the projection theorem in Hilbert space and the decomposition technique of time series, which include necessary and sufficient conditions for stationarity and invertibility, model parameter estimation, model testing and model forecasting.

Keywords:

• Matrix time series
• autoregressive moving average model
• bilinear model
• statistical inference

1. Introduction

Matrix time series is a time series whose cross-sectional data are matrices, which can be found in a variety of fields such as economics, business, ecology, psychology, meteorology, biology and fMRI (Samadi, 2014). For example, consider two stocks, $A_1$ and $A_2$, as potential investment products, whose prices and volumes are selected as two analysis factors. Denote the price and volume of stock $A_k$ at time t by $P_k\left(t\right)$ and $V_k\left(t\right)$, k = 1, 2, and then a $2\times 2$-dimensional matrix time series can be constructed as follows: $\{X_t\equiv \left[\begin{array}{cc}{P}_1\left(t\right)& P_2\left(t\right)\\ V_1\left(t\right)& V_2\left(t\right)\end{array}\right],t=1,2,\dots \}.$Matrix time series has attracted a few scholars' attention and research at the beginning of the century. Walden and Serroukh (2002) studied the construction of matrix-valued filters for multi-resolution analysis of matrix time series. Samadi (2014) brought forward and investigated a p-order autoregressive model for matrix time series, which is essentially a VAR(p) model in matrix form. D. Wang et al. (2019) proposed a novel factor model $X_t=R{F}_t{C}^{\mathrm{\top }}+{ϵ}_t,t=1,2,\dots ,$where $X_t$ and $F_t$ are matrix time series. Chen et al. (2021) first proposed one-order autoregressive model for matrix time series in the bilinear form, denoted by MAR(1), (1) $X_t=A{X}_{t-1}{B}^{\mathrm{\top }}+{ϵ}_t,t=1,2,\dots ,$(1) and investigated its stationarity, causality, method of parameter estimation, and asymptotics of statistic. Wu and Hua (2022) independently proposed the p-order autoregressive model for matrix time series in the bilinear form, denoted by MAR(p), (2) $X_t=\sum _{k=1}^p{A}_k{X}_{t-k}{B}_k^{\mathrm{\top }}+{ϵ}_t,t=1,2,\dots ,$(2) and presented parameter estimation, model identification criterion and model checking. For more literature studies on matrix time series, one can refer to H. Wang and West (2009), Zhou et al. (2018), Getmanov et al. (2021) and their references.

It is widely known that the autoregressive moving average model of time series (ARMA) plays a very important role in the theory and the application of one-dimensional time series, and we will show later that a bilinear model has its unique advantages for matrix time series. In the paper, autoregressive moving average models for matrix time series (MARMA) are first proposed and investigated. Necessary and sufficient conditions for stationarity of MARMA are provided, and parameter estimations are also considered by the conditional least squares method and the conditional maximum likelihood estimation method. At last, an example is presented to show the applications of the MARMA model.

2. Preliminaries

Let $(\mathrm{\Omega },\mathcal{F},P)$ be a probability space with a σ-filtration $\{{\mathcal{F}}_t,t\in \mathbb{N}\}$ in which the second moment of each variable exists, and $\mathbb{N}=\{1,2,3,\dots \}$.

Definition 2.1

For any given positive integers m and n, an $m\times n$-dimensional matrix time series refers to (3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) where $\left\{X_{ij}\right(t),t\in \mathbb{N}\}$ is a one-dimensional time series on a probability space $(\mathrm{\Omega },\mathcal{F},P)$ for any $i=1,2,\dots ,n$ and $j=1,2,\dots ,m$.

Definition 2.2

Let $X=\left\{X\right(t),t\in \mathbb{N}\}$ be an $m\times n$-dimensional matrix time series defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ), and then its mean function follows as (4) ${\mu }_X\left(t\right)\equiv E\left[X\right(t\left)\right]={\left(E\right[X_{ij}\left(t\right)\left]\right)}_{m\times n},t\in \mathbb{N}.$(4) Additionally, its autocovariance function follows as (5) ${\mathrm{\Gamma }}_X(t,s)\equiv {\mathrm{\Gamma }}_{\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right)}(t,s)={\left({\sigma }_{ij,k\ell }\right(t,s\left)\right)}_{mn\times mn},$(5) where ${\sigma }_{ij,k\ell }(t,s)=\mathrm{c}\mathrm{o}\mathrm{v}\left(X_{ij}\right(t),X_{k\ell }(s\left)\right)$, $i,k=1,2,\dots ,n$ and $j,\ell =1,2,\dots ,n$; $t,s\in \mathbb{N}$, and $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)$ is the vectorization of $X\left(t\right)$ by columns, that is, (6) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)={\left[X_{11}\right(t),X_{21}(t),\dots ,X_{m1}(t),X_{12}(t),X_{22}(t),\dots ,X_{mn}(t\left)\right]}^{\mathrm{\top }}.$(6)

Stationarity and matrix white noise play a very important role on time series analysis. Thus, we will introduce the concept of stationary matrix time series and matrix white noise in the following.

Definition 2.3

Let $\left\{X\right(t),t\in \mathbb{N}\}$ be a matrix time series defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ) and $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)$ be the vectorization of $X\left(t\right)$ defined by (6(6) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)={\left[X_{11}\right(t),X_{21}(t),\dots ,X_{m1}(t),X_{12}(t),X_{22}(t),\dots ,X_{mn}(t\left)\right]}^{\mathrm{\top }}.$(6) ). Then $\left\{X\right(t),t\in \mathbb{N}\}$ is a stationary matrix time series if and only if $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$ is stationary.

Definition 2.4

For any given positive integers m and n, denote an $m\times n$-dimensional matrix time series $ϵ=\left\{\right({ϵ}_{ij}\left(t\right){)}_{m\times n},t\in \mathbb{N}\}$, and then ε is called an $m\times n$-dimensional matrix white noise, if it satisfies the following conditions.

1. Its mean function ${\mu }_{ϵ}\left(t\right)=O_{m\times n}$ for all $t\in \mathbb{N}$, where $O_{m\times n}$ is the $m\times n$-dimensional zero matrix.

2. Its autocovariance function ${\mathrm{\Gamma }}_{ϵ}(t,s)$ defined by Definition 2.3 satisfies that ${\mathrm{\Gamma }}_{ϵ}(t,s)=\{\begin{array}{ll}{O}_{mn},& t\ne s,\\ {\mathrm{\Sigma }}_{mn},& t=s,\end{array}\mathrm{\forall }t,s\in \mathbb{N},$where $O_{mn}$ is the $mn\times mn$-dimensional zero matrix, and (7) ${\mathrm{\Sigma }}_{mn}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\sigma }_{11}^2,{\sigma }_{21}^2,\dots ,{\sigma }_{m1}^2,{\sigma }_{12}^2,{\sigma }_{22}^2,\dots ,{\sigma }_{(m-1)n}^2,{\sigma }_{mn}^2)$(7) is the $mn\times mn$-dimensional diagonal matrix with diagonal entries ${\sigma }_{11}^2$, ${\sigma }_{21}^2$, …, ${\sigma }_{m1}^2$, ${\sigma }_{12}^2$, ${\sigma }_{22}^2$, …, ${\sigma }_{(m-1)n}^2$, ${\sigma }_{mn}^2$.

For any matrix white noise $\left\{ϵ\right(t),t\in \mathbb{N}\}$, if its vectorization by columns, $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(ϵ\left(t\right)),t\in \mathbb{N}\}$, is Gaussian, then $\left\{ϵ\right(t),t\in \mathbb{N}\}$ is called a matrix Gaussian white noise.

Property 2.1

For any $m\times n$-dimensional matrix time series $\left\{ϵ\right(t),t\in \mathbb{N}\}$, it is an $m\times n$-dimensional matrix white noise if and only if $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(ϵ\left(t\right)),t\in \mathbb{N}\}$ is an mn-dimensional vector white noise, where $\mathbb{N}=\{1,2,3,\dots \}$.

The proof of Property 2.1 is not difficult, so we omit it.

When we investigate the autoregressive moving average model for matrix time series, we may use the Kronecker product, matrix reshape and derivative of matrix. Thus, we introduce them in the following.

Definition 2.5

Graham, 2018

Assume matrices $A=(a_{ij}{)}_{m\times n}$ and $C=(c_{ij}{)}_{p\times q}$, and then the $m\times n$ block matrix $(a_{ij}C{)}_{m\times n}$ is called the Kronecker product of A and C, denoted by $A\otimes C$, that is, $A\otimes C=\left[\begin{array}{cccc}{a}_{11}C& a_{12}C& \cdots & a_{1n}C\\ a_{21}C& a_{22}C& \cdots & a_{2n}C\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}C& a_{m2}C& \cdots & a_{mn}C\end{array}\right].$

Definition 2.6

For any $A=(a_{ij}{)}_{m\times n}$ and positive integers p, q satisfying pq = mn, the $(p,q)$-order reshaped matrix of A, denoted by $\mathrm{R}\mathrm{e}\mathrm{s}(A,p,q)$, is defined by $\mathrm{R}\mathrm{e}\mathrm{s}(A,p,q)=\left[\begin{array}{cccc}{a}_1& a_{p+1}& \cdots & a_{(p-1)q+1}\\ a_2& a_{p+2}& \cdots & a_{(p-1)q+2}\\ ⋮& ⋮& \ddots & ⋮\\ a_p& a_{2p}& \cdots & a_{pq}\end{array}\right],$where $a_k=a_{ij}$ for all $k=1,2,\dots ,pq$, $i=k-m\left[\right(k-1\left)/m\right]$ and $j=\left[\right(k-1\left)/m\right]+1$, where $[\cdot ]$ is the operator of taking the integer part.

Definition 2.7

Graham, 2018

Let $F=(F_{ij}{)}_{m\times n}$ and $X=(X_{ij}{)}_{p\times q}$ be two matrices, where m, n, p and q are natural numbers. The derivative of matrix F with respect to matrix X is defined by $\frac{\mathrm{\partial }F}{\mathrm{\partial }X}=\left[\begin{array}{cccc}{\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{11}}}& {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{12}}}& \cdots & {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{1q}}}\\ {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{21}}}& {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{22}}}& \cdots & {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{2q}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{p1}}}& {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{p2}}}& \cdots & {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{pq}}}\end{array}\right],$where the derivative of matrix F with respect to scalar $X_{ij}$ is defined by $\frac{\mathrm{\partial }F}{\mathrm{\partial }{X}_{ij}}=\left[\begin{array}{cccc}{\displaystyle \frac{\mathrm{\partial }{F}_{11}}{\mathrm{\partial }{X}_{ij}}}& {\displaystyle \frac{\mathrm{\partial }{F}_{12}}{\mathrm{\partial }{X}_{ij}}}& \cdots & {\displaystyle \frac{\mathrm{\partial }{F}_{1n}}{\mathrm{\partial }{X}_{ij}}}\\ {\displaystyle \frac{\mathrm{\partial }{F}_{21}}{\mathrm{\partial }{X}_{ij}}}& {\displaystyle \frac{\mathrm{\partial }{F}_{22}}{\mathrm{\partial }{X}_{ij}}}& \cdots & {\displaystyle \frac{\mathrm{\partial }{F}_{2n}}{\mathrm{\partial }{X}_{ij}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{\mathrm{\partial }{F}_{m1}}{\mathrm{\partial }{X}_{ij}}}& {\displaystyle \frac{\mathrm{\partial }{F}_{m2}}{\mathrm{\partial }{X}_{ij}}}& \cdots & {\displaystyle \frac{\mathrm{\partial }{F}_{mn}}{\mathrm{\partial }{X}_{ij}}}\end{array}\right],\begin{array}{l}i=1,2,\dots ,p,\\ j=1,2,\dots ,q.\end{array}$

For the derivative of matrix with respect to matrix, its product rule and two common formulas follow as Properties 2.2 and 2.3.

Property 2.2

Graham, 2018

For any $X=(x_{ij}{)}_{m\times n}$, $Y=(y_{ij}{)}_{n\times u}$ and $Z=(z_{ij}{)}_{p\times q}$, it follows that $\frac{\mathrm{\partial }\left(XY\right)}{\mathrm{\partial }Z}=\frac{\mathrm{\partial }X}{\mathrm{\partial }Z}(I_q\otimes Y)+(I_p\otimes X)\frac{\mathrm{\partial }Y}{\mathrm{\partial }Z},$where $I_q$ is the $q\times q$-dimensional identity matrix.

Taking $Y=(y_{ij}{)}_{n\times 1}$ and $X=Y^{\mathrm{\top }}$ into Property 2.2, we obtain Corollary 2.1.

Corollary 2.1

For any $Y=(y_{ij}{)}_{n\times 1}$ and $Z=(z_{ij}{)}_{p\times q}$, it follows that $\frac{\mathrm{\partial }\left(Y^{\mathrm{\top }}Y\right)}{\mathrm{\partial }Z}=2\frac{\mathrm{\partial }{Y}^{\mathrm{\top }}}{\mathrm{\partial }Z}(I_q\otimes Y).$

Property 2.3

Graham, 2018

For any $A=(a_{ij}{)}_{m\times n}$, $B=(b_{ij}{)}_{n\times u}$ and invertible $X=(X_{ij}{)}_{n\times n}$, it follows that $\frac{\mathrm{\partial }\mathrm{v}\mathrm{e}\mathrm{c}\left(A{X}^{-1}B\right)}{\mathrm{\partial }\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right)}=-\left(X^{-1}B\right)\otimes (A{X}^{-1}{)}^{\mathrm{\top }}$and $\frac{\mathrm{\partial }\ln \left(X\right)}{\mathrm{\partial }X}=(X^{-1}{)}^{\mathrm{\top }}.$

Taking $B=(b_{ij}{)}_{n\times 1}$ and $A=B^{\mathrm{\top }}$ into Property 2.3, we obtain Corollary 2.2.

Corollary 2.2

For any $B=(b_{ij}{)}_{n\times 1}$ and invertible $X=(X_{ij}{)}_{n\times n}$, it follows that $\frac{\mathrm{\partial }{B}^{\mathrm{\top }}{X}^{-1}B}{\mathrm{\partial }X}=-\mathrm{R}\mathrm{e}\mathrm{s}\left(\right(X^{-1}B)\otimes (A{X}^{-1}{)}^{\mathrm{\top }},n,n).$

3. Autoregressive moving average model for matrix time series

The autoregressive moving average model for matrix time series is an extension of the vector autoregressive moving average model (VARMA) to matrix time series. However, we cannot build the autoregressive moving average model for matrix time series like the VARMA model as follows: (8) $X\left(t\right)={\mathrm{\Phi }}_0+\sum _{i=1}^p{\mathrm{\Phi }}_{i}X(t-i)+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Psi }}_jϵ(t-j).$(8) The reason is that the form of (8(8) $X\left(t\right)={\mathrm{\Phi }}_0+\sum _{i=1}^p{\mathrm{\Phi }}_{i}X(t-i)+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Psi }}_jϵ(t-j).$(8) ) cannot describe the dependent relation between the different columns of $X\left(t\right)$ according to the rule of matrix multiplication. That is, the ℓth column of $X\left(t\right)$ will not be affected by the sth column of $X(t-1),X(t-2),\dots ,X(t-p)$ as $s\ne \ell$.

3.1. $\mathsf{M}\mathsf{A}\mathsf{R}\mathsf{M}\mathsf{A}(p,q)$ model

In this section, an autoregressive moving average model for matrix time series (MARMA) is first brought forward, whose degradation model, autoregressive model for matrix-valued time series (MAR), is just the model (2(2) $X_t=\sum _{k=1}^p{A}_k{X}_{t-k}{B}_k^{\mathrm{\top }}+{ϵ}_t,t=1,2,\dots ,$(2) ) proposed by Wu and Hua (2022) and the extension of model (1(1) $X_t=A{X}_{t-1}{B}^{\mathrm{\top }}+{ϵ}_t,t=1,2,\dots ,$(1) ) proposed by Chen et al. (2021).

Definition 3.1

Let $\left\{X\right(t),t\in \mathbb{N}\}$ be an $m\times n$-dimensional matrix time series. If X is stationary and for each $t\in \mathbb{N}$ it follows that (9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) where C is an $m\times n$-dimensional matrix; ${\mathrm{\Phi }}_k$ and ${\mathrm{\Theta }}_j$ are $m\times m$-dimensional matrices, and ${\mathrm{\Psi }}_k$ and ${\mathrm{\Xi }}_j$ are $n\times n$-dimensional matrices for each $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$, where p and q are two nonnegative integers; $\left\{ϵ\right(t),t\in \mathbb{N}\}$ is an $m\times n$-dimensional matrix white noise satisfying that $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)$ is independent with $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(s\left)\right)$ for all s<t, and then $\left\{X\right(t),t\in \mathbb{N}\}$ is said to follow a $(p,q)$-order autoregressive moving average model for matrix time series, denoted by MARMA$(p,q)$.

When q = 0, MARMA$(p,0)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) degenerates into the form (10) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right),$(10) which is a p-order autoregressive model for matrix time series, MAR(p).

When p = 0, MARMA$(0,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) degenerates into the form (11) $X\left(t\right)=C+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(11) which is called a q-order moving average model for matrix time series, denoted by MMA(q).

If $X=\left\{X\right(t),t\in \mathbb{N}\}$ is an $m\times n$-dimensional matrix time series defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ) and X is stationary, denote $\mu =E\left[X\right(t\left)\right],\mathrm{\forall }t\in \mathbb{N},$and then it follows from MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) that (12) $\mu =C+{\mathrm{\Phi }}_1\mu {\mathrm{\Psi }}_1+{\mathrm{\Phi }}_2\mu {\mathrm{\Psi }}_2+\cdots +{\mathrm{\Phi }}_p\mu {\mathrm{\Psi }}_p.$(12) Denote $Y\left(t\right)=X\left(t\right)-\mu .$It yields from (12(12) $\mu =C+{\mathrm{\Phi }}_1\mu {\mathrm{\Psi }}_1+{\mathrm{\Phi }}_2\mu {\mathrm{\Psi }}_2+\cdots +{\mathrm{\Phi }}_p\mu {\mathrm{\Psi }}_p.$(12) ) and MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) that (13) $Y\left(t\right)=\sum _{k=1}^p{\mathrm{\Phi }}_{k}Y(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j$(13) holds for all $t\in \mathbb{N}$, and then $Y=\left\{Y\right(t),t\in \mathbb{N}\}$ is said to follow a $(p,q)$-order centralized MARMA$(p,q)$ model.

Because every MARMA$(p,q)$ model can be changed into a centralized MARMA$(p,q)$ model and they have the same coefficient parameters. Thus, while estimating coefficient parameters of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) we will mainly study centralized MARMA$(p,q)$ model (13(13) $Y\left(t\right)=\sum _{k=1}^p{\mathrm{\Phi }}_{k}Y(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j$(13) ).

For any MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), and for any $c_k\ne 0$ and $d_j\ne 0$, $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$, it follows that $X\left(t\right)=C+\sum _{k=1}^p\left(c_k{\mathrm{\Phi }}_k\right)X(t-k)\left(\frac{1}{c_k}{\mathrm{\Psi }}_k\right)+ϵ\left(t\right)-\sum _{j=1}^q\left(d_j{\mathrm{\Theta }}_j\right)ϵ(t-j)\left(\frac{1}{d_j}{\mathrm{\Xi }}_j\right),$that is, coefficient parameters of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) are not unique! Thus, we present constraint conditions that (14) ${\mathrm{\Psi }}_k={\left({\psi }_{uv}\right)}_{n\times n}\text{ satisfies }\arg \underset{{\psi }_{ij}}{max}\{|{\psi }_{ij}|,i,j=1,2,\dots ,n\}=1$(14) and (15) ${\mathrm{\Xi }}_j={\left({\xi }_{uv}\right)}_{n\times n}\text{ satisfies }\arg \underset{{\xi }_{ij}}{max}\{|{\xi }_{ij}|,i,j=1,2,\dots ,n\}=1$(15) for all $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$.

3.2. Relationship between MARMA model and VARMA model

When the column number of matrix $X\left(t\right)$ equals one, i.e., n = 1, MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) degenerates into a $(p,q)$-order vector autoregressive moving average model, VARMA$(p,q)$, as follows: (16) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k)+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j),$(16) where $\left\{X\right(t),t\in \mathbb{N}\}$ is an m-dimensional vector time series, C is an m-dimensional vector, ${\mathrm{\Phi }}_k$ and ${\mathrm{\Theta }}_j$ are $m\times m$-dimensional matrices for all $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$, and $\left\{ϵ\right(t),t\in \mathbb{N}\}$ is a white noise of the m-dimensional vector time series satisfying that $ϵ\left(t\right)$ is independent with $X\left(s\right)$ for all s<t. Obviously, VARMA model (16(16) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k)+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j),$(16) ) is a special case of MARMA model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

On the other hand, for any $m\times n$-dimensional matrix time series $\left\{X\right(t),t\in \mathbb{N}\}$, its vectorization $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$ is an $mn\times 1$-dimensional time series, and the $(p,q)$-order vector autoregressive moving average model VARMA$(p,q)$ for $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$ follows as (17) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=A_0+\sum _{k=1}^p{A}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+ϵ\left(t\right)-\sum _{j=1}^q{B}_jϵ(t-j),$(17) where $A_0$ is an $mn\times 1$-dimensional vector; $A_k$ and $B_j$ are $mn\times mn$-dimensional matrices for $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$; and $\left\{ϵ\right(t),t\in \mathbb{N}\}$ is an $mn\times 1$-dimensional white noise satisfying that $ϵ\left(t\right)$ is independent with $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(s\left)\right)$ for all s<t.

A natural question is why the authors still bring forward MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) for $\left\{X\right(t),t\in \mathbb{N}\}$ but directly use VARMA$(p,q)$ model (17(17) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=A_0+\sum _{k=1}^p{A}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+ϵ\left(t\right)-\sum _{j=1}^q{B}_jϵ(t-j),$(17) ) for $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$.

In fact, there are two important reasons that the authors propose MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) for $\left\{X\right(t),t\in \mathbb{N}\}$. Firstly, MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) for $\left\{X\right(t),t\in \mathbb{N}\}$ can reveal the information structure of $\left\{X\right(t),t\in \mathbb{N}\}$ very clearly. Secondly, MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) for $\left\{X\right(t),t\in \mathbb{N}\}$ can reduce model parameters more greatly than VARMA$(p,q)$ model (17(17) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=A_0+\sum _{k=1}^p{A}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+ϵ\left(t\right)-\sum _{j=1}^q{B}_jϵ(t-j),$(17) ) for $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$. In fact, the parameter number of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) for $\left\{X\right(t),t\in \mathbb{N}\}$ is $2mn+(p+q)(m^2+n^2)$. However, the parameter number of VARMA$(p,q)$ model (17(17) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=A_0+\sum _{k=1}^p{A}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+ϵ\left(t\right)-\sum _{j=1}^q{B}_jϵ(t-j),$(17) ) for $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$ is $2mn+(p+q)m^2{n}^2$. Generally, $2mn+(p+q)(m^2+n^2)\ll 2mn+(p+q)m^2{n}^2.$For example, if p = q = 1 and m = n = 10, then $2mn+(p+q)(m^2+n^2)=600\ll 2mn+(p+q)m^2{n}^2=20200.$In today's big data era, m and n are often very large, taking m = n = 100 and p = q = 1 as an example, and then $2mn+(p+q)(m^2+n^2)=60000\ll 2mn+(p+q)m^2{n}^2=200020000.$

Remark 3.1

MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) greatly reduces model parameters compared with VARMA$(p,q)$ model (17(17) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=A_0+\sum _{k=1}^p{A}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+ϵ\left(t\right)-\sum _{j=1}^q{B}_jϵ(t-j),$(17) ).

Although it is not a good idea to replace MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) with VARMA$(p,q)$ model (17(17) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=A_0+\sum _{k=1}^p{A}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+ϵ\left(t\right)-\sum _{j=1}^q{B}_jϵ(t-j),$(17) ), in the following we will show there exists a special VARMA$(p,q)$ model equivalent to MARMA$(p,q)$ model, which will play a very important role in theoretical analysis of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

Theorem 3.1

MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) for $\left\{X\right(t),t\in \mathbb{N}\}$ is equivalent to VARMA$(p,q)$ model (18(18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) ) for $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X\left(t\right)),t\in \mathbb{N}\}$ as follows: (18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) where $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)$ and $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)$ represent the vectorization of matrices $X\left(t\right)$ and $ϵ\left(t\right)$ by columns, and ⊗ is the Kronecker product.

Theorem 3.1 can be proved by the following equivalence relation: for any matrices $Y_{m\times n}$, $A_{m\times m}$, $B_{m\times n}$ and $C_{n\times n}$, it follows that $Y_{m\times n}=A_{m\times m}{B}_{m\times n}{C}_{n\times n}⟺\mathrm{v}\mathrm{e}\mathrm{c}\left(Y_{m\times n}\right)=(C_{n\times n}^{\mathrm{\top }}\otimes A_{m\times m})\mathrm{v}\mathrm{e}\mathrm{c}\left(B_{m\times n}\right).$The equivalence relation is not difficult to prove, so we omit the proof and that of Theorem 3.1.

3.3. Stationary and invertible conditions for MARMA model

According to Theorem 3.1, any MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) can be converted into its corresponding VARMA$(p,q)$ model (18(18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) ). Furthermore, VARMA$(p,q)$ model (18(18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) ) can be rewritten as (19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) where (20) $\begin{array}{rl}P\left(B\right)& =I_{mn}-\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k{B}^k,\end{array}$(20) (21) $\begin{array}{rl}Q\left(B\right)& =I_{mn}-\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j{B}^j\end{array}$(21) and B is the delay operator, i.e., $BX\left(t\right)=X(t-1)$ holds for all $t\in \mathbb{N}$.

Theorem 3.2

For MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), the necessary and sufficient conditions for stationarity are that any root λ of (22(22) $|{\lambda }^p{I}_{mn}-{\lambda }^{p-1}{\mathrm{\Psi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_1-{\lambda }^{p-2}{\mathrm{\Psi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_2-\cdots -\lambda {\mathrm{\Psi }}_{p-1}^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{p-1}-{\mathrm{\Psi }}_p^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_p|=0.$(22) ) is in the unit circle, where (22) $|{\lambda }^p{I}_{mn}-{\lambda }^{p-1}{\mathrm{\Psi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_1-{\lambda }^{p-2}{\mathrm{\Psi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_2-\cdots -\lambda {\mathrm{\Psi }}_{p-1}^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{p-1}-{\mathrm{\Psi }}_p^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_p|=0.$(22) The necessary and sufficient conditions for invertibility are that any root λ of (23(23) $|{\lambda }^q{I}_{mn}-{\lambda }^{q-1}{\mathrm{\Xi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_1-{\lambda }^{q-2}{\mathrm{\Xi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_2-\cdots -\lambda {\mathrm{\Xi }}_{q-1}^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_{q-1}-{\mathrm{\Xi }}_q^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_q|=0.$(23) ) is in the unit circle, where (23) $|{\lambda }^q{I}_{mn}-{\lambda }^{q-1}{\mathrm{\Xi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_1-{\lambda }^{q-2}{\mathrm{\Xi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_2-\cdots -\lambda {\mathrm{\Xi }}_{q-1}^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_{q-1}-{\mathrm{\Xi }}_q^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_q|=0.$(23)

The proof of Theorem 3.2 is presented in Appendix 1.

Corollary 3.1

For MAR(p) model (10(10) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right),$(10) ), the necessary and sufficient conditions for stationarity are that any root λ of (22(22) $|{\lambda }^p{I}_{mn}-{\lambda }^{p-1}{\mathrm{\Psi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_1-{\lambda }^{p-2}{\mathrm{\Psi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_2-\cdots -\lambda {\mathrm{\Psi }}_{p-1}^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{p-1}-{\mathrm{\Psi }}_p^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_p|=0.$(22) ) is in the unit circle.

Remark 3.2

Corollary 3.1 expands Proposition 1 in Chen et al. (2021).

Corollary 3.2

For MMA(q) model (11(11) $X\left(t\right)=C+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(11) ), the necessary and sufficient conditions for invertibility are that any root λ of (23(23) $|{\lambda }^q{I}_{mn}-{\lambda }^{q-1}{\mathrm{\Xi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_1-{\lambda }^{q-2}{\mathrm{\Xi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_2-\cdots -\lambda {\mathrm{\Xi }}_{q-1}^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_{q-1}-{\mathrm{\Xi }}_q^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_q|=0.$(23) ) is in the unit circle.

3.4. Parameter estimation for MARMA model

In the section, we will present the conditional least square method and the conditional maximum likelihood estimation method for MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

Let $x_1,x_2,\dots ,x_N$ be a series of samples of the centralized matrix time series $X=\left\{X\right(t),t\in \mathbb{N}\}$ defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ) with $C=O_{m\times n}$, where (24) $x_t=\left[\begin{array}{cccc}{x}_{11}\left(t\right)& x_{12}\left(t\right)& \cdots & x_{1n}\left(t\right)\\ x_{21}\left(t\right)& x_{22}\left(t\right)& \cdots & x_{2n}\left(t\right)\\ ⋮& ⋮& \ddots & ⋮\\ x_{m1}\left(t\right)& x_{m2}\left(t\right)& \cdots & x_{mn}\left(t\right)\end{array}\right],t=1,2,\dots ,N,$(24) where the integer N is the sample length.

When the coefficient parameters of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) have been obtained, it follows from (12(12) $\mu =C+{\mathrm{\Phi }}_1\mu {\mathrm{\Psi }}_1+{\mathrm{\Phi }}_2\mu {\mathrm{\Psi }}_2+\cdots +{\mathrm{\Phi }}_p\mu {\mathrm{\Psi }}_p.$(12) ) that $C=\mu -{\mathrm{\Phi }}_1\mu {\mathrm{\Psi }}_1-{\mathrm{\Phi }}_2\mu {\mathrm{\Psi }}_2-\cdots -{\mathrm{\Phi }}_p\mu {\mathrm{\Psi }}_p,$and then the constant matrix C of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) can be estimated as follows: $\hat{C}=\overline{X}-{\mathrm{\Phi }}_1\overline{X}{\mathrm{\Psi }}_1-{\mathrm{\Phi }}_2\overline{X}{\mathrm{\Psi }}_2-\cdots -{\mathrm{\Phi }}_p\overline{X}{\mathrm{\Psi }}_p,$where $\overline{X}=\frac{1}{N}\sum _{t=1}^N{x}_t.$Thus, in the following we always assume the samples come from a centralized MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), i.e., $C=O_{m\times n}$.

We use VARMA$(p,q)$ model (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ) with $C=O_{m\times n}$, equivalent to centralized MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), to estimate the coefficient parameters of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) by the conditional least square method.

It yields from (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ) with $C=O_{m\times n}$ that (25) $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)=Q^{-1}\left(B\right)P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right),t\in \mathbb{N},$(25) where $Q^{-1}\left(B\right)$ is the inverse operator of $Q\left(B\right)$, and $P\left(B\right)=I_{mn}-\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k{B}^k,Q\left(B\right)=I_{mn}-\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j{B}^j.$For the sake of briefness, denote (26) $G\left(B\right)\stackrel{△}{=}\sum _{k=0}^{+\mathrm{\infty }}{G}_k{B}^k=Q^{-1}\left(B\right)P\left(B\right)$(26) and $P\left(B\right)=\sum _{i=0}^{+\mathrm{\infty }}{P}_i{B}^i,Q\left(B\right)=\sum _{j=0}^{+\mathrm{\infty }}{Q}_j{B}^j,$where we stipulate that (27) $P_i=\{\begin{array}{ll}{I}_{mn},& i=0,\\ -{\mathrm{\Psi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_i,& 1\le i\le p,\\ O_{mn},& i\ge p+1,\end{array}\text{ and }{Q}_j=\{\begin{array}{ll}{I}_{mn},& j=0,\\ -{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j{B}^j,& 1\le j\le q,\\ O_{mn},& j\ge q+1.\end{array}$(27) It follows from (26(26) $G\left(B\right)\stackrel{△}{=}\sum _{k=0}^{+\mathrm{\infty }}{G}_k{B}^k=Q^{-1}\left(B\right)P\left(B\right)$(26) ) that $Q\left(B\right)G\left(B\right)=P\left(B\right)$, which means that (28) $\sum _{i=0}^k{Q}_i{G}_{k-i}=P_k,k=0,1,2,\dots .$(28) It yields from (28(28) $\sum _{i=0}^k{Q}_i{G}_{k-i}=P_k,k=0,1,2,\dots .$(28) ) and (27(27) $P_i=\{\begin{array}{ll}{I}_{mn},& i=0,\\ -{\mathrm{\Psi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_i,& 1\le i\le p,\\ O_{mn},& i\ge p+1,\end{array}\text{ and }{Q}_j=\{\begin{array}{ll}{I}_{mn},& j=0,\\ -{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j{B}^j,& 1\le j\le q,\\ O_{mn},& j\ge q+1.\end{array}$(27) ) that (29) $G_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k+{\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& 1\le k\le p,\\ {\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& k\ge p+1,\end{array}$(29) where $k\wedge q=min\{k,q\}$.

In summary, centralized MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), i.e., $C=O_{m\times n}$, is equivalent to VARMA$(p,q)$ model (30(30) $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)=\sum _{k=0}^{+\mathrm{\infty }}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right),t\in \mathbb{N},$(30) ). (30) $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)=\sum _{k=0}^{+\mathrm{\infty }}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right),t\in \mathbb{N},$(30) where $G_k$, $k=0,1,2,\dots$, are given by (29(29) $G_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k+{\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& 1\le k\le p,\\ {\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& k\ge p+1,\end{array}$(29) ).

Theorem 3.3

According to the conditional least square method, the parameters of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) satisfy the following matrix differential equations: $\{\begin{array}{rl}& \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Phi }}_i}(I_m\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_m,i=1,2,\dots ,p,\\ & \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Psi }}_i}(I_n\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_n,i=1,2,\dots ,p,\\ & \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Theta }}_j}(I_m\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_m,j=1,2,\dots ,q,\\ & \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Xi }}_j}(I_n\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_n,j=1,2,\dots ,q,\end{array}$where $G_k$ is given by (29(29) $G_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k+{\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& 1\le k\le p,\\ {\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& k\ge p+1,\end{array}$(29) ).

The proof of Theorem 3.3 is presented in Appendix 2.

Corollary 3.3

According to the conditional least square method, the parameters of MAR(p) model (10(10) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right),$(10) ) satisfy the following matrix differential equations: $\{\begin{array}{c}{\displaystyle \sum _{t=p+1}^N\frac{\mathrm{\partial }\left({\tilde{x}}_{t-i}^{\mathrm{\top }}\right({\mathrm{\Psi }}_i\otimes {\mathrm{\Phi }}_i^{\mathrm{\top }}\left)\right)}{\mathrm{\partial }{\mathrm{\Phi }}_i}(I_m\otimes ({\tilde{x}}_t-\sum _{\ell =1}^p({\mathrm{\Psi }}_{\ell }^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{\ell }\left){\tilde{x}}_{t-\ell }\right))=O_m,}\\ {\displaystyle \sum _{t=p+1}^N\frac{\mathrm{\partial }\left({\tilde{x}}_{t-i}^{\mathrm{\top }}\right({\mathrm{\Psi }}_i\otimes {\mathrm{\Phi }}_i^{\mathrm{\top }}\left)\right)}{\mathrm{\partial }{\mathrm{\Psi }}_i}(I_n\otimes ({\tilde{x}}_t-\sum _{\ell =1}^p({\mathrm{\Psi }}_{\ell }^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{\ell }\left){\tilde{x}}_{t-\ell }\right))=O_n,}\end{array}i=1,2,\dots ,p.$

Theorem 3.4

Assume the innovations are Gaussian with the mean $O_{m\times n}$ and covariance matrix ${\mathrm{\Sigma }}_{mn}$. According to the conditional maximum likelihood estimation method, the parameters of centralized MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) satisfy the following matrix differential equations: $\{\begin{array}{rl}& \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Phi }}_i}(I_m\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_m,i=1,2,\dots ,p,\\ & \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Psi }}_i}(I_n\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_n,i=1,2,\dots ,p,\\ & \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Theta }}_j}(I_m\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_m,j=1,2,\dots ,q,\\ & \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Xi }}_j}(I_n\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_n,j=1,2,\dots ,q,\\ & \frac{1}{N}\sum _{t=1}^N\mathrm{R}\mathrm{e}\mathrm{s}(\left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right]\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right],mn,mn)={\mathrm{\Sigma }}_{mn}^{-1},\end{array}$where $H(t,G_{\cdot })={\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k}$ and $G_k$ is given by (29(29) $G_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k+{\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& 1\le k\le p,\\ {\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& k\ge p+1,\end{array}$(29) ).

The proof of Theorem 3.4 is presented in Appendix 3.

Corollary 3.4

Assume the innovations are Gaussian with the mean $O_{m\times n}$ and covariance matrix ${\mathrm{\Sigma }}_{mn}$. According to the conditional maximum likelihood estimation method, the parameters of centralized MAR(p) model (10(10) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right),$(10) ) satisfy the following matrix differential equations: $\{\begin{array}{rl}& \sum _{t=i+1}^N\frac{\mathrm{\partial }\left({\tilde{x}}_{t-i}^{\mathrm{\top }}\right({\mathrm{\Psi }}_i\otimes {\mathrm{\Phi }}_i^{\mathrm{\top }}\left)\right)}{\mathrm{\partial }{\mathrm{\Phi }}_i}(I_m\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,{\mathrm{\Phi }}_{\cdot },{\mathrm{\Psi }}_{\cdot }\left)\right])=O_m,i=1,2,\dots ,p,\\ & \sum _{t=i+1}^N\frac{\mathrm{\partial }\left({\tilde{x}}_{t-i}^{\mathrm{\top }}\right({\mathrm{\Psi }}_i\otimes {\mathrm{\Phi }}_i^{\mathrm{\top }}\left)\right)}{\mathrm{\partial }{\mathrm{\Psi }}_i}(I_n\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,{\mathrm{\Phi }}_{\cdot },{\mathrm{\Psi }}_{\cdot }\left)\right])=O_n,i=1,2,\dots ,p,\\ & \frac{1}{N}\sum _{t=1}^N\mathrm{R}\mathrm{e}\mathrm{s}(\left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,{\mathrm{\Phi }}_{\cdot },{\mathrm{\Psi }}_{\cdot }\left)\right]\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,{\mathrm{\Phi }}_{\cdot },{\mathrm{\Psi }}_{\cdot }\left)\right],mn,mn)={\mathrm{\Sigma }}_{mn}^{-1},\end{array}$where $H(t,{\mathrm{\Phi }}_{\cdot },{\mathrm{\Psi }}_{\cdot })={\tilde{x}}_t-\sum _{s=1}^{(t-1)\wedge p}({\mathrm{\Psi }}_s^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_s){\tilde{x}}_{t-s}$.

Remark 3.3

The matrix differential equations in Theorems 3.3 and 3.4 are very complex. Especially, the coefficients $G_k$ in (29(29) $G_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k+{\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& 1\le k\le p,\\ {\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& k\ge p+1,\end{array}$(29) ), $k=1,2,\dots$, are defined by a series of recursions, whose implied parameters are to be estimated. Thus, it is difficult to obtain its closed solution, but its approximate solutions can be obtained by the numerical computation method.

3.5. Hypothesis testing for the MARMA model

Let $x_1,x_2,\dots ,x_N$ be a series of samples of the centralized matrix time series $X=\left\{X\right(t),t\in \mathbb{N}\}$ defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ) with $C=O_{m\times n}$ and $x_t=(x_{ij}{)}_{m\times n}$ for all $t=1,2,\dots ,N$. Additionally assume $x_1,x_2,\dots ,x_N$ are Gaussian. In the section, we will test whether $x_1,x_2,\dots ,x_N$ follow MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

The null hypothesis and the alternative hypothesis follow as

• $H_0$: $X=\left\{X\right(t),t\in \mathbb{N}\}$ follows MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) );

• $H_1$: $X=\left\{X\right(t),t\in \mathbb{N}\}$ does not follow MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

When $H_0$ holds, denote (31) $\widetilde{ϵ}\left(t\right)=\sum _{k=0}^{t-1}{G}_k{\tilde{x}}_{t-k},t=1,2,\dots ,N,$(31) where $\widetilde{(\cdot )}=\mathrm{v}\mathrm{e}\mathrm{c}(\cdot )$. It follows from Corollary 5.3 (Karl & Simar, 2015) that $T^2=N(\overline{\widetilde{ϵ}}{)}^{\mathrm{\top }}{\left(S_{\widetilde{ϵ}}^2\right)}^{-1}\overline{\widetilde{ϵ}}\sim T^2(mn,N-1),$where (32) $\overline{\widetilde{ϵ}}=\frac{1}{N}\sum _{t=1}^N{\widetilde{ϵ}}_t\text{ and }{S}_{\widetilde{ϵ}}^2=\frac{1}{N-1}\sum _{t=1}^N{\widetilde{ϵ}}_t{\widetilde{ϵ}}_t^{\mathrm{\top }}.$(32) It follows from Theorem 5.9 (Karl & Simar, 2015) that $\frac{N-mn}{(N-1)mn}{T}^2(mn,N-1)\sim F(mn,N-mn),$that is, $F=\frac{N(N-mn)}{(N-1)mn}(\overline{\widetilde{ϵ}}{)}^{\mathrm{\top }}{\left(S_{\widetilde{ϵ}}^2\right)}^{-1}\overline{\widetilde{ϵ}}\sim F(mn,N-mn).$Summarize the above deduction and we obtain Theorem 3.5 for the hypothesis testing on MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

Theorem 3.5

For any given significance level $\alpha \in (0,1)$, if $F<F_{\frac{\alpha }{2}}(mn,N-mn)$ or $F>F_{1-\frac{\alpha }{2}}(mn,N-mn)$, then reject $\left\{X\right(t),t\in \mathbb{N}\}$ following MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ); otherwise, accept $\left\{X\right(t),t\in \mathbb{N}\}$ following MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), where $F=\frac{N(N-mn)}{(N-1)mn}(\overline{\widetilde{ϵ}}{)}^{\mathrm{\top }}{\left(S_{\widetilde{ϵ}}^2\right)}^{-1}\overline{\widetilde{ϵ}},$and $\overline{\widetilde{ϵ}}$, $S_{\widetilde{ϵ}}^2$, ${\widetilde{ϵ}}_t$ are given by (32(32) $\overline{\widetilde{ϵ}}=\frac{1}{N}\sum _{t=1}^N{\widetilde{ϵ}}_t\text{ and }{S}_{\widetilde{ϵ}}^2=\frac{1}{N-1}\sum _{t=1}^N{\widetilde{ϵ}}_t{\widetilde{ϵ}}_t^{\mathrm{\top }}.$(32) ) and (31(31) $\widetilde{ϵ}\left(t\right)=\sum _{k=0}^{t-1}{G}_k{\tilde{x}}_{t-k},t=1,2,\dots ,N,$(31) ).

3.6. Forecasting for the MARMA model

Let $\{X_t,t\in \mathbb{N}\}$ be an $m\times n$-dimensional matrix time series defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ) following MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), equivalently, $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X_t),t\in \mathbb{N}\}$ following VARMA$(p,q)$ model (18(18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) ), that is, (33) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X_t\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t-k}\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_t\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_{t-j}\right),\end{array}$(33) where $\{{ϵ}_t,t\in \mathbb{N}\}$ is an $m\times n$-dimensional matrix white noise.

Denote the forecasting for $X_{t+\ell }$ under the condition that $X_1,X_2,\dots ,X_t$ have been known by ${\hat{X}}_t\left(\ell \right)$, which refers to the ℓth step forecasting. It follows from (33(33) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X_t\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t-k}\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_t\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_{t-j}\right),\end{array}$(33) ) and the projection theorem in Hilbert space that (34) $\mathrm{v}\mathrm{e}\mathrm{c}\left({\hat{X}}_t\right(\ell \left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left({\tilde{X}}_t\right(\ell -k\left)\right)-\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left({\widetilde{ϵ}}_t\right(\ell -j\left)\right),$(34) where $\mathrm{v}\mathrm{e}\mathrm{c}\left({\tilde{X}}_t\right(k\left)\right)=\{\begin{array}{ll}\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t+k}\right),& k\le 0,\\ \mathrm{v}\mathrm{e}\mathrm{c}\left({\hat{X}}_t\right(k\left)\right),& k\ge 1,\end{array}\text{ and }\mathrm{v}\mathrm{e}\mathrm{c}\left({\widetilde{ϵ}}_t\right(k\left)\right)=\{\begin{array}{cl}\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_{t+k}\right),& k\le 0,\\ \mathrm{v}\mathrm{e}\mathrm{c}\left(O_{m\times n}\right),& k\ge 1.\end{array}$It yields from the equivalence relation of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) and VARMA$(p,q)$ model (18(18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) ) that (35) ${\hat{X}}_t\left(\ell \right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_k{\tilde{X}}_t(\ell -k){\mathrm{\Psi }}_k-\sum _{j=1}^q{\mathrm{\Theta }}_j{\widetilde{ϵ}}_t(\ell -j){\mathrm{\Xi }}_j,$(35) where ${\tilde{X}}_t\left(k\right)=\{\begin{array}{ll}{X}_{t+k},& k\le 0,\\ {\hat{X}}_t\left(k\right),& k\ge 1,\end{array}\text{ and }{\widetilde{ϵ}}_t\left(k\right)=\{\begin{array}{ll}{ϵ}_{t+k},& k\le 0,\\ O_{m\times n},& k\ge 1.\end{array}$In the following, we will study the interval estimation of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) and assume the innovations are Gaussian. Equivalently, $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(X_t),t\in \mathbb{N}\}$ follows VARMA$(p,q)$ model (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ), that is, $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X_t\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_t\right),$where $P\left(B\right)$ and $Q\left(B\right)$ are defined by (20(20) $\begin{array}{rl}P\left(B\right)& =I_{mn}-\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k{B}^k,\end{array}$(20) ) and (21(21) $\begin{array}{rl}Q\left(B\right)& =I_{mn}-\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j{B}^j\end{array}$(21) ), and $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right({ϵ}_t),t\in \mathbb{N}\}$ is a vector white noise.

Denote $\mathrm{\Pi }\left(B\right)\stackrel{△}{=}\sum _{k=0}^{+\mathrm{\infty }}{\mathrm{\Pi }}_k{B}^k=P^{-1}\left(B\right)Q\left(B\right),$and then (36) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X_t\right)=P^{-1}\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=0}^{+\mathrm{\infty }}{\mathrm{\Pi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_{t-k}\right),$(36) where (37) ${\mathrm{\Pi }}_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Xi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_k+{\displaystyle \sum _{i=1}^{k\wedge p}({\mathrm{\Psi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_i){\mathrm{\Pi }}_{k-i},}& 1\le k\le q,\\ {\displaystyle \sum _{i=1}^{k\wedge p}({\mathrm{\Psi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_i){\mathrm{\Pi }}_{k-i},}& k\ge q+1,\end{array}$(37) with $k\wedge p=min\{k,p\}$.

For any $\ell >0$, it follows from (36(36) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X_t\right)=P^{-1}\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=0}^{+\mathrm{\infty }}{\mathrm{\Pi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_{t-k}\right),$(36) ) and the estimation method of $\mathrm{v}\mathrm{e}\mathrm{c}\left({\hat{X}}_t\right(\ell \left)\right)$ that (38) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t+\ell }\right)-\mathrm{v}\mathrm{e}\mathrm{c}\left({\hat{X}}_t\right(\ell \left)\right)=\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left({ϵ}_{t+\ell -k}\right),$(38) and then (39) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t+\ell }\right)-\mathrm{v}\mathrm{e}\mathrm{c}\left({\hat{X}}_t\right(\ell \left)\right)\sim N(O_{mn\times 1},\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}).$(39) For any given $\alpha \in (0,1)$, it yields from (39(39) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t+\ell }\right)-\mathrm{v}\mathrm{e}\mathrm{c}\left({\hat{X}}_t\right(\ell \left)\right)\sim N(O_{mn\times 1},\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}).$(39) ) that the confidence interval of $\mathrm{v}\mathrm{e}\mathrm{c}\left(X_{t+\ell }\right)$ with confidence level $1-\alpha$ follows as $\left(\mathrm{v}\mathrm{e}\mathrm{c}\right({\hat{X}}_t\left(\ell \right))-U_{1-\frac{\alpha }{2}}\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}\right)},\mathrm{v}\mathrm{e}\mathrm{c}({\hat{X}}_t\left(\ell \right))+U_{1-\frac{\alpha }{2}}\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}\right)}),$where $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\cdot )$ refers to the vector composed by all main diagonal elements, and $\sqrt{\cdot }$ means taking the square roots of every elements. It yields from the equivalence relation of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) and VARMA$(p,q)$ model (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ) that the confidence interval of $X_{t+\ell }$ with confidence level $1-\alpha$ follows as $\begin{array}{rl}& \left({\hat{X}}_t\right(\ell )-U_{1-\frac{\alpha }{2}}\mathrm{R}\mathrm{e}\mathrm{s}(\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}\right)},m,n),{\hat{X}}_t(\ell )\\ & +U_{1-\frac{\alpha }{2}}\mathrm{R}\mathrm{e}\mathrm{s}(\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}\right)},m,n)).\end{array}$In summary, we can obtain the following results.

Theorem 3.6

Assume $\{X_t,t\in \mathbb{N}\}$ follows MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ).

(1) For any $\ell >0$, the ℓ-step point estimation follows as ${\hat{X}}_t\left(\ell \right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_k{\tilde{X}}_t(\ell -k){\mathrm{\Psi }}_k-\sum _{j=1}^q{\mathrm{\Theta }}_j{\widetilde{ϵ}}_t(\ell -j){\mathrm{\Xi }}_j,$where ${\tilde{X}}_t\left(k\right)=\{\begin{array}{ll}{X}_{t+k},& k\le 0,\\ {\hat{X}}_t\left(k\right),& k\ge 1,\end{array}\text{ and }{\widetilde{ϵ}}_t\left(k\right)=\{\begin{array}{ll}{ϵ}_{t+k},& k\le 0,\\ O_{m\times n},& k\ge 1.\end{array}$(2) For any $\ell >0$ and $\alpha \in (0,1)$, the ℓ-step interval estimation with confidence level $1-\alpha$ follows as $\begin{array}{rl}& \left({\hat{X}}_t\right(\ell )-U_{1-\frac{\alpha }{2}}\mathrm{R}\mathrm{e}\mathrm{s}(\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}\right)},m,n),{\hat{X}}_t(\ell )\\ & +U_{1-\frac{\alpha }{2}}\mathrm{R}\mathrm{e}\mathrm{s}(\sqrt{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\sum _{k=0}^{\ell -1}{\mathrm{\Pi }}_k{\mathrm{\Sigma }}_{mn}{\mathrm{\Pi }}_k^{\mathrm{\top }}\right)},m,n)),\end{array}$where $U_{1-\frac{\alpha }{2}}$ is the $1-\frac{\alpha }{2}$ level lower quantile of standard normal distribution, $\mathrm{R}\mathrm{e}\mathrm{s}(\cdot )$ the reshape function by Definition 2.6, $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\cdot )$ the vector composed by all main diagonal elements, $\sqrt{\cdot }$ takes the square roots of every elements, and ${\mathrm{\Pi }}_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ {\displaystyle -{\mathrm{\Xi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_k+\sum _{i=1}^{k\wedge p}({\mathrm{\Psi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_i){\mathrm{\Pi }}_{k-i},}& 1\le k\le q,\\ {\displaystyle \sum _{i=1}^{k\wedge p}({\mathrm{\Psi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_i){\mathrm{\Pi }}_{k-i},}& k\ge q+1.\end{array}$

3.7. Supplementary notes for the MARMA model

3.7.1. Model identification for the MARMA model

According to Theorem 3.1, MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) is equivalent to VARMA$(p,q)$ model (18(18) $\begin{array}{rl}\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)& =\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+\sum _{k=1}^p{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)\\ & -\sum _{j=1}^q{\mathrm{\Xi }}_j^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_j\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t-j\left)\right),\end{array}$(18) ). Thus, we can use the model identification method for the VARMA model to identify the order of MARMA model, such as $\begin{array}{rl}\mathrm{A}\mathrm{I}\mathrm{C}(p,q)=& \ln \left(|{\mathrm{\Sigma }}_{mn}\right(p,q\left)|\right)+\frac{2}{N}(p+q)(m^2+n^2),\\ \mathrm{B}\mathrm{I}\mathrm{C}(p,q)=& \ln \left(|{\mathrm{\Sigma }}_{mn}\right(p,q\left)|\right)+\frac{\ln \left(N\right)}{N}(p+q)(m^2+n^2),\end{array}$or alternatively, $\begin{array}{rl}\mathrm{A}\mathrm{I}\mathrm{C}(p,q)& =-\ln \left(L\right)+(p+q)(m^2+n^2),\\ \mathrm{B}\mathrm{I}\mathrm{C}(p,q)& =-2\ln \left(L\right)+\ln \left(N\right)(p+q)(m^2+n^2),\end{array}$where N is the length of observation sequence and $\ln \left(L\right)$ is the logarithm likelihood function.

3.7.2. MARIMA model

For any matrix time series $\left\{X\right(t)=(X_{ij}\left(t\right){)}_{m\times n},t\in \mathbb{N}\}$ defined by (3(3) $X=\{{\left(X_{ij}\right(t\left)\right)}_{m\times n},t\in \mathbb{N}\},$(3) ), the difference operator Δ for matrix time series follows as (40) $\begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =X\left(t\right)-X(t-1),\\ {\mathrm{\Delta }}^{k}X\left(t\right)& ={\mathrm{\Delta }}^{k-1}X\left(t\right)-{\mathrm{\Delta }}^{k-1}X(t-1),k=2,3,\dots ,\end{array}$(40) and Δ defined by (40(40) $\begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =X\left(t\right)-X(t-1),\\ {\mathrm{\Delta }}^{k}X\left(t\right)& ={\mathrm{\Delta }}^{k-1}X\left(t\right)-{\mathrm{\Delta }}^{k-1}X(t-1),k=2,3,\dots ,\end{array}$(40) ) has the same effect as the difference operator for vector time series. That is, if $\left\{X\right(t)=(X_{ij}\left(t\right){)}_{m\times n},t\in \mathbb{N}\}$ is nonstationary, then we can try to eliminate nonstationarity by Δ defined by (40(40) $\begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =X\left(t\right)-X(t-1),\\ {\mathrm{\Delta }}^{k}X\left(t\right)& ={\mathrm{\Delta }}^{k-1}X\left(t\right)-{\mathrm{\Delta }}^{k-1}X(t-1),k=2,3,\dots ,\end{array}$(40) ). If there exists a positive integer d such that $\left\{{\mathrm{\Delta }}^{d}X\right(t),t\in \mathbb{N}\}$ is stationary but $\left\{{\mathrm{\Delta }}^{d-1}X\right(t),t\in \mathbb{N}\}$ is nonstationary, and $\left\{{\mathrm{\Delta }}^{d}X\right(t),t\in \mathbb{N}\}$ follows a MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), then $\left\{X\right(t)=(X_{ij}\left(t\right){)}_{m\times n},t\in \mathbb{N}\}$ is called to follow a $(p,d,q)$-order autoregressive integrated moving average for matrix time series, and denoted by MARIMA$(p,d,q)$.

4. An application of the MARMA model

In this section, we will try to model the time series of daily closing prices and daily volumes of Haitong Securities Company Limited (Abbreviated as Haitong Securities; Stock code: 600837) and Ping An Insurance (Group) Company of China, Ltd. (Abbreviated as Ping An; Stock code: 601318). The data are downloaded from the China Stock Market & Accounting Research Database (CSMAR), and the time window is from January 2, 2018 to December 31, 2021, which includes 973 records every stock.

For the sake of clarity, we denote the time series by $\{\left[\begin{array}{cc}{P}_1\left(t\right)& P_2\left(t\right)\\ V_1\left(t\right)& V_2\left(t\right)\end{array}\right],t=1,2,3,\dots \},$where $P_1\left(t\right)$ and $V_1\left(t\right)$ are the daily closing price and daily volume of Haitong Securities, and $P_2\left(t\right)$ and $V_2\left(t\right)$ are the daily closing price and daily volume of Ping An.

4.1. Data preprocessing

We first conduct the Kwiatkowski, Phillips, Schmidt and Shin (KPSS) test, i.e., ‘kpsstest’ function in the software MATLAB R2020b, to test the stationarity of the daily closing prices and daily volumes of Haitong Securities and Ping An, and the results show that the daily closing prices and daily volumes of Haitong Securities and Ping An are nonstationary.

In the following we will consider the logarithmic rates (log rate) of daily closing prices and daily volumes of Haitong Securities and Ping An. Denote (41) $\left\{R\right(t)=\left[\begin{array}{cc}{R}_{11}\left(t\right)& R_{12}\left(t\right)\\ R_{21}\left(t\right)& R_{22}\left(t\right)\end{array}\right],t=2,3,4,\dots \},$(41) where $R_{1k}\left(t\right)=\ln \left(\frac{P_k\left(t\right)}{P_k(t-1)}\right)\text{ and }{R}_{2k}\left(t\right)=\ln \left(\frac{V_k\left(t\right)}{V_k(t-1)}\right),k=1,2.$That is, $R_{11}\left(t\right)$ is the logarithmic rate of daily closing price of Haitong Securities, $R_{21}\left(t\right)$ the logarithmic rate of daily volume of Haitong Securities, $R_{12}\left(t\right)$ the logarithmic rate of daily closing price of Ping An and $R_{22}\left(t\right)$ the logarithmic rate of daily volume of Ping An.

We conduct the Kwiatkowski, Phillips, Schmidt and Shin (KPSS) test, i.e., ‘kpsstest’ function in the software MATLAB R2020b, to test the stationarity of the logarithmic rates of daily closing prices and daily volumes of Haitong Securities and Ping An, and the results show that the logarithmic rates of daily closing prices and daily volumes of Haitong Securities and Ping An are stationary.

Additionally, we conduct a Ljung-Box Q test, i.e., ‘lbqtest’ function in the software MATLAB R2020b, to test the pure randomness of the logarithmic rates of daily closing prices and daily volumes of Haitong Securities and Ping An, and the results show that the logarithmic rates of daily closing prices or daily volumes of Haitong Securities and Ping An are not purely random.

In conclusion, for the stocks of Haitong Securities and Ping An, their daily closing prices and daily volumes are nonstationary, but their logarithmic rates of daily closing prices and daily volumes are stationary, and their logarithmic rates of daily closing prices or daily volumes are not purely random.

4.2. Modelling of $\mathsf{M}\mathsf{A}\mathsf{R}\mathsf{M}\mathsf{A}(p,q)$

We use the Bayesian information criterion (BIC) to select the model, and the results show that MARMA(4,0) is the best. Using the conditional least square method and MATLAB R2020b program, we establish MARMA(4,0) model for $\left\{R\right(t),t=1,2,3,\dots \}$ by (41(41) $\left\{R\right(t)=\left[\begin{array}{cc}{R}_{11}\left(t\right)& R_{12}\left(t\right)\\ R_{21}\left(t\right)& R_{22}\left(t\right)\end{array}\right],t=2,3,4,\dots \},$(41) ) as follows: (42) $R\left(t\right)=\hat{C}+{\hat{\mathrm{\Phi }}}_{1}R(t-1){\hat{\mathrm{\Psi }}}_1+{\hat{\mathrm{\Phi }}}_{2}R(t-2){\hat{\mathrm{\Psi }}}_2+{\hat{\mathrm{\Phi }}}_{3}R(t-3){\hat{\mathrm{\Psi }}}_3+{\hat{\mathrm{\Phi }}}_{4}R(t-4){\hat{\mathrm{\Psi }}}_4+ϵ\left(t\right),$(42) where $\begin{array}{rl}{\hat{\mathrm{\Phi }}}_1& =\left[\begin{array}{cc}2.11\times {10}^{-2}& 1.32\times {10}^{-3}\\ 3.4371& -0.3527\end{array}\right],{\hat{\mathrm{\Psi }}}_1=\left[\begin{array}{cc}1& 0.3708\\ -0.4014& 0.0354\end{array}\right],\\ {\hat{\mathrm{\Phi }}}_2& =\left[\begin{array}{cc}0.0339& -2.21\times {10}^{-3}\\ 1.1721& -0.1271\end{array}\right],{\hat{\mathrm{\Psi }}}_2=\left[\begin{array}{cc}0.4132& 0.4468\\ 1& 0.2409\end{array}\right],\\ {\hat{\mathrm{\Phi }}}_3& =\left[\begin{array}{cc}0.0101& -9.82\times {10}^{-4}\\ 0.3575& -0.1739\end{array}\right],{\hat{\mathrm{\Psi }}}_3=\left[\begin{array}{cc}1& -0.4715\\ -0.1625& 0.2516\end{array}\right],\\ {\hat{\mathrm{\Phi }}}_4& =\left[\begin{array}{cc}-0.0440& 8.75\times {10}^{-4}\\ -0.4732& -0.1012\end{array}\right],{\hat{\mathrm{\Psi }}}_4=\left[\begin{array}{cc}0.9537& 0.6853\\ 0.4385& 1\end{array}\right],\\ \hat{\mu }& =\left[\begin{array}{cc}-6.66\times {10}^{-5}& -3.75\times {10}^{-4}\\ 4.86\times {10}^{-4}& -1.64\times {10}^{-3}\end{array}\right],\end{array}$and $\hat{C}=\left[\begin{array}{cc}-3.27\times {10}^{-5}& -3.48\times {10}^{-4}\\ 1.27\times {10}^{-3}& -3.40\times {10}^{-3}\end{array}\right],$and then the covariance matrix of residuals $\left\{ϵ\right(t),t\in \mathbb{N}\}$ follows as (43) ${\mathrm{\Sigma }}_{ϵ}=\left[\begin{array}{cccc}4.40\times {10}^{-4}& 2.74\times {10}^{-3}& 2.17\times {10}^{-4}& 1.30\times {10}^{-3}\\ 2.74\times {10}^{-3}& 0.1386& 1.31\times {10}^{-3}& 5.60\times {10}^{-2}\\ 2.17\times {10}^{-4}& 1.31\times {10}^{-3}& 3.15\times {10}^{-4}& 1.31\times {10}^{-3}\\ 1.30\times {10}^{-3}& 5.60\times {10}^{-2}& 1.31\times {10}^{-3}& 0.1077\end{array}\right].$(43)

4.3. Evaluation on $\mathsf{M}\mathsf{A}\mathsf{R}\mathsf{M}\mathsf{A}(p,q)$

For the sake of saving space, we will not show the model test, model optimization or forecasting of MARMA(4,0) model (42(42) $R\left(t\right)=\hat{C}+{\hat{\mathrm{\Phi }}}_{1}R(t-1){\hat{\mathrm{\Psi }}}_1+{\hat{\mathrm{\Phi }}}_{2}R(t-2){\hat{\mathrm{\Psi }}}_2+{\hat{\mathrm{\Phi }}}_{3}R(t-3){\hat{\mathrm{\Psi }}}_3+{\hat{\mathrm{\Phi }}}_{4}R(t-4){\hat{\mathrm{\Psi }}}_4+ϵ\left(t\right),$(42) ), but present a comparison of the MARMA model and ARMA model in this subsection. We first establish ARMA$(p,q)$ model for $R_{11}\left(t\right),R_{21}\left(t\right),R_{12}\left(t\right)$ and $R_{22}\left(t\right),$ respectively, and obtain their models as follows: (44) $\begin{array}{rl}{R}_{11}\left(t\right)& =1.53\times {10}^{-8}+0.0204R_{11}(t-1)+0.0058R_{11}(t-2)\\ & +0.0625R_{11}(t-3)-0.0393R_{11}(t-4)+e_{11}\left(t\right),\\ R_{21}\left(t\right)& =-2.98\times {10}^{-4}-0.3994R_{21}(t-1)-0.2598R_{21}(t-2)\\ & -0.1918R_{21}(t-3)-0.1096R_{21}(t-4)+e_{21}\left(t\right),\\ R_{12}\left(t\right)& =5.17\times {10}^{-7}-0.0002R_{12}(t-1)-0.0041R_{12}(t-2)\\ & +0.0476R_{12}(t-3)-0.0699R_{12}(t-4)+e_{12}\left(t\right),\\ R_{22}\left(t\right)& =2.52\times {10}^{-4}-0.4842R_{22}(t-1)-0.3860R_{22}(t-2)\\ & -0.2660R_{22}(t-3)-0.1562R_{22}(t-4)+e_{22}\left(t\right),\end{array}$(44) where the covariance matrix of residuals $\left\{e\right(t)\stackrel{△}{=}(e_{11}\left(t\right),e_{21}\left(t\right),e_{12}\left(t\right),e_{22}\left(t\right){)}^{\mathrm{\top }},t\in \mathbb{N}\}$ follows as (45) ${\mathrm{\Sigma }}_e=\left[\begin{array}{cccc}4.49\times {10}^{-4}& 2.69\times {10}^{-3}& 2.22\times {10}^{-4}& 1.24\times {10}^{-3}\\ 2.69\times {10}^{-3}& 0.1692& 1.40\times {10}^{-3}& 7.22\times {10}^{-2}\\ 2.22\times {10}^{-4}& 1.40\times {10}^{-3}& 3.22\times {10}^{-4}& 1.32\times {10}^{-3}\\ 1.24\times {10}^{-3}& 7.22\times {10}^{-2}& 1.32\times {10}^{-3}& 0.1369\end{array}\right].$(45) It follows from (43(43) ${\mathrm{\Sigma }}_{ϵ}=\left[\begin{array}{cccc}4.40\times {10}^{-4}& 2.74\times {10}^{-3}& 2.17\times {10}^{-4}& 1.30\times {10}^{-3}\\ 2.74\times {10}^{-3}& 0.1386& 1.31\times {10}^{-3}& 5.60\times {10}^{-2}\\ 2.17\times {10}^{-4}& 1.31\times {10}^{-3}& 3.15\times {10}^{-4}& 1.31\times {10}^{-3}\\ 1.30\times {10}^{-3}& 5.60\times {10}^{-2}& 1.31\times {10}^{-3}& 0.1077\end{array}\right].$(43) ) and (45(45) ${\mathrm{\Sigma }}_e=\left[\begin{array}{cccc}4.49\times {10}^{-4}& 2.69\times {10}^{-3}& 2.22\times {10}^{-4}& 1.24\times {10}^{-3}\\ 2.69\times {10}^{-3}& 0.1692& 1.40\times {10}^{-3}& 7.22\times {10}^{-2}\\ 2.22\times {10}^{-4}& 1.40\times {10}^{-3}& 3.22\times {10}^{-4}& 1.32\times {10}^{-3}\\ 1.24\times {10}^{-3}& 7.22\times {10}^{-2}& 1.32\times {10}^{-3}& 0.1369\end{array}\right].$(45) ) that the residuals of MARMA(4,0) model (42(42) $R\left(t\right)=\hat{C}+{\hat{\mathrm{\Phi }}}_{1}R(t-1){\hat{\mathrm{\Psi }}}_1+{\hat{\mathrm{\Phi }}}_{2}R(t-2){\hat{\mathrm{\Psi }}}_2+{\hat{\mathrm{\Phi }}}_{3}R(t-3){\hat{\mathrm{\Psi }}}_3+{\hat{\mathrm{\Phi }}}_{4}R(t-4){\hat{\mathrm{\Psi }}}_4+ϵ\left(t\right),$(42) ) are almost consistently less than those of ARMA(4,0) model (44(44) $\begin{array}{rl}{R}_{11}\left(t\right)& =1.53\times {10}^{-8}+0.0204R_{11}(t-1)+0.0058R_{11}(t-2)\\ & +0.0625R_{11}(t-3)-0.0393R_{11}(t-4)+e_{11}\left(t\right),\\ R_{21}\left(t\right)& =-2.98\times {10}^{-4}-0.3994R_{21}(t-1)-0.2598R_{21}(t-2)\\ & -0.1918R_{21}(t-3)-0.1096R_{21}(t-4)+e_{21}\left(t\right),\\ R_{12}\left(t\right)& =5.17\times {10}^{-7}-0.0002R_{12}(t-1)-0.0041R_{12}(t-2)\\ & +0.0476R_{12}(t-3)-0.0699R_{12}(t-4)+e_{12}\left(t\right),\\ R_{22}\left(t\right)& =2.52\times {10}^{-4}-0.4842R_{22}(t-1)-0.3860R_{22}(t-2)\\ & -0.2660R_{22}(t-3)-0.1562R_{22}(t-4)+e_{22}\left(t\right),\end{array}$(44) ).

In practice, we are more concerned about the residual variance, i.e., the variance of every element of residual. Using (43(43) ${\mathrm{\Sigma }}_{ϵ}=\left[\begin{array}{cccc}4.40\times {10}^{-4}& 2.74\times {10}^{-3}& 2.17\times {10}^{-4}& 1.30\times {10}^{-3}\\ 2.74\times {10}^{-3}& 0.1386& 1.31\times {10}^{-3}& 5.60\times {10}^{-2}\\ 2.17\times {10}^{-4}& 1.31\times {10}^{-3}& 3.15\times {10}^{-4}& 1.31\times {10}^{-3}\\ 1.30\times {10}^{-3}& 5.60\times {10}^{-2}& 1.31\times {10}^{-3}& 0.1077\end{array}\right].$(43) ) and (45(45) ${\mathrm{\Sigma }}_e=\left[\begin{array}{cccc}4.49\times {10}^{-4}& 2.69\times {10}^{-3}& 2.22\times {10}^{-4}& 1.24\times {10}^{-3}\\ 2.69\times {10}^{-3}& 0.1692& 1.40\times {10}^{-3}& 7.22\times {10}^{-2}\\ 2.22\times {10}^{-4}& 1.40\times {10}^{-3}& 3.22\times {10}^{-4}& 1.32\times {10}^{-3}\\ 1.24\times {10}^{-3}& 7.22\times {10}^{-2}& 1.32\times {10}^{-3}& 0.1369\end{array}\right].$(45) ), we compute the relative change of the residual variance of MARMA(4,0) model (42(42) $R\left(t\right)=\hat{C}+{\hat{\mathrm{\Phi }}}_{1}R(t-1){\hat{\mathrm{\Psi }}}_1+{\hat{\mathrm{\Phi }}}_{2}R(t-2){\hat{\mathrm{\Psi }}}_2+{\hat{\mathrm{\Phi }}}_{3}R(t-3){\hat{\mathrm{\Psi }}}_3+{\hat{\mathrm{\Phi }}}_{4}R(t-4){\hat{\mathrm{\Psi }}}_4+ϵ\left(t\right),$(42) ) to the residual variance of ARMA(4,0) model (44(44) $\begin{array}{rl}{R}_{11}\left(t\right)& =1.53\times {10}^{-8}+0.0204R_{11}(t-1)+0.0058R_{11}(t-2)\\ & +0.0625R_{11}(t-3)-0.0393R_{11}(t-4)+e_{11}\left(t\right),\\ R_{21}\left(t\right)& =-2.98\times {10}^{-4}-0.3994R_{21}(t-1)-0.2598R_{21}(t-2)\\ & -0.1918R_{21}(t-3)-0.1096R_{21}(t-4)+e_{21}\left(t\right),\\ R_{12}\left(t\right)& =5.17\times {10}^{-7}-0.0002R_{12}(t-1)-0.0041R_{12}(t-2)\\ & +0.0476R_{12}(t-3)-0.0699R_{12}(t-4)+e_{12}\left(t\right),\\ R_{22}\left(t\right)& =2.52\times {10}^{-4}-0.4842R_{22}(t-1)-0.3860R_{22}(t-2)\\ & -0.2660R_{22}(t-3)-0.1562R_{22}(t-4)+e_{22}\left(t\right),\end{array}$(44) ) as follows: $\left[\begin{array}{cc}{\displaystyle \frac{\mathrm{v}\mathrm{a}\mathrm{r}\left({ϵ}_{11}\right(t\left)\right)-\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{11}\right(t\left)\right)}{\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{11}\right(t\left)\right)}}& {\displaystyle \frac{\mathrm{v}\mathrm{a}\mathrm{r}\left({ϵ}_{12}\right(t\left)\right)-\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{12}\right(t\left)\right)}{\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{12}\right(t\left)\right)}}\\ {\displaystyle \frac{\mathrm{v}\mathrm{a}\mathrm{r}\left({ϵ}_{21}\right(t\left)\right)-\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{21}\right(t\left)\right)}{\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{21}\right(t\left)\right)}}& {\displaystyle \frac{\mathrm{v}\mathrm{a}\mathrm{r}\left({ϵ}_{22}\right(t\left)\right)-\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{22}\right(t\left)\right)}{\mathrm{v}\mathrm{a}\mathrm{r}\left(e_{22}\right(t\left)\right)}}\end{array}\right]=\left[\begin{array}{cc}-1.93\mathrm{\%}& -2.24\mathrm{\%}\\ -18.10\mathrm{\%}& -21.31\mathrm{\%}\end{array}\right].$That is, MARMA(4,0) model (42(42) $R\left(t\right)=\hat{C}+{\hat{\mathrm{\Phi }}}_{1}R(t-1){\hat{\mathrm{\Psi }}}_1+{\hat{\mathrm{\Phi }}}_{2}R(t-2){\hat{\mathrm{\Psi }}}_2+{\hat{\mathrm{\Phi }}}_{3}R(t-3){\hat{\mathrm{\Psi }}}_3+{\hat{\mathrm{\Phi }}}_{4}R(t-4){\hat{\mathrm{\Psi }}}_4+ϵ\left(t\right),$(42) ) reduces all residual variance relative to ARMA(4,0) model (44(44) $\begin{array}{rl}{R}_{11}\left(t\right)& =1.53\times {10}^{-8}+0.0204R_{11}(t-1)+0.0058R_{11}(t-2)\\ & +0.0625R_{11}(t-3)-0.0393R_{11}(t-4)+e_{11}\left(t\right),\\ R_{21}\left(t\right)& =-2.98\times {10}^{-4}-0.3994R_{21}(t-1)-0.2598R_{21}(t-2)\\ & -0.1918R_{21}(t-3)-0.1096R_{21}(t-4)+e_{21}\left(t\right),\\ R_{12}\left(t\right)& =5.17\times {10}^{-7}-0.0002R_{12}(t-1)-0.0041R_{12}(t-2)\\ & +0.0476R_{12}(t-3)-0.0699R_{12}(t-4)+e_{12}\left(t\right),\\ R_{22}\left(t\right)& =2.52\times {10}^{-4}-0.4842R_{22}(t-1)-0.3860R_{22}(t-2)\\ & -0.2660R_{22}(t-3)-0.1562R_{22}(t-4)+e_{22}\left(t\right),\end{array}$(44) ). Especially, the relative change of volume's residual variance exceeds $10\mathrm{\%}$ by MARMA(4,0) model (42(42) $R\left(t\right)=\hat{C}+{\hat{\mathrm{\Phi }}}_{1}R(t-1){\hat{\mathrm{\Psi }}}_1+{\hat{\mathrm{\Phi }}}_{2}R(t-2){\hat{\mathrm{\Psi }}}_2+{\hat{\mathrm{\Phi }}}_{3}R(t-3){\hat{\mathrm{\Psi }}}_3+{\hat{\mathrm{\Phi }}}_{4}R(t-4){\hat{\mathrm{\Psi }}}_4+ϵ\left(t\right),$(42) ) relative to ARMA(4,0) model (44(44) $\begin{array}{rl}{R}_{11}\left(t\right)& =1.53\times {10}^{-8}+0.0204R_{11}(t-1)+0.0058R_{11}(t-2)\\ & +0.0625R_{11}(t-3)-0.0393R_{11}(t-4)+e_{11}\left(t\right),\\ R_{21}\left(t\right)& =-2.98\times {10}^{-4}-0.3994R_{21}(t-1)-0.2598R_{21}(t-2)\\ & -0.1918R_{21}(t-3)-0.1096R_{21}(t-4)+e_{21}\left(t\right),\\ R_{12}\left(t\right)& =5.17\times {10}^{-7}-0.0002R_{12}(t-1)-0.0041R_{12}(t-2)\\ & +0.0476R_{12}(t-3)-0.0699R_{12}(t-4)+e_{12}\left(t\right),\\ R_{22}\left(t\right)& =2.52\times {10}^{-4}-0.4842R_{22}(t-1)-0.3860R_{22}(t-2)\\ & -0.2660R_{22}(t-3)-0.1562R_{22}(t-4)+e_{22}\left(t\right),\end{array}$(44) ), which means the MARMA model could really improve the prediction accuracy.

5. Conclusion

We proposed an autoregressive moving average model for matrix time series (MARMA), which is an extension of the autoregressive model for matrix time series (MAR). Like the MAR model, the MARMA model retains the original matrix structure, and provides a much more parsimonious model, compared with the approach of the vector autoregressive model for vectorizing the matrix into a long vector. Compared with MAR model, MARMA models are capable of modelling the unknown process with the minimum number of parameters.

As for MARMA model, the necessary and sufficient conditions for stationarity and invertibility are established. Parameter estimation methods are investigated for the conditional least square method and the conditional maximum likelihood estimation method. Point forecasting and interval forecasting are presented by using the projection theorem in the Hilbert space and the decomposition technique of time series. Additionally, model identification, model testing and possible extensions are discussed.

There are many directions to extend the scope of the MARMA model. Random environment such as the Markov environment might be imposed on the MARMA model to depict the impact of environmental change. Additionally, sparsity or group sparsity might be imposed on coefficient matrices to reach a further dimension reduction. Furthermore, the idea of MARMA can be applied for yield modelling, volatility modelling, weather forecast modelling and animal migration modelling.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This paper is partially supported by the basic scientific research business expenses of Universities in Xinjiang, China [Grant Number XQZX20230057] and the National Natural Science Foundation of China [Grant Number 11671142].

Appendices

Appendix 1. Proof of Theorem 3.2

In order to obtain stationary conditions and invertible conditions for MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), we first give a lemma as follows.

Lemma A.1

For any square matrices $A_1,A_2,\dots ,A_k$, the operator $G\left(B\right)\stackrel{△}{=}I-A_{1}B-\cdots -A_{k-1}{B}^{k-1}-A_k{B}^k$is invertible if and only if any root λ of (A1(A1) $|{\lambda }^{k}I-{\lambda }^{k-1}{A}_1-\cdots -\lambda A_{k-1}-A_k|=0.$(A1) ) satisfies $|\lambda |<1$, where k is a natural number and B is the delay operator. (A1) $|{\lambda }^{k}I-{\lambda }^{k-1}{A}_1-\cdots -\lambda A_{k-1}-A_k|=0.$(A1)

Proof.

For the polynomial with k degree and matrix coefficients $G\left(z\right)=I-A_{1}z-\cdots -A_{k-1}{z}^{k-1}-A_k{z}^k,$it can be factorized into k linear polynomials with the matrix coefficient in the complex field as follows: $G\left(z\right)=(I-C_{1}z)(I-C_{2}z)\cdots (I-C_{k}z),$where $C_1,C_2,\dots ,C_k$ are determined by (A2) $\sum _{1\le j_1<j_2<\cdots <j_u\le k}{C}_{j_1}{C}_{j_2}\cdots C_{j_u}=(-1{)}^{u-1}{A}_u,u=1,2,\dots ,k.$(A2) Thus, (A3) $G\left(B\right)=(I-C_{1}B)(I-C_{2}B)\cdots (I-C_{k}B).$(A3) For any $i=1,2,\dots ,k$, it is easy to prove that $I-C_{i}B$ is invertible if and only if $\rho \left(C_i\right)<1$, that is, all roots of $|\lambda I-C_i|=0$ are in the unit circle. It follows from (A3(A3) $G\left(B\right)=(I-C_{1}B)(I-C_{2}B)\cdots (I-C_{k}B).$(A3) ) that $G\left(B\right)$ is invertible if and only if all $I-C_{i}B$, $i=1,2,\dots ,k$, are invertible. Thus, $G\left(B\right)$ is invertible if and only if all roots of $|\lambda I-C_i|=0$ are in the unit circle for all $i=1,2,\dots ,k$. According to determinant properties, $G\left(B\right)$ is invertible if and only if all roots of (A4) $\left|\right(\lambda I-C_1\left)\right(\lambda I-C_2)\cdots (\lambda I-C_k\left)\right|=0$(A4) are in the unit circle. It yields from (A2(A2) $\sum _{1\le j_1<j_2<\cdots <j_u\le k}{C}_{j_1}{C}_{j_2}\cdots C_{j_u}=(-1{)}^{u-1}{A}_u,u=1,2,\dots ,k.$(A2) ) that $(\lambda I-C_1)(\lambda I-C_2)\cdots (\lambda I-C_k)={\lambda }^{k}I-{\lambda }^{k-1}{A}_1-\cdots -\lambda A_{k-1}-A_k.$Thus, $G\left(B\right)$ is invertible if and only if all roots of $|{\lambda }^{k}I-{\lambda }^{k-1}{A}_1-\cdots -\lambda A_{k-1}-A_k|=0$are in the unit circle.

Proof

Proof of Theorem 3.2

For VARMA$(p,q)$ model (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ), $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N}.$It follows from the concept of stationarity that the necessary and sufficient conditions of stationarity are that the operator $P\left(B\right)$ is invertible. According to Lemma A.1, the operator $P\left(B\right)$ is invertible if and only if any root λ of (22(22) $|{\lambda }^p{I}_{mn}-{\lambda }^{p-1}{\mathrm{\Psi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_1-{\lambda }^{p-2}{\mathrm{\Psi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_2-\cdots -\lambda {\mathrm{\Psi }}_{p-1}^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{p-1}-{\mathrm{\Psi }}_p^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_p|=0.$(22) ) satisfies $|\lambda |<1$. Thus, VARMA$(p,q)$ model (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ) is stationary if and only if any root λ of (22(22) $|{\lambda }^p{I}_{mn}-{\lambda }^{p-1}{\mathrm{\Psi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_1-{\lambda }^{p-2}{\mathrm{\Psi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_2-\cdots -\lambda {\mathrm{\Psi }}_{p-1}^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{p-1}-{\mathrm{\Psi }}_p^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_p|=0.$(22) ) satisfies $|\lambda |<1$. Note that VARMA$(p,q)$ model (19(19) $P\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=\mathrm{v}\mathrm{e}\mathrm{c}\left(C\right)+Q\left(B\right)\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N},$(19) ) is equivalent to MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ), so MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) is stationary if and only if any root λ of (22(22) $|{\lambda }^p{I}_{mn}-{\lambda }^{p-1}{\mathrm{\Psi }}_1^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_1-{\lambda }^{p-2}{\mathrm{\Psi }}_2^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_2-\cdots -\lambda {\mathrm{\Psi }}_{p-1}^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_{p-1}-{\mathrm{\Psi }}_p^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_p|=0.$(22) ) satisfies $|\lambda |<1$.

The necessary and sufficient conditions for invertibility can be obtained by the similar method to obtain the necessary and sufficient conditions for stationarity, so we omit it.

Appendix 2. Proof of Theorem 3.3

Noting that $\left\{\mathrm{v}\mathrm{e}\mathrm{c}\right(ϵ\left(t\right)),t\in \mathbb{N}\}$ is an $mn\times 1$-dimensional white noise, and the objective function of VARMA$(p,q)$ model (30(30) $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)=\sum _{k=0}^{+\mathrm{\infty }}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right),t\in \mathbb{N},$(30) ) using the conditional least square method follows as (A5) $\begin{array}{rl}& J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1\dots ,{\mathrm{\Xi }}_q)\\ & =\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right),\end{array}$(A5) where we take $x_t=O_{m\times n}$ for all $t\le 0$.

Lemma A.2

$J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)$ defined by (A5(A5) $\begin{array}{rl}& J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1\dots ,{\mathrm{\Xi }}_q)\\ & =\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right),\end{array}$(A5) ) has the minimum value about ${\mathrm{\Phi }}_k$, ${\mathrm{\Psi }}_k$, ${\mathrm{\Theta }}_j$ and ${\mathrm{\Xi }}_j$ for all $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$.

Proof.

It yields from analysing (29(29) $G_k=\{\begin{array}{ll}{I}_{mn},& k=0,\\ -{\mathrm{\Psi }}_k^{\mathrm{\top }}\otimes {\mathrm{\Phi }}_k+{\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& 1\le k\le p,\\ {\displaystyle \sum _{i=1}^{k\wedge q}({\mathrm{\Xi }}_i^{\mathrm{\top }}\otimes {\mathrm{\Theta }}_i)G_{k-i},}& k\ge p+1,\end{array}$(29) ) that $J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)$ by (A5(A5) $\begin{array}{rl}& J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1\dots ,{\mathrm{\Xi }}_q)\\ & =\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right),\end{array}$(A5) ) is a multivariate polynomial of ${\mathrm{\Phi }}_k$, ${\mathrm{\Psi }}_k$, ${\mathrm{\Theta }}_j$ and ${\mathrm{\Xi }}_j$ for all $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$. And it is obvious that $J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)$ by (A5(A5) $\begin{array}{rl}& J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1\dots ,{\mathrm{\Xi }}_q)\\ & =\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right),\end{array}$(A5) ) is greater than or equal to zero, which means that $J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)$ by (A5(A5) $\begin{array}{rl}& J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1\dots ,{\mathrm{\Xi }}_q)\\ & =\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right),\end{array}$(A5) ) has lower bound. Thus, $J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)$ by (A5(A5) $\begin{array}{rl}& J({\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1\dots ,{\mathrm{\Xi }}_q)\\ & =\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right),\end{array}$(A5) ) has the minimum value about ${\mathrm{\Phi }}_k$, ${\mathrm{\Psi }}_k$, ${\mathrm{\Theta }}_j$ and ${\mathrm{\Xi }}_j$ for all $k=1,2,\dots ,p$ and $j=1,2,\dots ,q$.

Proof

Proof of Theorem 3.3.

It follows from Lemma A.2 that, according to the conditional least square method, the parameters of MARMA$(p,q)$ model (9(9) $X\left(t\right)=C+\sum _{k=1}^p{\mathrm{\Phi }}_{k}X(t-k){\mathrm{\Psi }}_k+ϵ\left(t\right)-\sum _{j=1}^q{\mathrm{\Theta }}_jϵ(t-j){\mathrm{\Xi }}_j,$(9) ) satisfy the following matrix differential equations: $\{\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Phi }}_i}\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{\ell =1}^{t-1}{G}_{\ell }\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-\ell }\left)\right)=O_m,i=1,2,\dots ,p,\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Psi }}_i}\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{\ell =1}^{t-1}{G}_{\ell }\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-\ell }\left)\right)=O_n,i=1,2,\dots ,p,\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Theta }}_j}\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{\ell =1}^{t-1}{G}_{\ell }\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-\ell }\left)\right)=O_m,j=1,2,\dots ,q,\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Xi }}_j}\sum _{t=p+1}^N{\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{k=1}^{t-1}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-k}\left)\right)}^{\mathrm{\top }}\left(\mathrm{v}\mathrm{e}\mathrm{c}\right(x_t)+\sum _{\ell =1}^{t-1}{G}_{\ell }\mathrm{v}\mathrm{e}\mathrm{c}(x_{t-\ell }\left)\right)=O_n,j=1,2,\dots ,q.\end{array}$Using the derivative of scalar by matrix, it yields from Corollary 2.1 that $\{\begin{array}{rl}& \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Phi }}_i}(I_m\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_m,i=1,2,\dots ,p,\\ & \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Psi }}_i}(I_n\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_n,i=1,2,\dots ,p,\\ & \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Theta }}_j}(I_m\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_m,j=1,2,\dots ,q,\\ & \sum _{t=p+1}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Xi }}_j}(I_n\otimes ({\tilde{x}}_t+\sum _{\ell =1}^{t-1}{G}_{\ell }{\tilde{x}}_{t-\ell }))=O_n,j=1,2,\dots ,q.\end{array}\text{■}$

Appendix 3. Proof of Theorem 3.4

It yields from (30(30) $\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right)=\sum _{k=0}^{+\mathrm{\infty }}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right),t\in \mathbb{N},$(30) ) that (A6) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=-\sum _{k=1}^{+\mathrm{\infty }}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N}.$(A6) For the sake of briefness, we denote ${\tilde{X}}_t=\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right),t\in \mathbb{N}\text{ and }{\tilde{x}}_k=\mathrm{v}\mathrm{e}\mathrm{c}\left(x_k\right),k=1,2,\dots ,N.$It yields from (A6(A6) $\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t\left)\right)=-\sum _{k=1}^{+\mathrm{\infty }}{G}_k\mathrm{v}\mathrm{e}\mathrm{c}\left(X\right(t-k\left)\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(ϵ\right(t\left)\right),t\in \mathbb{N}.$(A6) ) that (A7) ${\tilde{X}}_t|\{{\tilde{X}}_{t-1},{\tilde{X}}_{t-2},\dots \}\sim N(-\sum _{k=1}^{+\mathrm{\infty }}{G}_k{\tilde{X}}_{t-k},{\mathrm{\Sigma }}_{mn}),t\in \mathbb{N},$(A7) where ${\mathrm{\Sigma }}_{mn}$ is defined by (7(7) ${\mathrm{\Sigma }}_{mn}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\sigma }_{11}^2,{\sigma }_{21}^2,\dots ,{\sigma }_{m1}^2,{\sigma }_{12}^2,{\sigma }_{22}^2,\dots ,{\sigma }_{(m-1)n}^2,{\sigma }_{mn}^2)$(7) ).

Let $X\left(t\right)=O_{m\times n}$ for all $t\le 0$. It follows from (A7(A7) ${\tilde{X}}_t|\{{\tilde{X}}_{t-1},{\tilde{X}}_{t-2},\dots \}\sim N(-\sum _{k=1}^{+\mathrm{\infty }}{G}_k{\tilde{X}}_{t-k},{\mathrm{\Sigma }}_{mn}),t\in \mathbb{N},$(A7) ) that (A8) ${\tilde{X}}_1\sim N(O_{mn\times 1},{\mathrm{\Sigma }}_{mn})$(A8) and (A9) ${\tilde{X}}_t|\{{\tilde{X}}_{t-1},{\tilde{X}}_{t-2},\dots ,{\tilde{X}}_1\}\sim N(-\sum _{k=1}^{t-1}{G}_k{\tilde{X}}_{t-k},{\mathrm{\Sigma }}_{mn}),t\in \mathbb{N}.$(A9) Thus, the maximum likelihood function of $x_1,x_2,\dots ,x_V$ follows as $\begin{array}{rl}& L(x_1,x_2,\dots ,x_N;{\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)\\ & =L({\tilde{x}}_1,{\tilde{x}}_2,\dots ,{\tilde{x}}_N;{\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)\\ & =f\left({\tilde{x}}_1\right)f\left({\tilde{x}}_2\left|\left\{{\tilde{x}}_1\right\}\right)f\left({\tilde{x}}_3\right|\{{\tilde{x}}_2,{\tilde{x}}_1\}\right)\cdots f\left({\tilde{x}}_N|\{{\tilde{x}}_{N-1},{\tilde{x}}_{N-2},\dots ,{\tilde{x}}_1\}\right)\\ & =(2\text{π}{)}^{-\frac{Nmn}{2}}|{\mathrm{\Sigma }}_{mn}{|}^{-\frac{N}{2}}\exp \{-\frac{1}{2}\sum _{t=1}^N{({\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k})}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}({\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k})\},\end{array}$where $f(\cdot )$ means the probability density function, and we stipulate $\sum _{k=1}^0(\cdot )$ equals zero vector or zero matrix as needed. Therefore, the logarithm maximum likelihood function of $x_1,x_2,\dots ,x_N$ follows as (A10) $\begin{array}{rl}& \ell (x_1,x_2,\dots ,x_N;{\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)\\ & =\ln \left(L(x_1,x_2,\dots ,x_N;{\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)\right)\\ & =-\frac{Nmn}{2}\ln \left(2\text{π}\right)-\frac{N}{2}\ln \left(|{\mathrm{\Sigma }}_{mn}|\right)-\frac{1}{2}\sum _{t=1}^N{({\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k})}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}({\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k}).\end{array}$(A10) Using the derivative of scalar by matrix, it yields from (A10(A10) $\begin{array}{rl}& \ell (x_1,x_2,\dots ,x_N;{\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)\\ & =\ln \left(L(x_1,x_2,\dots ,x_N;{\mathrm{\Phi }}_1,\dots ,{\mathrm{\Phi }}_p,{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_p,{\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_q,{\mathrm{\Xi }}_1,\dots ,{\mathrm{\Xi }}_q)\right)\\ & =-\frac{Nmn}{2}\ln \left(2\text{π}\right)-\frac{N}{2}\ln \left(|{\mathrm{\Sigma }}_{mn}|\right)-\frac{1}{2}\sum _{t=1}^N{({\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k})}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}({\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k}).\end{array}$(A10) ) that (A11) $\{\begin{array}{rl}& \sum _{t=1}^N\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Phi }}_i}\left(H\right(t,G_{\cdot }{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}H(t,G_{\cdot }))=O_m,i=1,2,\dots ,p,\\ & \sum _{t=1}^N\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Psi }}_i}\left(H\right(t,G_{\cdot }{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}H(t,G_{\cdot }))=O_n,i=1,2,\dots ,p,\\ & \sum _{t=1}^N\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Theta }}_j}\left(H\right(t,G_{\cdot }{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}H(t,G_{\cdot }))=O_m,j=1,2,\dots ,q,\\ & \sum _{t=1}^N\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Xi }}_j}\left(H\right(t,G_{\cdot }{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}H(t,G_{\cdot }))=O_n,j=1,2,\dots ,q,\\ & N\frac{\mathrm{\partial }\ln \left(|{\mathrm{\Sigma }}_{mn}|\right)}{\mathrm{\partial }{\mathrm{\Sigma }}_{mn}}+\sum _{t=1}^N\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathrm{\Sigma }}_{mn}}\left(H\right(t,G_{\cdot }{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_{mn}^{-1}H(t,G_{\cdot }))=O_{mn},\end{array}$(A11) where $H(t,G_{\cdot })={\tilde{x}}_t+\sum _{k=1}^{t-1}{G}_k{\tilde{x}}_{t-k}$. It yields from Corollary 2.2 and Property 2.3 that $\{\begin{array}{rl}& \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Phi }}_i}(I_m\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_m,i=1,2,\dots ,p,\\ & \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Psi }}_i}(I_n\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_n,i=1,2,\dots ,p,\\ & \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Theta }}_j}(I_m\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_m,j=1,2,\dots ,q,\\ & \sum _{t=2}^N\sum _{k=1}^{t-1}\frac{\mathrm{\partial }{\left(G_k{\tilde{x}}_{t-k}\right)}^{\mathrm{\top }}}{\mathrm{\partial }{\mathrm{\Xi }}_j}(I_n\otimes \left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right])=O_n,j=1,2,\dots ,q,\\ & \frac{1}{N}\sum _{t=1}^N\mathrm{R}\mathrm{e}\mathrm{s}(\left[{\mathrm{\Sigma }}_{mn}^{-1}H\right(t,G_{\cdot }\left)\right]\otimes \left[{\left({\mathrm{\Sigma }}_{mn}^{-1}\right)}^{\mathrm{\top }}H\right(t,G_{\cdot }\left)\right],mn,mn)={\left({\mathrm{\Sigma }}_{mn}^{-1}\right)}^{\mathrm{\top }}.\end{array}$