# Identification and validation of periodic autoregressive model with additive
noise: finite-variance case

Wojciech Żuławiński<EMAIL_ADDRESS>Aleksandra Grzesiek
Radosław Zimroz  Agnieszka Wyłomańska  Faculty of Pure and Applied
Mathematics, Hugo Steinhaus Center, Wroclaw University of Science and
Technology, Wyspianskiego 27, 50-370 Wroclaw, Poland  Faculty of
Geoengineering, Mining and Geology, Wroclaw University of Science and
Technology, Na Grobli 15, 50-421 Wroclaw, Poland

###### Abstract

In this paper, we address the problem of modeling data with periodic
autoregressive (PAR) time series and additive noise. In most cases, the data
are processed assuming a noise-free model (i.e., without additive noise),
which is not a realistic assumption in real life. The first two steps in PAR
model identification are order selection and period estimation, so the main
focus is on these issues. Finally, the model should be validated, so a
procedure for analyzing the residuals, which are considered here as
multidimensional vectors, is proposed. Both order and period selection, as
well as model validation, are addressed by using the characteristic function
(CF) of the residual series. The CF is used to obtain the probability density
function, which is utilized in the information criterion and for residuals
distribution testing. To complete the PAR model analysis, the procedure for
estimating the coefficients is necessary. However, this issue is only
mentioned here as it is a separate task (under consideration in parallel). The
presented methodology can be considered as the general framework for analyzing
data with periodically non-stationary characteristics disturbed by finite-
variance external noise. The original contribution is in the selection of the
optimal model order and period identification, as well as the analysis of
residuals. All these findings have been inspired by our previous work on
machine condition monitoring that used PAR modeling

###### keywords:

periodic autoregressive model , additive noise , model identification , model
validation , autocovariance function , Monte Carlo simulations

††journal: arXiv

##  1 Introduction

In this paper, we propose methods for identifying periodic autoregressive
(PAR) models in data that consist of pure PAR time series and additive noise.
Our focus is on the case of finite-variance, meaning that both the innovation
series and the additive noise have finite variance. The PAR time series is
considered one of the most common models with periodically correlated (PC)
property, which is often observed in systems that have a periodically non-
stationary random structure [1, 2]. It has found many interesting
applications; see, e.g., hydrology [3], climatology and meteorology [4],
economics [5], and condition monitoring [6, 7, 8, 9]. In the classical
literature, mostly the pure (noise-free) PAR model is discussed. However, this
situation seldom occurs in real-world applications, and the measured signal is
almost always corrupted by noise. In practice, the additive noise can be
Gaussian, contain outliers, or even be leptokurtic. This requires specific
procedures that differ from those dedicated to a pure model. The research
studies on autoregressive (AR) model identification with additive noise with
different types of additional disturbances are well-known in the literature
[10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. In recent research, one can also
find studies on the PAR model with additive noise. However, most of the
attention is focused on the case where the additive noise consists of a
sequence of additive outliers, which are independent and identically
distributed (i.i.d.) random variables with large values that occur with a
certain probability. We recall here the bibliography references where new
estimation techniques are proposed for such models [20, 21, 22, 23, 24].
However, these algorithms do not account for the presence of additive
disturbances and instead use robust estimators of classical statistics to
minimize the effect of large observations that can be seen in the signal. The
special case, namely the PAR(1) model with general finite-variance additive
noise, is discussed in [25]. Therefore, it can be concluded that the
identification of the PAR model with general additive noise is rarely
addressed in the literature and there is a space for developing new and
effective approaches. This aspect has significant practical importance

The procedure of identifying a PAR model can be structured into several
stages: selecting the model order, determining the period, estimating the
model coefficients, and conducting a residual analysis to validate the model.
In certain situations, the model order and period may be advised or provided
by other sources. The most challenging aspect is finding the values of the
coefficients to complete the model, which is a serious problem from both a
theoretical and algorithmic perspective and has been discussed separately. In
this paper, we briefly recall the modified errors-in-variables method (see D).
The final step is to validate the model by analyzing its residuals, which
could also be considered as identifying an unknown distribution.

As mentioned above, the initial step in analyzing real data is selecting the
optimal order of the model and determining its period (relevant for PC
models). For pure models, the model order identification problem can be solved
using commonly known information criteria, such as $BIC$ and $AIC$ (see [26,
27]), which are calculated under the assumption of a specific distribution of
the residual series. However, when additive noise is present, these classical
approaches may not be sufficient. In this paper, the proposed methods for
optimal model order and period identification are also based on information
criteria. However, unlike classical approaches, the residuals of the analyzed
noise-corrupted PAR model are dependent. To overcome this issue, we analyze
the residuals as $T$-dimensional vectors, where $T$ is the period of the PAR
model. This arrangement in vectors allows us to consider the residuals as
i.i.d $T$-dimensional random variables with dependent components. Then, when
we assume that both the innovation series of the PAR model and the additive
noise are Gaussian, it can be deduced that the residuals also have a
$T$-dimensional Gaussian distribution with appropriate parameters. However,
the situation becomes much more complicated when the additive noise is not
Gaussian, as is the case in some real-world applications.

In such a case, we propose to apply the inverse of the characteristic function
(CF) for the residuals to obtain their probability density function (PDF),
which will be used in the information criterion. The theoretical CF is
provided in this paper for $T$-dimensional vectors representing the residuals
under the general assumption of a finite-variance PAR model with additive
noise. The form of the CF is also utilized in the final step of the data
analysis, which is the validation of the model. We propose an approach to test
the distribution of the residuals by comparing the empirical CF of the
residual series with the theoretical CF of the tested distribution. A similar
methodology has been applied, for instance, in [28] for sub-Gaussian
distributions.

The proposed algorithms were confirmed to be efficient through extensive Monte
Carlo simulations for two cases of additive noise distribution: Gaussian and a
mixture of Gaussian distributions (which represents the leptokurtic class of
distributions). The optimal order identification was considered in two
scenarios, with known and unknown period $T$. In both cases, the performed
simulations clearly indicate the efficiency of the proposed methodology. We
note that the presented methodology for order and period identification,
estimation of the parameters, and validation of the model by residual
analysis, can be considered as a general framework for analysis of the data
with periodically non-stationary characteristics disturbed by general finite-
variance external noise. We believe that this methodology will find
applications in many areas of interest.

The rest of the paper is organized as follows. In Section 2, we introduce the
noise-corrupted finite-variance PAR model. In Section 3, we describe the
procedure for the identification of the optimal order and period for the
considered model. Next, in Section 4, we introduce a procedure for the
residual analysis to validate the fitted model. In Section 5, we present a
simulation study for the optimal $p$ and $T$ value selection and for testing
the distribution of residuals. The last section concludes the paper. Note that
in A we demonstrate the equivalent matrix representation of the residual
series. In B, we recall the class of mixture of Gaussian distributions, while
in C we present the alternative proof for residuals distribution in the case
where the additive noise is a mixture of Gaussian distributions. Finally, in
D, we provide more details connected with the third step of the mentioned
procedure, i.e., estimating the model coefficients. However, a more detailed
study regarding this subject is prepared as a separate material.

##  2 Finite-variance periodic autoregressive model with additive noise

The periodic autoregressive model of order $p$ with additive noise (also
called noise-corrupted PAR($p$) model) is defined as follows:

$\displaystyle Y_{t}=X_{t}+Z_{t},$ (1)

where $\\{X_{t}\\}_{t\in\mathbb{Z}}$ is a periodic autoregressive (PAR) time
series of order $p\in\mathbb{N}$ and period $T\in\mathbb{N}$ defined below.
For simplicity, we assume that $p<T$. However, the case $p\geq T$ is briefly
treated as well. We assume that the innovations of the PAR($p$) model (denoted
further as $\\{\xi_{t}\\}_{t\in\mathbb{Z}}$) and the additive disturbances
$\\{Z_{t}\\}_{t\in\mathbb{Z}}$ are independent sequences of i.i.d. random
variables with mean zero and variances $\sigma_{\xi}^{2}$ and
$\sigma_{Z}^{2}$, respectively. In this paper, we consider the additive noise
as a sample of i.i.d. random variables from a continuous distribution with PDF
$f_{Z}(\cdot)$. However, the presented methodology is universal and can also
be applied in the case where the sequence $\\{Z_{t}\\}$ is derived from a
discrete distribution. In this case, the additive disturbances are considered
as additive outliers [20].

######  Definition 1.

[29] A sequence $\\{X_{t}\\}$ is a finite-variance PAR($p$) model with period
$T\in\mathbb{N}$ if it satisfies the following equation:

$\displaystyle X_{t}-\phi_{1}(t)X_{t-1}-\cdots-\phi_{p}(t)X_{t-p}=\xi_{t}.$
(2)

In Eq. (2), it is assumed that the parameters
$\phi_{1}(t),\phi_{2}(t),\cdots,\phi_{p}(t)$ are periodic in $t$ with the same
period $T$. It is common to consider $\\{\xi_{t}\\}$ to be a sequence of
Gaussian distributed random variables. In simulation studies, we also consider
this case. However, the presented methodology is also valid for any finite-
variance distribution of the innovation series. We mention that, similarly to
the PAR($p$) time series, the noise-corrupted model $\\{Y_{t}\\}$ exhibits the
periodically correlated (PC) property, which means that its mean and
autocovariance functions are periodic with respect to $t$ with the same period
$T$ [2]. More precisely, we can show that

$\displaystyle\mathbb{E}Y_{t}=\mathbb{E}Y_{t+T},\leavevmode\nobreak\
\leavevmode\nobreak\
cov(Y_{t},Y_{t+h})=cov(Y_{t+T},Y_{t+h+T}),\leavevmode\nobreak\
\leavevmode\nobreak\ h\in\mathbb{Z}.$ (3)

For a detailed proof of the special case, namely for the model with $p=1$, we
refer the readers to [25].

One can easily show that for each $t\in\mathbb{Z}$ the time sequence
$\\{Y_{t}\\}$ satisfies the following equation:

$\displaystyle
Y_{t}-\phi_{1}(t)Y_{t-1}-\cdots-\phi_{p}(t)Y_{t-p}=\xi_{t}+Z_{t}-\phi_{1}(t)Z_{t-1}-\cdots-\phi_{p}(t)Z_{t-p}.$
(4)

The residuals $\\{R_{t}\\}$ of the model given in Eq. (1) are as follows:

$\displaystyle R_{t}=Y_{t}-\phi_{1}(t)Y_{t-1}-\cdots-\phi_{p}(t)Y_{t-p}.$ (5)

Obviously, in the considered case, the residuals are not independent. However,
by defining the $T$-dimensional time series

$\displaystyle\mathbf{R}_{n}=[R_{nT+1},\cdots,R_{(n+1)T}],\leavevmode\nobreak\
\leavevmode\nobreak\ n=1,2,\cdots,N-1,$ (6)

one obtains the vectors of $T$-dimensional i.i.d. random variables (with
dependent components). In A, the matrix representation of the residual series
is presented. This representation is equivalent to the one given in this
section and can be useful e.g. for the identification of the residuals
distribution. One can show that
$\mathbb{E}(\mathbf{R}_{n}\mathbf{R}_{m}^{\prime})=0$, $n\neq m$ ($A^{\prime}$
denotes the transposition of the matrix $A$). Moreover, the covariance matrix
${\mathbf{\Gamma}^{R}}=[{\gamma^{R}}(k,l)]_{k,l=1}^{T}$ of $\mathbf{R}_{n}$ is
independent on $n$ and the coefficients
${\gamma^{R}}(k,l)=\mathbb{E}(R_{nT+k}R_{nT+l})$ are given by

$\displaystyle{\gamma}^{R}(k,l)=\left\\{\begin{array}[]{ll}\sigma_{\xi}^{2}+\sigma_{Z}^{2}\sum_{j=0}^{p}\phi^{2}_{j}(k)&\textrm{if
$k=l$},\\\
\sigma_{Z}^{2}\sum_{j=0}^{p+k-l}\phi_{j}(k)\phi_{j+l-k}(l)&\textrm{if
$k<l$},\\\
\sigma_{Z}^{2}\sum_{j=0}^{p+l-k}\phi_{j}(l)\phi_{j+k-l}(k)&\textrm{if
$k>l$},\end{array}\right.$ (10)

assuming that $\phi_{0}(k)=-1$ for any $k=1,2,\cdots,T$. In the case where the
sequences $\\{\xi_{t}\\}$ and $\\{Z_{t}\\}$ are Gaussian distributed, the PDF
of $\mathbf{R}_{n}$ for any $n$ is given by

$\displaystyle
f_{\mathbf{R_{n}}}(\mathbf{r})=\frac{\exp\left\\{-\frac{1}{2}\mathbf{r}(\mathbf{\Gamma}^{R})^{-1}\mathbf{r}^{\prime}\right\\}}{\sqrt{(2\pi)^{T}|{\mathbf{\Gamma}^{R}|}}},$
(11)

where $\mathbf{r}=[r_{1},\cdots,r_{T}]$. The fact that in this case the
$\mathbf{R}_{n}$ is multivariate Gaussian with the PDF given in Eq. (11) can
also be derived from a matrix representation presented in A. In the case where
$p\geq T$, one can consider the same methodology as above, but for
$KT$-dimensional residual blocks $\mathbf{R}_{n}$, where $K$ is the least
natural number for which $p<KT$. In Section 4 we also present how to find the
distribution of the residual series of the model (1) if the sequence
$\\{Z_{t}\\}$ is not Gaussian distributed. Let us note that the distribution
of the $T$-dimensional sequence $\\{\mathbf{R}_{n}\\}$ is crucial for the
validation of the model.

##  3 Identification of the optimal $p$ and $T$ parameters

For identifying the optimal order $p$ and period $T$ of the PAR($p$) model
with additive noise given in Eq. (1), we propose to use the Bayesian
information criterion ($BIC$), which is the estimator of the prediction error
and, thereby, the relative quality of statistical models for a given set of
data. However, other well-known criteria (such as $AIC$) could also be used in
a similar manner. Given a collection of fitted models for real data, $BIC$
estimates the quality of each model relative to other models. First, we
present the case where the period $T$ of the PAR model is known. The idea of
selecting the optimal order of the model $p^{*}_{opt}$ in this case is as
follows. For each order $p^{*}=1,2,\cdots,p_{max}$ ($p_{max}<T$), we estimate
the parameters of the model given in Eq. (1). Here, we propose to use the
estimation algorithm presented in D. Then we calculate the $BIC$ statistic,
defined as

$\displaystyle BIC=-2log(\mathcal{L})+log(\mathcal{N})\mathcal{K},$ (12)

where $\mathcal{L}$ is the likelihood function, $\mathcal{K}$ is the number of
the estimated parameters, and $\mathcal{N}$ is the number of observations (see
[26, 27]). Finally, the optimal order is equal to the $p_{opt}^{*}$ that
minimizes the $BIC$ statistic in Eq. (12).

For the considered case, the residuals of the model, see Eq. (5), do not
constitute a sample of independent random variables. Thus, the $BIC$ statistic
is defined here for the series of $T$-dimensional i.i.d. random variables
defined in Eq. (6). If we assume that for each $n$, $\mathbf{R}_{n}$ has PDF
$f_{\mathbf{R}_{n}}(\cdot)$, then the $BIC$ statistic takes the form

$\displaystyle
BIC(p^{*},T)=-2\sum_{n=1}^{N-1}log(f_{\mathbf{R}_{n}}(\mathbf{r}_{n}))+log(NT)(Tp^{*}+2),$
(13)

where $\mathbf{r}_{n}$ is the $T$-dimensional residual vector corresponding to
$\mathbf{R}_{n}$ and $NT$ is the number of observations of the model (1). Note
that starting from $n=1$ means that we actually take into account the
truncated sequence of residuals $[R_{T+1},\cdots,R_{NT}]$. If $\\{Z_{t}\\}$ is
a sequence of Gaussian distributed random variables, then the PDF used in Eq.
(13) is given by Eq. (11). If the sequence $\\{Z_{t}\\}$ is not Gaussian, we
propose to use the methodology based on the inverse CF, see Section 4 for more
details.

In case where the period $T$ is unknown, we can also apply the $BIC$ statistic
defined in Eq. (13); however, we consider it as a function of two arguments,
$p^{*}$ and $T^{*}$. We calculate its value for all combinations of
$p^{*}=1,\cdots,p_{max}$ and $T^{*}=1,\cdots,T_{max}$ such that $p^{*}<T^{*}$.
In the final stage, we select the optimal $p_{opt}^{*}$ and $T_{opt}^{*}$ such
that they minimize the $BIC$ statistic. The other approach is to apply one of
the known methods for selecting the period $T$ and then apply the above-
described methodology assuming that the period is known. For other methods of
$T$ selection, we refer the readers, for instance, to [5, 2]. In the mentioned
bibliography positions, the authors discuss efficient methods for period
identification for PC models without additive disturbances. However, when the
variance of $\\{Z_{t}\\}$ is relatively small with respect to the pure model,
the classical approaches may also be used. The methodology for identifying the
period for the general case will be considered in future studies. Let us note
that the methodology described above can also be used for the simultaneous
order/period identification for the pure PAR model given in Eq. (2).

##  4 Validation of the model - residual analysis

In time series analysis, one of the most common techniques for model
validation is residual analysis. We recall that the residuals of the model
given by Eq. (5) are not independent. Thus, we examine the $T$-dimensional
residual series corresponding to independent $T$-dimensional random vectors
$\mathbf{R}_{n}$ given by Eq. (6). The first step of the analysis is to
confirm whether the $T$-dimensional blocks of the residuals constitute a
sample of independent random vectors. For this purpose, we propose to apply
one of the tests verifying the independence of multidimensional vectors; see,
e.g., [30, 31, 32]. The next step is to test whether the residual series comes
from the assumed $T$-dimensional distribution. The situation is relatively
easy if both the innovation sequence $\\{\xi_{t}\\}$ and the additive noise
$\\{Z_{t}\\}$ are Gaussian distributed, as then the residual series is also
$T$-dimensional Gaussian. In this case, to test the $T$-dimensional Gaussian
distribution, we refer the reader to the review of the new developments in
affine-invariant tests for multivariate Gaussianity [33]. The situation is
less standard when the additive noise is not Gaussian, so that the residual
series is not $T$-dimensional Gaussian as well. To overcome this issue, in the
following part of this section, we propose a testing procedure that can be
used in case of any distribution of the additive noise. It is based on the CF
of the residual series, specified in Theorem 1.

######  Theorem 1.

The characteristic function of $T$-dimensional residual series
$\mathbf{R}_{n}=[R_{nT+1},\cdots,R_{(n+1)T}]$ (see Eq. (5) and Eq. (6)) of the
noise-corrupted PAR($p$) model with period $T$, takes the form

$\Phi_{\mathbf{R}_{n}}(t_{1},\cdots,t_{T})=\prod_{j=1}^{T}\Phi_{\xi}(t_{j})\,\prod_{j=1}^{p+T}\Phi_{Z}\left(-\sum_{i=\max\\{1,j-p\\}}^{\min\\{j,T\\}}t_{i}\phi_{p-(j-i)}(i)\right),$
(14)

where the $\Phi_{\xi}(\cdot)$ and $\Phi_{Z}(\cdot)$ are the CFs of the
$\\{\xi_{t}\\}$ and $\\{Z_{t}\\}$, respectively, and $\phi_{0}(k)=-1$ for
$k=1,\cdots,T$.

######  Proof.

Using Eq. (5) and the fact that the parameters of the PAR(p) model are
periodic with period $T$ we can write that

$\displaystyle\Phi_{\mathbf{R}_{n}}(t_{1},\cdots,t_{T})=$
$\displaystyle\mathbb{E}\left[\exp\\{it_{1}R_{nT+1}+it_{2}R_{nT+2}+\cdots+it_{T}R_{nT+T}\\}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\exp\\{it_{1}(\xi_{nT+1}+Z_{nT+1}-\phi_{1}(nT+1)Z_{(nT+1)-1}-\cdots-\phi_{p}(nT+1)Z_{(nT+1)-p})+$
$\displaystyle\quad\quad\quad
it_{2}(\xi_{nT+2}+Z_{nT+2}-\phi_{1}(nT+2)Z_{(nT+2)-1}-\cdots-\phi_{p}(nT+2)Z_{(nT+2)-p})+$
$\displaystyle\quad\quad\quad\cdots+$ $\displaystyle\quad\quad\quad
it_{T}(\xi_{nT+T}+Z_{nT+T}-\phi_{1}(nT+T)Z_{(nT+T)-1}-\cdots-\phi_{p}(nT+T)Z_{(nT+T)-p})\\}]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\exp\\{it_{1}(\xi_{nT+1}+Z_{nT+1}-\phi_{1}(1)Z_{nT}-\cdots-\phi_{p}(1)Z_{nT+1-p})+$
$\displaystyle\quad\quad\quad
it_{2}(\xi_{nT+2}+Z_{nT+2}-\phi_{1}(2)Z_{nT+1}-\cdots-\phi_{p}(2)Z_{nT+2-p})+$
$\displaystyle\quad\quad\quad\cdots+$ $\displaystyle\quad\quad\quad
it_{T}(\xi_{nT+T}+Z_{nT+T}-\phi_{1}(T)Z_{nT+T-1}-\cdots-\phi_{p}(nT+T)Z_{nT+T-p})\\}].$

Let us note that the argument of the exponential function under expected value
can be written as a linear combination of $2T+p$ elements, namely
$\xi_{nT+1},\xi_{nT+2},\cdots,\xi_{nT+T}$ and
$Z_{nT+1-p},Z_{nT+2-p},\cdots,Z_{nT+T}$, which are all independent random
variables. Using the above remark and the fact that

$\xi_{nT+1}\stackrel{{\scriptstyle d}}{{=}}\xi_{nT+2}\stackrel{{\scriptstyle
d}}{{=}}\cdots\stackrel{{\scriptstyle d}}{{=}}\xi_{nT+T}$

and

$Z_{nT+1-p}\stackrel{{\scriptstyle d}}{{=}}Z_{nT+2-p}\stackrel{{\scriptstyle
d}}{{=}}\cdots\stackrel{{\scriptstyle d}}{{=}}Z_{nT+T}$

one can write the formula given in Eq. (14). ∎

######  Example 1.

Let $\\{Y_{t}\\}$ denote the PAR model with additive noise defined in Eq. (1)
with Gaussian innovations $\\{\xi_{t}\\}$ and Gaussian additive noise
$\\{Z_{t}\\}$. Then, assuming $T=2$ and $p=1$ and using Eq. (14), we obtain
the following characteristic function for the residual vector
$\mathbf{R}_{n}=[R_{2n+1},R_{2n+2}]$:

$\displaystyle\Phi_{\mathbf{R}_{n}}(t_{1},t_{2})=$
$\displaystyle\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{Z}^{2}+\sigma_{Z}^{2})-2t_{1}t_{2}\phi(2)\sigma_{Z}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{Z}^{2}+\sigma_{Z}^{2})\right)\right),$
(15)

which is simply the characteristic function of a two-dimensional zero-mean
Gaussian random vector with the elements of the covariance matrix
$\bm{\Gamma}^{R}=[\gamma^{R}(k,l)]_{k,l=1}^{2}$ given in Eq. (10).

######  Example 2.

Let $\\{Y_{t}\\}$ denote the PAR model with additive noise defined in Eq. (1).
Moreover, let us assume that the innovations $\\{\xi_{t}\\}$ are Gaussian,
whereas the additive disturbances $\\{Z_{t}\\}$ are mixture of two Gaussian
distributions with parameters $a_{1},a_{2}$, $\sigma_{1}^{2},\sigma_{2}^{2}$
and $\mu_{1}=\mu_{2}=0$, see B for more details. Then, using Eq. (14) for
$T=2$ and $p=1$, we obtain the following characteristic function for the
residual vector $\mathbf{R}_{n}=[R_{2n+1},R_{2n+2}]$:

$\displaystyle\Phi_{\mathbf{R}_{n}}(t_{1},t_{2})$
$\displaystyle=a_{1}^{3}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{1}^{2}+\sigma_{1}^{2})-2t_{1}t_{2}\phi(2)\sigma_{1}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{1}^{2}+\sigma_{1}^{2})\right)\right)$
(16)
$\displaystyle+a_{1}^{2}a_{2}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{1}^{2}+\sigma_{1}^{2})-2t_{1}t_{2}\phi(2)\sigma_{1}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{1}^{2}+\sigma_{2}^{2})\right)\right)$
$\displaystyle+a_{1}^{2}a_{2}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{1}^{2}+\sigma_{2}^{2})-2t_{1}t_{2}\phi(2)\sigma_{2}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{2}^{2}+\sigma_{1}^{2})\right)\right)$
$\displaystyle+a_{1}a_{2}^{2}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{1}^{2}+\sigma_{2}^{2})-2t_{1}t_{2}\phi(2)\sigma_{2}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{2}^{2}+\sigma_{2}^{2})\right)\right)$
$\displaystyle+a_{1}^{2}a_{2}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{2}^{2}+\sigma_{1}^{2})-2t_{1}t_{2}\phi(2)\sigma_{1}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{1}^{2}+\sigma_{1}^{2})\right)\right)$
$\displaystyle+a_{1}a_{2}^{2}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{2}^{2}+\sigma_{1}^{2})-2t_{1}t_{2}\phi(2)\sigma_{1}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{1}^{2}+\sigma_{2}^{2})\right)\right)$
$\displaystyle+a_{1}a_{2}^{2}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{2}^{2}+\sigma_{2}^{2})-2t_{1}t_{2}\phi(2)\sigma_{2}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{2}^{2}+\sigma_{1}^{2})\right)\right)$
$\displaystyle+a_{2}^{3}\exp\left(-\frac{1}{2}\left(t_{1}^{2}(\sigma_{\xi}^{2}+\phi^{2}(1)\sigma_{2}^{2}+\sigma_{2}^{2})-2t_{1}t_{2}\phi(2)\sigma_{2}^{2}+t_{2}^{2}(\sigma_{\xi}^{2}+\phi^{2}(2)\sigma_{2}^{2}+\sigma_{2}^{2})\right)\right).$

It should be noted that the function above is the CF of a mixture of eight
two-dimensional Gaussian distributions. The distribution of the residual
vector in this case can also be obtained using the matrix representation
provided in A. Therefore, the PDF of the residual vector can be expressed as a
linear combination of eight zero-mean two-dimensional Gaussian probability
density functions. For further explanation, please refer to C.

The CF presented in Theorem 1 can be applied to the residuals distribution
testing that is crucial for the validation of the considered model. The
procedure is similar to the methodology presented in [28], where the authors
introduced a testing algorithm based on the distance between the empirical
(for a given vector of observations) and theoretical (from tested
distribution) CFs. The proposed test is especially useful when the tested
distribution does not rely on a closed-form PDF. In the mentioned bibliography
position, the methodology was applied to the bivariate sub-Gaussian
distributed random vectors; however, it is universal and can be extended to
any distribution with a closed-form CF, see also [34].

A sketch of the proposed test is as follows. We consider the following test
statistic:

$D=\sup_{t\in\mathbb{R}^{n}}\left|\hat{\Phi}_{\mathbf{R}_{n}}(\mathbf{t})-\Phi_{\mathbf{R}_{n}}(\mathbf{t})\right|,$
(17)

where $\hat{\Phi}_{\mathbf{R}_{n}}(\cdot)$ and
${\Phi}_{\mathbf{R}_{n}}(\cdot)$ denote the empirical and theoretical
characteristic functions of the residual series ${\mathbf{R}_{n}}$,
respectively. Moreover, let us formulate the null hypothesis ($H_{0}$) that
the vector of observations corresponds to the model defined in Eq. (1) with
specified distributions of the innovations and additive noise for specified
set of all parameters. In practice, those parameters are estimated using the
procedures presented in D and, in consequence, the specification of the null
hypothesis is done based on the parameters obtained from the numerical search.
The alternative hypothesis ($H_{a}$) states that the above model does not fit
the data considered. The testing procedure can be briefly described as
follows: (a) For the data under consideration estimate the parameters of PAR
model with additive noise and calculate the test statistic $D$ given in Eq.
(17); (b) Assume that the model fitted in (a) is the model under $H_{0}$; (c)
Simulate the data from the model specified in $H_{0}$; (d) For the simulated
data, calculate the residuals corresponding to $\mathbf{R}_{n}$ under $H_{0}$;
(e) Calculate the test statistic given in Eq. (17) - denote it as $D_{1}$; (f)
Repeat the points (c)-(e) $M$ times to obtain $D_{1},D_{2},\cdots,D_{M}$
values; (g) Calculate the p-value as a fraction of $D_{i}$’s, $i=1,\cdots,M$,
that exceed $D$; (h) If the p-value is less that the chosen significance
level, reject the $H_{0}$ hypothesis.

We mention here that formula (14) can be used to find the PDF of the residual
series by applying the inverse Fourier transform of CF [35]. Obviously, the
inverse Fourier transform is calculated numerically. In this paper, we use the
multivariate version of the fast Fourier transform-based algorithm, proposed
in [36]. The PDF is approximated on the grid consisting of $(2^{5})^{T}$
points, and its value for in-between points is computed using linear
interpolation. The resulting PDF is further used for the $BIC$ statistic (13)
to obtain the optimal $p$ and $T$ parameters. More details on that aspect are
presented in Section 5.1. 111 All algorithms used in this paper were prepared
in MATLAB R2021b and are available upon the special request of the reader.

##  5 Simulation study

###  5.1 Identification of the optimal $p$ and $T$ parameters

In this part, the methodology for the selection of model order $p$ and period
$T$ presented in Section 3 is assessed using Monte Carlo simulations with
sample trajectories generated from the noise-corrupted PAR model $\\{Y_{t}\\}$
given in Eq. (1). Here, we consider models with $p=1,2,3$ and $T=4$, with the
following matrices of PAR coefficients $\Phi^{(p)}$:

$\displaystyle\Phi^{(1)}=\begin{bmatrix}-0.1208\\\ -0.5773\\\ -0.0362\\\
-0.3254\\\ \end{bmatrix}\quad\Phi^{(2)}=\begin{bmatrix}-0.1208&-0.0878\\\
-0.5773&-0.9798\\\ -0.0362&0.9196\\\ -0.3254&-0.5802\\\
\end{bmatrix}\quad\Phi^{(3)}=\begin{bmatrix}-0.1208&-0.0878&0.6605\\\
-0.5773&-0.9798&-0.6826\\\ -0.0362&0.9196&0.6555\\\ -0.3254&-0.5802&-0.5313\\\
\end{bmatrix},$ (18)

where the element in the $v$-th row and $i$-th column is the $\phi_{i}(v)$
coefficient, for $v=1,\cdots,T$ and $i=1,\cdots,p$. In the whole simulation
study, we assume standard Gaussian innovations of PAR model $\\{\xi_{t}\\}\sim
N(0,\sigma_{\xi}^{2}=1)$.

The first experiment performed concerns the selection of the order $p$ in the
situation with the known period $T=4$ and Gaussian-distributed additive noise.
The procedure is as follows. We simulate $M=1000$ trajectories of length
$NT=1200$ from the PAR model of order $p$ (i.e. with coefficients
$\Phi^{(p)}$, see Eq. (18)) with additive noise $\\{Z_{t}\\}\sim
N(0,\sigma_{Z}^{2})$. Then, for each generated trajectory, the $BIC(p^{*},T)$
statistic is calculated for $p^{*}=1,2,3$ (see Eq. (13)). The $p^{*}_{opt}$ is
the order $p^{*}$ that minimizes $BIC(p^{*},T)$. In calculations, to obtain
residual blocks $\mathbf{r}_{n}$ for $n=1,\cdots,N-1$, we estimate the
coefficients of the PAR model using the modified errors-in-variables method
described in D. This experiment is performed for $p=1,\,2,\,3$ and for
$\sigma_{Z}^{2}=0.2,\,1,\,2$. The time series sample trajectories considered
in these simulations are presented in Fig. 1. The same sample of pure model
$\\{X_{t}\\}$ is taken for each $\sigma_{Z}^{2}$. Moreover, on the mentioned
figure, the empirical variances (calculated from 100000 trajectories) are also
plotted. Due to space limitations, we present only the samples for $p=1$, but
for $p=2,3$ the observed behaviour is similar. It can be seen that for
increasing additive noise variance, the underlying PAR trajectory becomes more
hidden. This is confirmed by empirical signal-to-noise ratios defined as

$\displaystyle SNR(t)=\frac{Var(X_{t})}{Var(Z_{t})},$ (19)

calculated for 100000 trajectories and illustrated in Fig. 2, now for each
$p=1,2,3$. Because with larger $\sigma_{Z}^{2}$ the pure PAR process is more
disturbed, the main challenge here is to provide reliable results also for low
signal-to-noise ratio samples.

Figure 1: Sample trajectories of $\\{Y_{t}\\}$ time series and empirical
variances of $Y_{t}$ for each $t$ calculated from 100000 trajectories for
different values of Gaussian additive noise variance in the case of
$\Phi^{(1)}$ coefficients ($p=1$, $T=4$, see Eq. (18)).

Figure 2: Empirical signal-to-noise ratios (see Eq. (19)) calculated for
100000 trajectories of $\\{Y_{t}\\}$ time series with coefficients
$\Phi^{(p)}$ (see Eq. (18)) and different values of $\sigma_{Z}^{2}$.

The boxplots of all $BIC$ values in these simulations set are presented in
Fig. 3. Let us first analyze the results from the top panel, i.e. for the
lowest considered additive noise variance $\sigma_{Z}^{2}=0.2$. One can see
that for each $p$ the results for $p^{*}<p$ (if any) are visibly larger than
the remaining ones. Furthermore, starting from $p^{*}=p$, the values of $BIC$
seem to stabilize. Both of these characteristics are present in the case of
$p=2$. This behaviour is caused by the fact that for $p^{*}>p$ the model can
also be adjusted properly, so that the estimated coefficients
$\\{\phi_{i}(v)\\}$ for $i=1,\cdots,p$ are close to the true ones and for
$i=p+1,\cdots,p^{*}$ are nearly zero. From the presented results, one can
conclude that instead of the $BIC$ minimization, one can consider another
approach of order identification, i.e., detecting the order $p^{*}$ for which
the described stabilization behaviour starts to occur. For the case of
$\sigma_{Z}^{2}=1$ (middle panel), although the results for all $p^{*}$ are
now relatively closer to each other, one can see a pattern similar to that
before. In particular, one can still distinguish between the results for
$p^{*}<p$ and $p^{*}\geq p$. However, when additive noise variance is even
larger (see the bottom panel of Fig. 3, where $\sigma_{Z}^{2}=2$), the
difference between boxplots for different $p^{*}$ becomes much less visible.
Hence, in this case, the order identification becomes more challenging, as
expected, considering the lower signal-to-noise ratio.

Figure 3: Boxplots of $BIC(p^{*},T)$ values for checked orders $p^{*}=1,2,3$
calculated in case of known $T=4$ for true model orders $p=1,2,3$ and
$\\{Z_{t}\\}\sim N(0,\sigma_{Z}^{2})$ with $\sigma_{Z}^{2}=0.2,1,2$. The
shaded area on each subplot refers to the boxplot for $p^{*}=p$.

We also calculate the fraction of correctly identified orders by $BIC$
minimization, i.e. the percentage of cases where $p^{*}_{opt}=p$, for each
considered combination of $p$ and $\sigma_{Z}^{2}$. The results are presented
in Table 1. Most of all, we can see that the identification efficiency indeed
drops for a larger additive noise variance. However, in general, the results
show that in the vast majority of cases the proposed criterion managed to
identify the order correctly, even for $\sigma_{Z}^{2}=2$. Let us recall that
the order selection here was done by the minimization of $BIC$, despite the
stabilization behaviour mentioned above, which could result in a wrong
selection (i.e., $BIC$ would tend to choose larger order than the correct
one).

| $\\{Z_{t}\\}\sim N(0,\sigma_{Z}^{2})$ | $\\{Z_{t}\\}=\sigma_{Z}\cdot\\{\tilde{Z}_{t}\\}$
---|---|---
| $p=1$ | $p=2$ | $p=3$ | $p=1$ | $p=2$ | $p=3$
$\sigma_{Z}^{2}=0.2$ | 98.7% | 99.9% | 100% | 98.5% | 100% | 100%
$\sigma_{Z}^{2}=1$ | 96% | 96.8% | 99.2% | 93.9% | 97.8% | 99.2%
$\sigma_{Z}^{2}=2$ | 91.8% | 86.1% | 73.0% | 90.2% | 85.8% | 76.6%
Table 1: Fraction of correctly identified orders in case of known $T=4$ for
different $p$, $\sigma_{Z}^{2}$ and both considered additive noise
distributions. Here $\\{\tilde{Z}_{t}\\}$ is a sequence of mixture Gaussian
distributed random variables with $a_{1}=a_{2}=0.5$, $\mu_{1}=\mu_{2}=0$ and
$\sigma_{1}^{2}=0.5$, $\sigma_{2}^{2}=1.5$.

Next, let us present the same experiment as the one described above but for
another type of additive noise distribution. Let us consider a sequence
$\\{\tilde{Z}_{t}\\}$ of i.i.d. random variables from a mixture of $m=2$
Gaussian distributions (see B) with $a_{1}=a_{2}=0.5$, $\mu_{1}=\mu_{2}=0$ and
$\sigma_{1}^{2}=0.5$, $\sigma_{2}^{2}=1.5$. These parameters are tuned in such
a way that $\tilde{Z}_{t}$ for each $t$ has a unit variance. Furthermore, this
mixture yields a heavy-tailed distribution as its excess kurtosis (see B) is
positive. As the additive noise sequences used in simulations described below
we use scaled versions of $\\{\tilde{Z}_{t}\\}$, setting
$Z_{t}=\sigma_{Z}\cdot\tilde{Z}_{t}$, again taking $\sigma_{Z}^{2}=0.2,1,2$.
For signal-to-noise ratios for this case, see again Fig. 19. In this case, to
present the proposed methodology, we use the procedure presented in Section 4
based on the CF inversion. However, it should be noted that the procedure of
inverting the CF could be ignored here, since the PDF of the residuals is
given in the explicit form, see C for more details. In this part, the whole
procedure is presented for illustration of the introduced methodology for the
case where the additive noise has different than the Gaussian distribution.

Let us note that in this case the estimation stage provides information about
$\sigma_{Z}^{2}$ (in the form of its estimate $\hat{\sigma}_{Z}^{2}$), but not
about the structure of the mixture (i.e., the form of $\tilde{Z}_{t}$) that is
also needed in the calculation of the CF. Hence, here we assume that this
structure is known and only the additive noise variance may vary dependent on
estimation results. In other words, during the CF computation for a trajectory
with estimated additive noise variance $\hat{\sigma}_{Z}^{2}$, we assume that
$Z_{t}=\hat{\sigma}_{Z}\cdot\tilde{Z}_{t}$.

The $BIC$ statistic calculated for trajectories with additive noise from
mixture Gaussian distribution (all other details of these simulations are the
same as for the Gaussian case) are presented in Fig. 4. One can see that
observed characteristics of the results are similar to those obtained for the
Gaussian case described above. Most importantly, once again the aforementioned
stabilization for $p^{*}\geq p$ occurs, and the results for different $p^{*}$
become less distinguishable with increasing additive noise variance. The
latter can be again seen in fractions of correctly identified orders which are
presented in Table 1. Let us note that these results are at a similar level as
the proportions obtained for samples with Gaussian additive noise, even though
here we only use an approximation of the PDF of the residuals. Hence, one can
see that the presented approach shows an acceptable efficiency even for non-
Gaussian cases.

Figure 4: Boxplots of $BIC(p^{*},T)$ values for checked orders $p^{*}=1,2,3$
calculated in case of known $T=4$ for true model orders $p=1,2,3$ and
$\\{Z_{t}\\}=\sigma_{Z}\cdot\\{\tilde{Z}_{t}\\}$ with
$\sigma_{Z}^{2}=0.2,1,2$, where $\\{\tilde{Z}_{t}\\}$ is a sequence of mixture
Gaussian distributed random variables with $a_{1}=a_{2}=0.5$,
$\mu_{1}=\mu_{2}=0$ and $\sigma_{1}^{2}=0.5$, $\sigma_{2}^{2}=1.5$. The shaded
area on each subplot refers to the boxplot for $p^{*}=p$.

In the last part of the simulation study, let us consider the problem of
identification of both, the order $p$ and the period $T$, simultaneously,
using the $BIC(p^{*},T^{*})$ statistic. Here, we only analyze the case with
Gaussian additive noise. The simulations procedure performed here, similar to
the ones presented before, is as follows. We generate $M=1000$ trajectories of
length $1205$ of the PAR model with $\Phi^{(2)}$ coefficients (thus, $p=2$ and
$T=4$, see Eq. (18)) and additive noise $\\{Z_{t}\\}\sim N(0,\sigma_{Z}^{2})$.
Next, for each simulated sample path, we calculate the $BIC(p^{*},T^{*})$
statistic for all pairs $p^{*},T^{*}\in\\{1,2,3,4,5\\}$ such that
$p^{*}<T^{*}$. We also identify $p^{*}_{opt}$ and $T^{*}_{opt}$, i.e. the pair
that minimizes the computed criterion. Such choice of sample length is made to
ensure that for all checked potential periods the sequence of residuals
considered in $BIC$ calculation has the same length. Here, for each $T^{*}$,
we use the sequence $[R_{T^{*}_{max}+1},\cdots,R_{L}]$, where $T^{*}_{max}=5$,
which here is of length 1200 (note that it is divisible by each considered
$T^{*}$). As before, we consider three levels of additive noise variance,
namely $\sigma_{Z}^{2}=0.2,1,2$. The calculated $BIC$ results are illustrated
in Fig. 5. First, let us analyze the case where $\sigma_{Z}^{2}=0.2$. One can
clearly see that the lowest values were obtained for $p^{*}=2$ and $T^{*}=4$,
which is the actual order/period pair of the simulated model, and for
$p^{*}=3$, $T^{*}=4$, so for the true period and the order larger than the
true one. Let us note that this is consistent with the stabilization behaviour
observed in the case of known $T$. Moreover, since the results for $T^{*}\neq
T$ are at a significantly high level so that the $BIC$ statistic would not
select them as the period of a model, one can assume that this criterion can
serve as a period identification method for noise-corrupted PAR models.
However, as before, for larger values of additive noise variance, the $BIC$
results for all considered $p^{*},T^{*}$ pairs become relatively closer to
each other so that the proper detection of order and period is then more
challenging. Although for $\sigma_{Z}^{2}=1$ the $p^{*}=2$, $T^{*}=4$ pair
still stands out, for $\sigma_{Z}^{2}=2$ all boxplots are in fact at a similar
level.

Figure 5: Boxplots of $BIC(p^{*},T^{*})$ values for all checked pairs
$p^{*},T^{*}$ calculated for $\\{Z_{t}\\}\sim N(0,\sigma_{Z}^{2})$ with
$\sigma_{Z}^{2}=0.2,1,2$ in case of unknown $T$. The shaded area on each
subplot refers to the boxplot for $p^{*}=p=2$, $T^{*}=T=4$.

We calculate the fraction of correctly identified values of $p$ and $T$, i.e.
the percentage of cases where both $p^{*}_{opt}=p$ and $T^{*}_{opt}=T$.
Moreover, we also consider the proportion of trajectories for which just the
latter equality is satisfied (regardless of the selected order). Let us note
that this result indicates the performance of the presented methodology in
period detection. The percentages considered for different values
$\sigma_{Z}^{2}$ are presented in Table 2. As before, the model identification
efficiency drops with increasing additive noise variance from an almost
perfect performance for $\sigma_{Z}^{2}=0.2$ to a significantly worse one for
$\sigma_{Z}^{2}=2$. This can be seen for both simultaneous selections of $p$
and $T$ selection as well as for just the identification of the period.
However, let us note that the considered $BIC$ statistic may sometimes be able
to find the correct $T$ even if the order obtained is wrong. This is most
clearly seen on the difference between both obtained fractions for
$\sigma_{Z}^{2}=1$. In this case, even though the performance of both $p,T$
selection significantly decreases in comparison to the $\sigma_{Z}^{2}=0.2$
situation, for the period identification it is still at a very high level.

|  $p^{*}_{opt}=p$, $T^{*}_{opt}=T$ | $T^{*}_{opt}=T$
---|---|---
$\sigma_{Z}^{2}=0.2$ | 98.9% | 100%
$\sigma_{Z}^{2}=1$ | 74.1% | 98.6%
$\sigma_{Z}^{2}=2$ | 37.5% | 50.7%
Table 2: Fraction of correctly identified parameters for different
$\sigma_{Z}^{2}$ with Gaussian additive noise $\\{Z_{t}\\}$ in case of unknown
$T$.

###  5.2 Testing the distribution of residuals

In this part, we present the results of the residual series testing
methodology outlined in detail in Section 4. In the simulations conducted, we
focus on assessing the test power, which is considered here as a function of
the additive noise variance. Below, we provide a brief overview of the
verification process carried out in this section. The procedure is as follows:
under the null hypothesis $H_{0}$, we assume a certain PAR model of order $1$
with $T=2$, zero-mean, unit variance Gaussian innovations, and a certain type
of additive noise. We then simulate $M=1000$ trajectories from that model for
different levels of additive noise variance (not necessarily equal to the
level assumed in $H_{0}$, but the other parameters do not change). For these
data, we perform the proposed test and calculate the percentage of repetitions
where $H_{0}$ is rejected. The above-mentioned simulations are conducted for a
PAR(1) model with $T=2$ ($\phi(1)=0.4$, $\phi(2)=-0.6$), Gaussian innovations
($\sigma^{2}_{\xi}=1$), and two types of additive noise distributions:
Gaussian ($\sigma_{Z}^{2}=1$ for $H_{0}$) and a mixture of two Gaussians
($a_{1}=a_{2}=0.5$, $\mu_{1}=\mu_{2}=0$, $\sigma_{1}^{2}=0.5$, and
$\sigma_{2}^{2}=1.5$ for $H_{0}$). The results are shown in Fig. 6, where on
the $y$-axis we present the power of the test for different values of the
additive noise variance given on the $x$-axis. It should be noted that
$\sigma^{2}_{\xi}=1$ corresponds to the null hypothesis $H_{0}$, and in this
case, the power of the test is close to the assumed significance level (here
equal to $0.05$). As can be seen, the test power takes slightly higher values
for the same level of additive noise variance in the Gaussian case than for
the mixture of Gaussian distributions. The test has very low power when the
variance of the additional disturbances is close to the one assumed in
$H_{0}$. We can see that the power of the test is close to the significance
level when the variance of the additional noise differs by a maximum of $\pm
0.2$ (for the Gaussian) or $\pm 0.3$ (for the mixture Gaussian) from the
variance assumed in $H_{0}$. The power of the test increases as the distance
of the additional noise variance from the level assumed in $H_{0}$ becomes
greater. However, we can see that the test power functions shown in Fig. 6 are
not symmetrical around 1, and the power is greater for the smaller values of
the actual variance of the additional disturbances, so we conclude that it is
easier to recognize the wrong model in this case. It is important to note that
the results presented in this section are calculated using a bivariate
argument $t$ for the test statistic, as defined by Eq. (17), that was chosen
from a square grid of values ranging from $-10$ to $10$ in both dimensions
with a step size of $0.01$.

Figure 6: The power of the test for different levels of additive noise
variance. In $H_{0}$ we assume PAR(1) model with $T=2$ ($\phi(1)=0.4$,
$\phi(2)=-0.6$), Gaussian innovations ($\sigma^{2}_{\xi}$=1) and the additive
noise on the left panel is Gaussian ($\sigma_{Z}^{2}=1$) and on the right
panel is the mixture of two Gaussian distributions ($a_{1}=a_{2}=0.5$,
$\mu_{1}=\mu_{2}=0$, $\sigma_{1}^{2}=0.5$, and $\sigma_{2}^{2}=1.5$).

##  6 Conclusions

In this article, we present a new framework for identifying and validating the
noise-corrupted PAR model where the data are a mixture of pure PAR time series
and additive noise. This scenario has been observed in some real-world data,
particularly in the context of machine condition monitoring. We have focused
on the procedures for selecting the model order and estimating the period for
the assumed model of observation. Finally, a novel method for analyzing
residuals is presented. The proposed methodology is universal because it is
independent of the distribution of the innovation series of the PAR model and
the noise distribution. The efficiency is demonstrated for two types of
additive noise: Gaussian and finite mixtures of Gaussian distributions. The
common feature of the introduced methods is the form of the CF of the residual
series (presented here as a $T$-dimensional random vector). We believe that
the proposed procedures can be applied to various data with expected
periodically non-stationary characteristics that are naturally affected by
external sources.

## Acknowledgments

The work of WŻ, RZ and AW is supported by National Center of Science under
Sheng2 project No. UMO-2021/40/Q/ST8/00024 "NonGauMech - New methods of
processing non-stationary signals (identification, segmentation, extraction,
modeling) with non-Gaussian characteristics for the purpose of monitoring
complex mechanical structures".

## Conflicts of interest

Declarations of interest: none

## References

  * [1] A. Napolitano, Cyclostationarity: New trends and applications, Signal Process. 120 (2016) 385–408.
  * [2] H. L. Hurd, A. Miamee, Periodically correlated random sequences: Spectral theory and practice, Vol. 355, John Wiley & Sons, 2007.
  * [3] D. Bukofzer, Optimum and suboptimum detector performance for signals in cyclostationary noise, IEEE J. Ocean. Eng. 12 (1987) 97–115.
  * [4] P. Bloomfield, H. L. Hurd, R. B. Lund, Periodic correlation in stratospheric ozone time series, J. Time Ser. Anal. 15 (2) (1994) 127–150.
  * [5] E. Broszkiewicz-Suwaj, A. Makagon, R. Weron, A. Wyłomańska, On detecting and modeling periodic correlation in financial data, Physica A 336 (1-2) (2004) 196–205.
  * [6] J. Antoni, J. Daniere, F. Guillet, Effective vibration analysis of IC engines using cyclostationarity. Part I - A methodology for condition monitoring, J. Sound Vib. 257 (5) (2002) 815–837.
  * [7] J. Wodecki, A. Michalak, R. Zimroz, T. Barszcz, A. Wyłomańska, Impulsive source separation using combination of Nonnegative Matrix Factorization of bi-frequency map, spatial denoising and Monte Carlo simulation, Mech. Syst. Signal Process. 127 (2019) 89–101.
  * [8] R.-B. Sun, Z.-B. Yang, K. Gryllias, X.-F. Chen, Cyclostationary modeling for local fault diagnosis of planetary gear vibration signals, J. Sound Vib. 471 (2020) 115175.
  * [9] J. Wodecki, A. Michalak, A. Wyłomańska, R. Zimroz, Influence of non-Gaussian noise on the effectiveness of cyclostationary analysis – simulations and real data analysis, Measurement 171 (2021) 108814.
  * [10] M. Esfandiari, S. A. Vorobyov, M. Karimi, New estimation methods for autoregressive process in the presence of white observation noise, Signal Process. 171 (2020) 107480.
  * [11] R. Diversi, R. Guidorzi, U. Soverini, Identification of autoregressive models in the presence of additive noise, Int. J. Adapt. Control Signal Process. 22 (5) (2008) 465–481.
  * [12] A. Mahmoudi, M. Karimi, Parameter estimation of autoregressive signals from observations corrupted with colored noise, Signal Process. 90 (1) (2010) 157–164.
  * [13] O. Çayır, Ç. Candan, Maximum likelihood autoregressive model parameter estimation with noise corrupted independent snapshots, Signal Process. 186 (2021) 108118.
  * [14] A. Mahmoudi, M. Karimi, H. Amindavar, Parameter estimation of autoregressive signals in presence of colored AR(1) noise as a quadratic eigenvalue problem, Signal Process. 92 (4) (2012) 1151–1156.
  * [15] A. Mahmoudi, M. Karimi, Inverse filtering based method for estimation of noisy autoregressive signals, Signal Process. 91 (7) (2011) 1659–1664.
  * [16] M. Hasan, M. Hossain, M. Haque, Parameter estimation of multichannel autoregressive processes in noise, Signal Process. 83 (3) (2003) 603–610.
  * [17] A. Mahmoudi, M. Karimi, Estimation of the parameters of multichannel autoregressive signals from noisy observations, Signal Process. 88 (11) (2008) 2777–2783.
  * [18] R. Diversi, R. Guidorzi, Fast filtering of noisy autoregressive signals, Signal Process. 87 (11) (2007) 2843–2849.
  * [19] D. Labarre, E. Grivel, Y. Berthoumieu, E. Todini, M. Najim, Consistent estimation of autoregressive parameters from noisy observations based on two interacting Kalman filters, Signal Process. 86 (10) (2006) 2863–2876, Special Section: Fractional Calculus Applications in Signals and Systems.
  * [20] A. J. Q. Sarnaglia, V. A. Reisen, C. Lévy-Leduc, Robust estimation of periodic autoregressive processes in the presence of additive outliers, J. Multivar. Anal. 101 (9) (2010) 2168–2183.
  * [21] A. J. Q. Sarnaglia, V. A. Reisen, P. Bondou, C. Lévy-Leduc, A robust estimation approach for fitting a PARMA model to real data, in: 2016 IEEE Statistical Signal Processing Workshop (SSP), 2016, pp. 1–5.
  * [22] A. A. Samadi, A. M. Al-Quraam, Estimation of the seasonal ACF of PAR(1) model in the presence of additive outliers, J. Appl. Stat. Sci. 19 (2) (2011) 169–182.
  * [23] A. J. Q. Sarnaglia, V. A. Reisen, P. Bondon, C. Lévy-Leduc, M-regression spectral estimator for periodic ARMA models. An empirical investigation, Stoch. Environ. Res. Risk Assess. 35 (3) (2021) 653–664.
  * [24] Q. Shao, Robust estimation for periodic autoregressive time series, J. Time Ser. Anal. 29 (2) (2008) 251–263.
  * [25] W. Żuławiński, A. Wyłomańska, New estimation method for periodic autoregressive time series of order 1 with additive noise, Int. J. Adv. Eng. Sci. Appl. Math. 13 (2021) 163–176.
  * [26] K. P. Burnham, D. R. Anderson, Model Selection and Multimodel Inference: A practical information-theoretic approach (2nd ed.), Springer-Verlag, 2002.
  * [27] D. R. Anderson, Model Based Inference in the Life Sciences, Springer, 2008.
  * [28] M. Jabłońska-Sabuka, M. Teuerle, A. Wyłomańska, Bivariate sub-Gaussian model for stock index returns, Physica A 486 (2017) 628–637.
  * [29] A. Vecchia, Periodic autoregressive-moving average (PARMA) modeling with applications to water resources, J. Am. Water Resour. Assoc. 21 (5) (1985) 721–730.
  * [30] D. Najarzadeh, Testing independence in high-dimensional multivariate normal data, Commun. Stat. Theory Methods 50 (14) (2021) 3421–3435.
  * [31] D. He, H. Liu, K. Xu, M. Cao, Generalized Schott type tests for complete independence in high dimensions, J. Multivar. Anal. 183 (2021) 104731.
  * [32] S. Chang, Y. Qi, On Schott’s and Mao’s test statistics for independence of normal random vectors, Stat. Probab. Lett. 140 (2018) 132–141.
  * [33] B. Ebner, N. Henze, Tests for multivariate normality—a critical review with emphasis on weighted $l^{2}$-statistics, TEST 29 (2020) 845–892.
  * [34] M. Arnold, D. Iskander, A. Zoubir, Testing Gaussianity with the characteristic function, in: 1995 International Conference on Acoustics, Speech, and Signal Processing, Vol. 3, 1995, pp. 2012–2015 vol.3.
  * [35] G. W. Inverarity, Fast computation of multidimensional Fourier integrals, SIAM J. Sci. Comput. 24 (2) (2002) 645–651.
  * [36] V. Witkovský, Numerical inversion of a characteristic function: An alternative tool to form the probability distribution of output quantity in linear measurement models, ACTA IMEKO 5 (3) (2016) 32–44.

##  Appendix A Matrix representation of residuals for noise-corrupted PAR($p$)
model

In this part, the matrix representation of the residual series defined in Eq.
(6) is provided. Its form is as follows:

$\mathbf{R}_{n}=\bm{\xi}_{n}+\mathbf{Z}_{n}\mathbf{A},$ (20)

where $\bm{\xi}_{n}=[\xi_{nT+1},\ldots,\xi_{nT+T}]$ is a $1\times T$ vector of
i.i.d. zero mean components with variance $\sigma_{\xi}^{2}$,
$\mathbf{Z}_{n}=[Z_{nT+T},\ldots,Z_{nT+1-p}]$ is a $1\times(p+T)$ vector of
i.i.d. zero mean components with variance $\sigma_{Z}^{2}$, independent of
$\bm{\xi}_{n}$, and $\mathbf{A}$ is a $(p+T)\times T$ matrix of elements
$a_{kl}$ defined as

$a_{kl}=-\phi_{k+l-T-1}(l),\quad k=1,\cdots,p+T,\,\,\,l=1,\cdots,T$ (21)

with $\phi_{0}(\cdot)=-1$ and $\phi_{j}(\cdot)=0$ for
$j\notin\\{0,\cdots,p\\}$. Let us note that the representation given in Eq.
(20) immediately produces the covariance matrix of $\mathbf{R}_{n}$ to be on
the form

$\mathbf{\Gamma}^{R}=\sigma_{\xi}^{2}\mathbf{I}_{T}+\sigma_{Z}^{2}\mathbf{A^{\prime}A},$
(22)

where $\mathbf{I}_{T}$ is a $T\times T$ identity matrix, leading directly to
the formulas in Eq. (10).

##  Appendix B Mixture of Gaussian distributions

Let $a_{1},a_{2},\cdots,a_{m}$ denote a series of non-negative weights
satisfying $\sum_{i=1}^{m}a_{i}=1$. Let
$F_{1}(\cdot),F_{2}(\cdot),\cdots,F_{m}(\cdot)$ denote an arbitrary sequence
of Gaussian cumulative distribution functions (CDFs) and let
$f_{1}(\cdot),f_{2}(\cdot),\cdots,f_{m}(\cdot)$ be the PDFs with means
$\mu_{1},\mu_{2},\cdots,\mu_{m}$ and variances
$\sigma_{1}^{2},\sigma_{2}^{2},\cdots,\sigma_{m}^{2}$. A random variable $X$
with the following CDF and PDF:

$\displaystyle F(z)=\sum_{i=1}^{m}a_{i}F_{i}(z),\leavevmode\nobreak\
\leavevmode\nobreak\ f(z)=\sum_{i=1}^{m}a_{i}f_{i}(z)$ (23)

is called a mixture of Gaussian distributions. In this paper for the
distribution of the additive noise we consider two types of mixture of
Gaussian distributions, namely with $m=1$ and $m=2$. In both cases, the
considered distributions are centered. For $m=1$ this distribution reduces to
the Gaussian distribution. In that case the excess kurtosis is equal to zero.
For $m=2$ and $\mu_{1}=\mu_{2}=0$ the excess kurtosis is equal to
$Kur=\frac{3(a_{1}\sigma_{1}^{4}+a_{2}\sigma_{2}^{4})}{(a_{1}\sigma_{1}^{2}+a_{2}\sigma_{2}^{2})^{2}}-3$.

##  Appendix C Distribution of the residuals for PAR($p$) model disturbed by
noise with mixture of two Gaussian distributions

To show that the distribution of the residual vector in Example 2 is a mixture
of eight zero-mean bivariate Gaussian distributions, one can use the matrix
representation given in A. Let us notice that in that case using Eq. (20) we
have

$\mathbf{R}_{2}\overset{d}{=}\bm{\xi}+\mathbf{ZA},$ (24)

where $\bm{\xi}=[\xi_{1},\xi_{2}]$ is a $1\times 2$ vector of i.i.d. zero mean
Gaussian components with variance $\sigma_{\xi}^{2}$,
$\mathbf{Z}=[Z_{1},Z_{2},Z_{3}]$ is a $1\times 3$ vector of i.i.d. zero mean
components with variance $\sigma_{Z}^{2}$, independent of $\bm{\xi}$, and
$\mathbf{A}$ is a $3\times 2$ matrix of constant coefficients of the form

$\mathbf{A}=\begin{pmatrix}0&1\\\ 1&-\phi_{1}(2)\\\
-\phi_{1}(1)&0\end{pmatrix}.$ (25)

Now, since the $\\{Z_{i}\\}$ are two-component mixtures of Gaussian
distributions, they can be represented as

$Z_{i}\overset{d}{=}I_{i}X_{1}^{(i)}+(1-I_{i})X_{2}^{(i)},\quad i=1,2,3,$ (26)

where the $\\{X_{1}^{(i)}\\}$ are i.i.d. zero-mean Gaussian variables with
variance $\sigma_{1}^{2}$, the $\\{X_{2}^{(i)}\\}$ are i.i.d. zero-mean
Gaussian variables with variance $\sigma_{2}^{2}$, and the $\\{I_{i}\\}$ are
i.i.d. Bernoulli variables with $\mathbb{P}(I_{i}=1)=a_{1}$ and
$\mathbb{P}(I_{i}=0)=a_{2}=1-a_{1}$, with all the variables in Eq. (26)
mutually independent. One can note that there are eight different combinations
of the values of $I_{i}$, and conditionally on these values, the vector
$\mathbf{Z}=[Z_{1},Z_{2},Z_{3}]$ in Eq. (24) is zero-mean multivariate
Gaussian, with the diagonal covariance matrix
$\bm{\Omega}=\text{diag}(\omega_{1},\omega_{2},\omega_{3})$ where

$\omega_{i}=I_{i}\sigma_{1}^{2}+(1-I_{i})\sigma_{2}^{2},\quad i=1,2,3.$ (27)

Thus, conditionally on the values of $I_{i}$, the random vector
$\mathbf{R}_{2}$ in Eq. (24) is also multivariate Gaussian, with the following
covariance matrix:

$\sigma_{\xi}^{2}\mathbf{I}_{2}+\mathbf{A}^{\prime}\bm{\Omega}\mathbf{A}.$
(28)

This shows that the distribution of $\mathbf{R}_{2}$ is a mixture of eight
zero-mean bivariate Gaussian distributions, with mixing probabilities given by
$\pi_{ijk}=\mathbb{P}(I_{1}=i,\,I_{2}=j,\,I_{3}=k),\,i,\,j,\,k=1,2$. One can
notice that these probabilities reduce to
$a_{1}^{3},\,a_{1}^{2}a_{2},\,a_{1}a_{2}^{2}$, and $a_{2}^{3}$, used in Eq.
(16) provided in Example 2. It is also worth mentioning that the above
findings could be written in a more general case, assuming that the additional
noise is an $m$-component mixture of zero-mean Gaussian distributions with
variances $\sigma_{1}^{2},\cdots,\sigma_{m}^{2}$. However, as the above
example is used only for the purpose of illustrating the general methodology,
we choose to omit that part of the presentation.

##  Appendix D Modified errors-in-variables method of estimation for PAR model
with additive noise

In this part we describe the estimation method used in this paper. It can be
considered as a generalization of the algorithm proposed in [11] for
autoregressive model with additive noise, utilizing both low- and high-order
Yule-Walker equations based on the periodic autocovariance function
${\gamma}^{Y}(w,k)=\mathbb{E}Y_{nT+w}Y_{nT+w-k}$ (of a zero-mean process). For
the process $\\{Y_{t}\\}$, for $v=1,\cdots,T$, we define low-order Yule-Walker
equations as

$\displaystyle(\mathbf{\Gamma}^{Y}_{v}-\sigma^{2}_{Z}\mathbf{I}_{p})\mathbf{\Phi}_{v}=\mathbf{\gamma}^{Y}_{v},$
(29)

where

$\displaystyle(\mathbf{\Gamma}^{Y}_{v})_{i,j}=\gamma^{Y}(v-i,j-i),\leavevmode\nobreak\
\mathbf{\gamma}^{Y}_{v}=[\gamma^{Y}(v,1),\cdots,\gamma^{Y}(v,p)]^{\prime},\leavevmode\nobreak\
\leavevmode\nobreak\
\mathbf{\Phi}_{v}=[\phi_{1}(v),\cdots,\phi_{p}(v)]^{\prime}$ (30)

and $i,j=1,2,\cdots,p$. On the other hand, we define the system of $s$ high-
order Yule-Walker equations for the process $\\{Y_{t}\\}$, again for
$v=1,\cdots,T$, in the following way:

${}_{s}\tilde{\mathbf{\Gamma}}^{Y}_{v}\mathbf{\Phi}_{v}={}_{s}\tilde{\mathbf{\gamma}}^{Y}_{v}.$
(31)

where ${}_{s}\tilde{\mathbf{\Gamma}}^{Y}_{v}$ is a $s\times p$ matrix of
elements

$\displaystyle({}_{s}\tilde{\mathbf{\Gamma}}^{Y}_{v})_{i,j}=\gamma^{Y}(v-j,p+i-j)\leavevmode\nobreak\
\text{ while}\leavevmode\nobreak\
_{s}\tilde{\mathbf{\gamma}}^{Y}_{v}=[\gamma^{Y}(v,p+1),\cdots,\gamma^{Y}(v,p+s)]^{\prime}.$
(32)

In the estimation, the ${\gamma}^{Y}(w,k)$ expression is replaced by its
empirical version defined for a sample $y_{1},\cdots,y_{NT}$ as

$\displaystyle\hat{\gamma}^{Y}(w,k)=\frac{1}{N}\sum_{n=l}^{r}y_{nT+w}y_{nT+w-k},$
(33)

where

$\displaystyle
l=\max\left(\left\lceil\frac{1-w}{T}\right\rceil,\left\lceil\frac{1-(w-k)}{T}\right\rceil\right),\quad
r=\min\left(\left\lfloor\frac{NT-w}{T}\right\rfloor,\left\lfloor\frac{NT-(w-k)}{T}\right\rfloor\right).$
(34)

The main idea of the presented method is finding an estimate of additive noise
variance using high-order Yule-Walker equations and then utilizing it in low-
order Yule-Walker equations for estimation of model’s coefficients and
innovations’ variance. Let us define the functions
$\mathbf{\Phi}^{*}_{v}(\cdot)$, $\sigma^{2*}_{\xi}(v)(\cdot)$ of the parameter
$\sigma_{Z}^{2*}$ based on low-order Yule-Walker equations for $\\{Y_{t}\\}$:

$\displaystyle\hat{\mathbf{\Phi}}^{*}_{v}(\sigma^{2*}_{Z})=(\hat{\mathbf{\Gamma}}^{Y}_{v}-\sigma^{2*}_{Z}\mathbf{I}_{p})^{-1}\hat{\mathbf{\gamma}}^{Y}_{v},\leavevmode\nobreak\
\leavevmode\nobreak\
\hat{\sigma}^{2*}_{\xi}(v)(\sigma^{2*}_{Z})=\hat{\gamma}^{Y}(v,0)-\hat{\mathbf{\Phi}}^{*}_{v}(\sigma^{2*}_{Z})^{\prime}\hat{\mathbf{\gamma}}^{Y}_{v}-\sigma^{2*}_{Z},$
(35)

where $\hat{\mathbf{\Gamma}}^{Y}_{v}$ and $\hat{\mathbf{\gamma}}^{Y}_{v}$ are
constructed using Eq. (30) by replacing all the terms with their empirical
counterparts,

$\displaystyle(\hat{\mathbf{\Gamma}}^{Y}_{v})_{i,j}=\hat{\gamma}^{Y}(v-i,j-i),\leavevmode\nobreak\
\leavevmode\nobreak\
\hat{\mathbf{\gamma}}^{Y}_{v}=[\hat{\gamma}^{Y}(v,1),\cdots,\hat{\gamma}^{Y}(v,p)]^{\prime},$
(36)

for $i,j=1,2,\cdots,p$. To construct an upper bound for $\sigma^{2*}_{Z}$,
similarly as in [11], let us also consider the following matrices
$\hat{G}^{Y}_{v}$ for each $v=1,\cdots,T$:

$\displaystyle\hat{\mathbf{G}}_{v}^{Y}=\begin{bmatrix}\hat{\gamma}^{Y}(v,0)&\hat{\mathbf{\gamma}}_{v}^{Y\prime}\\\
\hat{\mathbf{\gamma}}^{Y}_{v}&\hat{\mathbf{\Gamma}}^{Y}_{v}\end{bmatrix}.$
(37)

As the estimated additive noise variance $\hat{\sigma}^{2*}_{Z}$ we consider
the value which minimizes the following high-order Yule-Walker equations-based
cost function

$\displaystyle
J_{\text{total}}(\sigma^{2*}_{Z})=\sum_{v=1}^{T}J_{v}(\sigma^{2*}_{Z}),$ (38)

over an interval $\sigma^{2*}_{Z}\in[0,\zeta]$, for
$\zeta=\min\\{\min\text{eig}(\hat{\mathbf{G}}^{Y}_{v}),v=1,\cdots,T\\}$, where

$\displaystyle
J_{v}({\sigma}^{2*}_{Z})=||{}_{s}\hat{\tilde{\mathbf{\Gamma}}}^{Y}_{v}\hat{\mathbf{\Phi}}^{*}_{v}({\sigma}^{2*}_{Z})-{}_{s}\hat{\tilde{\mathbf{\gamma}}}^{Y}_{v}||^{2}_{2},$
(39)
$\displaystyle({}_{s}\hat{\tilde{\mathbf{\Gamma}}}^{Y}_{v})_{i,j}=\hat{\gamma}^{Y}(v-j,p+i-j),\leavevmode\nobreak\
\leavevmode\nobreak\
_{s}\hat{\tilde{\mathbf{\gamma}}}^{Y}_{v}=[\hat{\gamma}^{Y}(v,p+1),\cdots,\hat{\gamma}^{Y}(v,p+s)]^{\prime}$
(40)

for $i=1,2,\cdots,s$, $j=1,2,\cdots,p$. At the end, having the optimal
$\hat{\sigma}^{2*}_{Z}$, for each $v=1,\cdots,T$ we calculate

$\displaystyle\hat{\mathbf{\Phi}}_{v}=(\hat{\mathbf{\Gamma}}^{Y}_{v}-\hat{\sigma}^{2}_{Z}\mathbf{I}_{p})^{-1}\hat{\mathbf{\gamma}}^{Y}_{v},\leavevmode\nobreak\
\hat{\sigma}^{2}_{\xi}(v)=\hat{\gamma}^{Y}(v,0)-\hat{\mathbf{\Phi}}_{v}^{\prime}\hat{\mathbf{\gamma}}^{Y}_{v}-\hat{\sigma}^{2}_{Z},$
(41)

and estimate the variance of the innovations $\hat{\sigma}^{2}_{\xi}$, which
is the mean of all the $\hat{\sigma}^{2}_{\xi}(v)$ across $v=1,\cdots,T$. The
entire procedure is presented in Algorithm 1.

Algorithm 1 Modified errors-in-variables method

1:Set value of $s$ (where $s\geq p$). In this paper, we set $s=p$.

2:For each $v=1,\cdots,T$:

  1. i:

Construct $\hat{\mathbf{\Gamma}}^{Y}_{v}$ (Eq. (36)),
$\hat{\mathbf{\gamma}}^{Y}_{v}$ (Eq. (36)),
${}_{s}\hat{\tilde{\mathbf{\Gamma}}}^{Y}_{v}$ (Eq. (40)) and
${}_{s}\hat{\tilde{\mathbf{\gamma}}}^{Y}_{v}$ (Eq. (40)).

  2. ii:

Construct $\hat{\mathbf{G}}^{Y}_{v}$ (Eq. (37)) and compute
$\min\text{eig}(\hat{\mathbf{G}}^{Y}_{v})$.

3:Compute
$\zeta=\min\left\\{\min\text{eig}(\hat{\mathbf{G}}^{Y}_{v}),\>:\>v=1,\cdots,T\right\\}$

4:Determine $\hat{\sigma}^{2}_{Z}$ – a value which minimizes
$J_{\text{total}}(\sigma^{2*}_{Z})$ (Eq. (38)) over interval
$\sigma^{2*}_{Z}\in[0,\zeta]$.

5:For each $v=1,\cdots,T$ :

  1. i:

Compute
$\hat{\mathbf{\Phi}}_{v}=(\hat{\mathbf{\Gamma}}^{Y}_{v}-\hat{\sigma}^{2}_{Z}\mathbf{I}_{p})^{-1}\hat{\mathbf{\gamma}}^{Y}_{v}$.

  2. ii:

Compute
$\hat{\sigma}^{2}_{\xi}(v)=\hat{\gamma}^{Y}(v,0)-\hat{\mathbf{\Phi}}_{v}^{\prime}\hat{\mathbf{\gamma}}^{Y}_{v}-\hat{\sigma}^{2}_{Z}$.

6:Compute $\hat{\sigma}^{2}_{\xi}=1/T\sum_{v=1}^{T}\hat{\sigma}^{2}_{\xi}(v)$