# Quickest Change Detection in Autoregressive Models

Zhongchang Sun Shaofeng Zou Zhongchang Sun and Shaofeng Zou are with the
Department of Electrical Engineering, University at Buffalo, Buffalo, NY 14228
USA (e-mail<EMAIL_ADDRESS>szou3@buffalo.edu).

###### Abstract

The problem of quickest change detection (QCD) in autoregressive (AR) models
is investigated. A system is being monitored with sequentially observed
samples. At some unknown time, a disturbance signal occurs and changes the
distribution of the observations. The disturbance signal follows an AR model,
which is dependent over time. Before the change, observations only consist of
measurement noise, and are independent and identically distributed (i.i.d.).
After the change, observations consist of the disturbance signal and the
measurement noise, are dependent over time, which essentially follow a
continuous-state hidden Markov model (HMM). The goal is to design a stopping
time to detect the disturbance signal as quickly as possible subject to false
alarm constraints. Existing approaches for general non-i.i.d. settings and
discrete-state HMMs cannot be applied due to their high computational
complexity and memory consumption, and they usually assume some asymptotic
stability condition. In this paper, the asymptotic stability condition is
firstly theoretically proved for the AR model by a novel design of forward
variable and auxiliary Markov chain. A computationally efficient Ergodic CuSum
algorithm that can be updated recursively is then constructed and is further
shown to be asymptotically optimal. The data-driven setting where the
disturbance signal parameters are unknown is further investigated, and an
online and computationally efficient gradient ascent CuSum algorithm is
designed. The algorithm is constructed by iteratively updating the estimate of
the unknown parameters based on the maximum likelihood principle and the
gradient ascent approach. The lower bound on its average running length to
false alarm is also derived for practical false alarm control. Simulation
results are provided to demonstrate the performance of the proposed
algorithms.

###### Index Terms:

Hidden Markov model, forward variable, sequential change detection, asymptotic
optimality, non-i.i.d..

##  I Introduction

The problem of quickest change detection (QCD) has been widely studied in the
literature [1, 2, 3, 4], where the goal is to detect an abrupt change in the
data-generating distribution as quickly as possible, subject to false alarm
constraints. For the basic setting where the observations are independent over
time, algorithms, e.g., cumulative sum (CuSum) algorithm and Shiryaev-Roberts
algorithm have been proposed and have been shown to be optimal [5, 6, 7, 8,
9]. However, for a wide range of practical applications, observations may not
be independent over time. For example, in power systems, the faulty signal is
dependent in time [10] and such dependency is usually exploited by an
autoregressive (AR) model [11, 12, 13, 14, 15].

In this paper, we investigate the problem of QCD in AR models. Specifically,
before the change, the observed signal only consists of Gaussian measurement
noise. After the change, a non-i.i.d. disturbance signal occurs in the system,
and the observation consists of the disturbance signal and the measurement
noise. The disturbance signal follows an AR model, and thus the observation
equivalently follows a hidden Markov model (HMM) [16]. The goal is to detect
the occurrence of the disturbance signal as quickly as possible subject to
false alarm constraints.

###  I-A Related Works

The problem in this paper is closely related to QCD in Markov chains [17, 18,
19, 20, 21, 22, 23, 24, 25], where the pre- and/or post-change samples follow
the Markov model. In [17], optimal change detection algorithms were proposed
for finite-state Markov chains under the Bayesian setting and the minimax
setting. The continuous state Markov chain was studied in [24], where the
CuSum algorithm was proved to be asymptotically optimal and its performance
was characterized. In [19], the change detection in sensor networks under the
Bayesian setting was studied where the change propagates across sensors and
its propagation was modeled as a Markov process. It was shown that an exactly
optimal algorithm for QCD of Markov chains under the Bayesian setting is a
threshold test on the posterior of no change has happened in [23]. In [21],
two CuSum type schemes were proposed for the QCD of finite-state Markov chain
with unknown post-change transition probabilities, and bounds on the average
detection delay (ADD) and the average running length (ARL) were further given
for the proposed schemes. A maximum mean discrepancy (MMD) based method for
change detection of Markov kernels with unknown post-change kernel was
proposed in [22, 25] for general state Markov chain. However, the algorithms
and analyses for Markov chains cannot be applied to our problem directly. This
is because with measurement noise, the post-change samples in this paper
essentially follow an HMM, and the hidden state is not directly observable.

The problem of QCD in discrete-state HMMs has been studied in [26, 27, 28, 29,
30, 31]. In [26], Page’s test was extended to the QCD of discrete-state HMMs
and the ARL under both hypotheses were approximated. In [28], the detection of
the change point at which the discrete-state HMM enters the absorbing state
and the identification of the absorbing set was investigated and optimal
decision rule was provided. The QCD problem in two-state HMM was studied in
[29], where two computationally efficient schemes were developed. In [30], the
Shiryaev algorithm was proved to be asymptotically optimal under some
regularity conditions for the QCD in discrete-state HMMs under the Bayesian
setting. The exact optimal algorithm for QCD in discrete-state HMM under the
Bayesian setting was established in [31]. These works mainly focus on the
discrete-state HMM. The AR model in our paper is an HMM model with a
continuous state space. To the best of the author’s knowledge, research on QCD
in continuous HMMs is rather limited. The performance characterization of
Shiryaev algorithm for the discrete-state HMM relies on the limiting theory
for products of random matrices [32] and the computational complexity scales
with the size of the state space. For the continuous-state Markov chain, the
theory of random matrices is not applicable anymore since the number of state
is uncountably infinite.

The AR models are commonly used in power system to characterize the
disturbance signal and inter-area oscillations [11, 12, 13, 14, 15]. In [33,
34, 35, 12, 36, 37, 38], the change detection in the AR model was studied
where the observation is an AR time series. Our problem is fundamentally
different from [33, 34, 35, 12, 36, 37, 38] since the disturbance signal which
has AR structure is unobservable in our problem. The observation is a noised
version of the AR time series which makes our problem more challenging. The
quickest detection of cyber-attacks in discrete-time linear dynamic system,
which is an AR model, was studied in [14]. The Kalman filter [39] was first
performed to estimate the state and the generalized CuSum algorithm was
further designed using the state estimation. Though the Kalman filter is
effective in state estimation, the performance of the detection algorithm is
difficult to characterize due to the estimation error. In [40], a robust
algorithm was proposed for distributed change detection in AR models with
noised observations. However, no theoretical guarantee was provided. In our
paper, we design an Ergodic CuSum algorithm directly and show that it is
asymptotically optimal under the Lorden’s criterion [41].

In our paper, we consider both the model-based setting where parameters of the
disturbance signal are known and the data-driven setting where parameters of
the disturbance signal are unknown. For the data-driven setting, our problem
can be viewed as a composite QCD problem [42, 24, 43, 44, 45, 46, 47, 48].
Existing works mainly assume that the samples are i.i.d., while in the AR
models, the observations are dependent in time. In [15], the QCD in AR model
was studied where the properties of the post-change signal are unknown. Based
on the assumption that the parameters of the post-change signal were small,
the generalized likelihood ratio test (GLRT) was proposed. Our work is
different from the works in [15] since we do not assume that the parameter of
the post-change signal is small. Moreover, simulation results demonstrate that
our proposed algorithm outperforms the GLRT in [15].

###  I-B Contributions and Major Challenges

For the model-based setting, we design an Ergodic CuSum algorithm which can be
updated recursively and thus is computationally efficient, and we
theoretically prove that it is asymptotically optimal. While the general
theory for non-i.i.d. setting were developed in e.g., [24, 49, 50, 51], it
relies on some asymptotic stability condition that the normalized log
likelihood ratio between the post- and pre-change distributions converges to
some finite and positive number. However, there are very few studies that try
to verify this condition, e.g., discrete HMMs in [30] and Markov model in
[49]. For the AR model, whether such stability condition holds remains
unexplored. In this paper, we show that the normalized log likelihood ratio
converges to some $\mathcal{K}>0$ almost surely under the post-change
distribution, i.e., satisfies the stability condition.

The difficulty of analyzing the convergence of the log likelihood ratio lies
in that the log likelihood function for the AR model has a non-additive form
due to the non-i.i.d. nature of the observation. For the discrete-state HMM in
[30], the likelihood ratio function is represented by the $L_{1}$-norm of
products of Markov random matrices and thus the log likelihood ratio has an
additive form. However, the Markov random matrices techniques can not be
applied to our problem since our AR model has uncountably infinite hidden
states. To overcome this difficulty, we represent the likelihood function as
the integration of the product of a sequence of functions based on the hidden
Markov structure of the observation. We further show that the likelihood
function can be represented as the integration of a scaled Gaussian density
function. We define this scaled Gaussian density function as the forward
variable and show that it can be updated recursively. Since the scaled
Gaussian density function has parameterized representation, we develop a novel
approach to write the conditional likelihood ratio using the parameter of the
forward variable and the observation. The log likelihood ratio can then be
written as the sum of the conditional log likelihood ratio and thus has an
additive form. To characterize the convergence of the log likelihood ratio, we
design an auxiliary Markov chain using the parameter of the forward variable
and the observation and represent the likelihood function as a function of the
auxiliary Markov chain. We further show that the stationary distribution of
this auxiliary Markov chain exists and then show that $\mathcal{K}$ exists by
applying the ergodic theorem [52, 53] to the auxiliary Markov chain. With the
convergence of the log likelihood ratio, we derive the universal lower bound
on the worst-case average detection delay (WADD).

For the CuSum algorithm designed for the general non-i.i.d. setting, e.g.,
[24], it is computationally expensive for our AR model since the likelihood
ratio for our AR model depends on the change point and we need to compute the
likelihood ratio for every possible change point. Motivated by the dependence
of the WADD lower bound on ${\mathcal{K}}$, we aim to find a computationally
efficient test statistic that has a positive drift of $\mathcal{K}$ under the
post-change distribution and has a negative drift under the pre-change
distribution. We propose a novel Ergodic CuSum algorithm by exploiting the
ergodicity of the underlying HMM and employing the likelihood ratio of the
observations with change point at time $1$. Since the likelihood function can
be written as a function of the observation and the forward variable, our
Ergodic CuSum algorithm can then be updated recursively, and thus is
computationally efficient. From the ergodic theorem [52, 53] of Markov chain,
the long term statistical performance of our test statistic is independent of
the state when the change occurs. Based on this fact, the test statistic of
our Ergodic CuSum algorithm converges to $\mathcal{K}$ in the post-change
regime despite of the fact that time 1 is not the true change point. Under the
pre-change distribution, we show that our test statistic has a negative drift.
Moreover, since our test statistic is a likelihood ratio function of the
current observation, the ARL lower bound can be derived following the proof of
the ARL lower bound for general non-i.i.d. case. Therefore, we show that when
satisfying the constraint on the average running length (ARL), the WADD upper
bound of our Ergodic CuSum algorithm matches with the universal lower bound
and thus the asymptotic optimality of our Ergodic CuSum algorithm follows.

For the data-driven setting, the generalized likelihood ratio test (GLRT)
which replaces the unknown post-change parameter with its maximum likelihood
estimate (MLE) [24] is not computationally efficient. We propose an online
gradient ascent CuSum algorithm (OGA-CuSum) based on the online convex
optimization (OCO) algorithm [54, 48]. We iteratively update the estimate of
the unknown parameters based on the maximum likelihood principle using
gradient ascent. Specifically, at each time, we compute the gradient of the
log likelihood ratio function with respect to the unknown parameters based on
the current observation and update the estimate in the direction of the
gradient. We then replace the unknown parameter in the Ergodic CuSum with its
estimate to construct the OGA-CuSum. The OGA-CuSum can be updated recursively
using only the most recent sample and thus is memory and computationally
efficient. We derive a lower bound on its ARL so that a threshold can be
chosen analytically to control the false alarm. We also provide simulation
results to demonstrate the performance of our algorithms.

###  I-C Paper Organization

The remainder of this paper is organized as follows. In Section II, we present
the problem formulation. In Section III, we present the design of a forward
variable and establish the universal lower bound on WADD. In Section IV, we
develop a computationally efficient Ergodic CuSum algorithm and show that it
is asymptotically optimal under Lorden’s criterion. In Section V, we design an
online gradient ascent CuSum algorithm for the case where the post-change
parameters are unknown. Numerical results are provided to demonstrate the
performance of our algorithms in VI. In Section VII, we present some
concluding remarks.

##  II Problem formulation

Consider a system being monitored by a sequentially observed signal
$\bm{y}_{t}\in\mathbb{R}^{K},t=1,2,\cdots$. At some unknown time $t_{0}$, a
disturbance signal occurs in the system and changes the distribution of the
observed signal. Specifically, before the change point $t_{0}$, the observed
signal consists of the measurement noise $\bm{\nu}_{t}$ only, which is
Gaussian distributed and independent over time:

$\displaystyle\bm{y}_{t}=\bm{\nu}_{t}\sim\mathcal{N}(\bm{0},\bm{I}),\
t<t_{0},$ (1)

where $\bm{I}$ is the $K\times K$ identity matrix. The post-change signal
consists of the disturbance signal $\bm{x}_{t}$ and the measurement noise:

$\displaystyle\bm{y}_{t}=\bm{x}_{t}+\bm{\nu}_{t},\ t\geq t_{0},$ (2)

Note that the identity covariance matrix of $\bm{\nu}_{t}$ can be satisfied by
applying a linear transformation to whiten the noise in the observation
$\bm{y}_{t}$. The disturbance signal follows an AR model [11, 12, 13, 14, 15]:

$\displaystyle\bm{x}_{t}=\sum_{i=1}^{q}\bm{A}_{i}\bm{x}_{t-i}+\bm{\omega}_{t},$
(3)

where $\bm{A}_{i}\in\mathbb{R}^{K\times K}$ is the matrix coefficient and is
invertible, $\bm{\omega}_{t}\in\mathbb{R}^{K}$ is the innovation noise vector
and $\bm{\omega}_{t}\sim\mathcal{N}(\bm{0},\bm{R}_{\omega})$ and $q$ denotes
the order of the AR model. The disturbance signal $\bm{x}_{t}$ is not directly
observable. The goal is to detect the change at time $t_{0}$ as soon as
possible subject to false alarm constraints.

Denote by $\mathbb{P}_{\infty}$ the probability measure when there is no
change, and denote by $p_{\infty}$ the corresponding probability density. For
any $t_{0}>0$, denote by $\mathbb{P}_{t_{0}}$ the probability measure when the
change happens at $t_{0}$, and denote by $p_{t_{0}}$ the corresponding
probability density. It is clear that

$\displaystyle p_{\infty}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})$
$\displaystyle=p_{\infty}(\bm{y}_{t}),$ $\displaystyle
p_{t_{0}}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})$
$\displaystyle=p_{\infty}(\bm{y}_{t}),\text{ if }t<t_{0},$ $\displaystyle
p_{t_{0}}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})$
$\displaystyle=p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1}),\text{
if }t\geq t_{0}.$ (4)

The AR model is equivalent to an HMM. Let
$f_{1}(\bm{x}_{t_{0}},\bm{x}_{t_{0}+1},\cdots,\bm{x}_{t_{0}+q-1})$ be the
joint probability density of the initial state
$(\bm{x}_{t_{0}},\bm{x}_{t_{0}+1},\cdots,\bm{x}_{t_{0}+q-1})$. Let
$f(\bm{x}_{t}|\bm{x}_{t-q},\cdots,\bm{x}_{t-1})$ denote the transitional
probability density of the HMM and $g(\bm{y}_{t}|\bm{x}_{t})$ be the
conditional probability density of the observation $\bm{y}_{t}$ given the
hidden state $\bm{x}_{t}$. The probability density $p_{t_{0}}$ then satisfies
that for $t\geq t_{0}+q$

$\displaystyle p_{t_{0}}(\bm{y}_{1},\cdots,\bm{y}_{t})=$ $\displaystyle
p_{\infty}(\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})\cdot\int
f_{1}(\bm{x}_{t_{0}},\bm{x}_{t_{0}+1},\cdots,\bm{x}_{t_{0}+q-1})g(\bm{y}_{t_{0}}|\bm{x}_{t_{0}})\cdots
g(\bm{y}_{t_{0}+q-1}|\bm{x}_{t_{0}+q-1})$ $\displaystyle
f(\bm{x}_{t_{0}+q}|\bm{x}_{t_{0}},\cdots,\bm{x}_{t_{0}+q-1})\cdots
g(\bm{y}_{t}|\bm{x}_{t})d\bm{x}_{t_{0}}\bm{x}_{t_{0}+1}\cdots\bm{x}_{t}.$ (5)

In this paper, we consider a deterministic but unknown change point $t_{0}$.
The goal is to detect the change as quickly as possible subject to false alarm
constraints based on the sequentially observed samples. Let $\mathcal{F}_{t}$
be the $\sigma$-algebra generated by the first $t$ samples
$\bm{y}_{1},\bm{y}_{2},\cdots,\bm{y}_{t}$. A stopping time $\tau$ is a random
variable with the property that for each $t$, the event
$\\{\tau=t\\}\in\mathcal{F}_{t}$. We define the worse-case average detection
delay (WADD) under Lorden’s criterion [41] and the average running length
(ARL) for any stopping time $\tau$ as follows:

$\displaystyle\text{WADD}(\tau)$ $\displaystyle\triangleq\sup_{t_{0}\geq
1}\text{esssup}\mathbb{E}_{t_{0}}\left[(\tau-
t_{0})^{+}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}\right],$
$\displaystyle\text{ARL}(\tau)$
$\displaystyle\triangleq\mathbb{E}_{\infty}[\tau],$ (6)

where $\mathbb{E}_{t_{0}}$ ($\mathbb{E}_{\infty}$) denotes the expectation
under the probability measure $\mathbb{P}_{t_{0}}$ ($\mathbb{P}_{\infty}$).
The goal is to design a stopping time that minimizes the WADD subject to the
constraint on the ARL:

$\displaystyle\inf_{\tau:\text{ARL}(\tau)\geq\gamma}\text{WADD}(\tau),$ (7)

where $\gamma>0$ is a pre-specified threshold.

##  III Universal Lower Bound on WADD

For convenience, we first present results for the first-order AR model to
illustrate our approach. We will then show the generalization to any $q$-th
order AR model. We note that any $q$-th order AR model can be equivalently
converted to a first-order AR model (see details in Appendix A).

Specifically, the first-order AR model is defined as:

$\displaystyle\bm{x}_{t}=\bm{A}\bm{x}_{t-1}+\bm{\omega}_{t},$ (8)

where $\bm{A}\in\mathbb{R}^{K\times K}$ is the matrix coefficient and is
invertible. The initial disturbance signal $\bm{x}_{t_{0}}$ is assumed to be
Gaussian with probability density $f_{1}(\bm{x}_{t_{0}})$ and is independent
of the pre-change observations. We assume that the operator norm of $\bm{A}$
is strictly less than 1. This assumption guarantees the stability of this
system [55].

In the following, we derive a universal lower bound on the WADD for any $\tau$
satisfying the false alarm constraint: $\mathbb{E}_{\infty}[\tau]\geq\gamma$.
To derive the universal lower bound, we first prove the following stability
condition [24, 49, 50, 51]: convergence of the log likelihood ratio

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}.$
(9)

Note that for the AR model, the log likelihood ratio in (9) does not admit an
additive form, and thus its convergence analysis is challenging. To overcome
this difficulty, we first introduce a forward variable $\alpha_{t}$, which is
a scaled Gaussian density function and can be updated recursively. We show
that the conditional density
$p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})$ can be represented
using $\bm{y}_{t}$ and the parameters of the forward variable $\alpha_{t}$.
The log likelihood ratio in (9) can then be written in an additive form. To
further prove the convergence of the log likelihood ratio, we construct an
auxiliary Markov chain using $\bm{y}_{t}$ and the parameter of the forward
variable $\alpha_{t}$. We show that the stationary distribution of this
auxiliary Markov chain exists and then show the convergence of the log
likelihood ratio in (9) by applying the ergodic theorem [52, 53]. The
universal lower bound on the WADD can then be derived.

###  III-A Forward Variable

We first introduce a forward variable, which plays an important role in
deriving the universal lower bound on the WADD. Moreover, our proposed
algorithm in Section IV can be recursively updated using the forward variable.

Note that for the i.i.d. setting, the log likelihood ratio of
$\bm{y}_{t_{0}},\cdots,\bm{y}_{t}$ can be written as the sum of the log
likelihood ratio of each individual sample. For the discrete state HMM in
[30], the likelihood ratio function is represented by the $L_{1}$-norm of
products of Markov random matrices and thus the log likelihood ratio has an
additive form. However, the Markov random matrices techniques can not be
applied to our problem since our AR model has uncountably infinite hidden
states.

Observe that $p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})=\int
f_{1}(\bm{x}_{t_{0}})g(\bm{y}_{t_{0}}|\bm{x}_{t_{0}})f(\bm{x}_{t_{0}+1}|\bm{x}_{t_{0}})\cdots
g(\bm{y}_{t}|\bm{x}_{t})d\bm{x}_{t_{0}}\bm{x}_{t_{0}+1}\cdots\bm{x}_{t}$ is
the integration of the products of a sequence of functions. This motivates us
to represent the likelihood ratio function by replacing the Markov random
matrices and the $L_{1}$-norm of matrices in [30] with functions and
integration of functions, respectively. Moreover, we leverage the Gaussian
property of the innovation noise and the measurement noise in our AR model to
further show that these functions can be represented using the parameters of
the Gaussian density function. Therefore, the likelihood ratio function admits
an additive form and can be updated efficiently.

We define the following forward variable to compute the likelihood function
after the change point ($t\geq t_{0}$),

$\displaystyle\alpha_{t}(\bm{x}_{t})=p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t},\bm{x}_{t}).$
(10)

Note that $\alpha_{t}(\bm{x}_{t})$ also depends on
$\bm{y}_{t_{0}},\cdots,\bm{y}_{t}$. For notational convenience, we treated
$\bm{y}_{t_{0}},\cdots,\bm{y}_{t}$ as known parameters here and write
$\alpha_{t}$ as only a function of $\bm{x}_{t}$. We then have that

$\displaystyle
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})=\int\alpha_{t}(\bm{x}_{t})d\bm{x}_{t}.$
(11)

It can be easily verified that the following recursion holds for the forward
variable:

$\displaystyle\alpha_{t+1}(\bm{x}_{t+1})=\int\alpha_{t}(\bm{x}_{t})f(\bm{x}_{t+1}|\bm{x}_{t})g(\bm{y}_{t+1}|\bm{x}_{t+1})d\bm{x}_{t}.$

We first provide a formal definition for the Gaussian function.

######  Definition 1.

A Gaussian function $f^{\prime}(\cdot):\mathbb{R}^{K}\rightarrow\mathbb{R}$ is
a function of the form

$\displaystyle
f^{\prime}(\bm{x})=a\exp\Big{(}-\frac{1}{2}(\bm{x}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})\Big{)},$
(12)

where $a$ is a constant, $\bm{\mu}\in\mathbb{R}^{K}$ is an arbitrary vector
and $\bm{\Sigma}\in\mathbb{R}^{K\times K}$ is a positive definite matrix.

In the following lemma, we show that $\alpha_{t}(\bm{x}_{t})$ is a Gaussian
function of $\bm{x}_{t}$.

######  Lemma 1.

$\alpha_{t}(\bm{x}_{t})$ is a Gaussian function of $\bm{x}_{t}$.

###### Proof.

We prove this by induction. If the distribution of $\bm{x}_{t_{0}}$ is
Gaussian,
$\alpha_{1}(\bm{x}_{t_{0}})=f_{1}(\bm{x}_{t_{0}})g(\bm{y}_{t_{0}}|\bm{x}_{t_{0}})$
is a Gaussian distribution times a coefficient, which is a Gaussian function.
This can be proved by computing the product of two Gaussian density functions.
Assume that for some $c_{t-1},\bm{\Sigma}_{t-1}$ and $\bm{\mu}_{t-1}$,

$\displaystyle\alpha_{t-1}(\bm{x}_{t-1})=\frac{c_{t-1}}{\sqrt{(2\pi)^{K}\det(\bm{\Sigma}_{t-1})}}\exp\big{(}-\frac{1}{2}(\bm{x}_{t-1}-\bm{\mu}_{t-1})^{\top}\bm{\Sigma}_{t-1}^{-1}(\bm{x}_{t-1}-\bm{\mu}_{t-1})\big{)},$
(13)

which is a Gaussian function of $\bm{x}_{t-1}$. We will show that
$\alpha_{t}(\bm{x}_{t})$ is also a Gaussian function of $\bm{x}_{t}$. From the
product and the convolution of Gaussian density functions, we have that

$\displaystyle\alpha_{t}(\bm{x}_{t})=\int\alpha_{t}(\bm{x}_{t-1})f(\bm{x}_{t}|\bm{x}_{t-1})g(\bm{y}_{t}|\bm{x}_{t})d\bm{x}_{t-1}$
$\displaystyle=\frac{c_{t}}{\sqrt{(2\pi)^{K}\det(\bm{\Sigma}_{t})}}\exp\big{(}-\frac{1}{2}(\bm{x}_{t}-\bm{\mu}_{t})^{\top}\bm{\Sigma}_{t}^{-1}(\bm{x}_{t}-\bm{\mu}_{t})\big{)},$
(14)

where

$\displaystyle\bm{\Sigma}_{t}=(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1},$
$\displaystyle\bm{\mu}_{t}=(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\mu}_{t-1}+(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{y}_{t},$
(15)

and

$\displaystyle\frac{c_{t}}{c_{t-1}}=\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})}}\frac{1}{\sqrt{\det((\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I})}}$
$\displaystyle\hskip
14.22636pt\cdot\exp\bigg{(}-\frac{1}{2}\Big{(}(\bm{A}\bm{\mu}_{t-1})^{\top}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1})$
$\displaystyle\hskip
14.22636pt+\bm{y}_{t}^{\top}\bm{y}_{t}-\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1})+\bm{y}_{t}\big{)}^{\top}$
$\displaystyle\hskip
14.22636pt\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1})+\bm{y}_{t}\big{)}\Big{)}\bigg{)}.$
(16)

Moreover, from the update rule of $\bm{\Sigma}_{t}$, it can be seen that
$\bm{\Sigma}_{t}$ is positive definite. Therefore, $\alpha_{t}(\bm{x}_{t})$ is
a scaled Gaussian distribution of $\bm{x}_{t}$ with coefficient $c_{t}$. ∎

With the forward variable, the conditional density
$p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})$ can be written as
follows:

$\displaystyle
p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})=\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}$
$\displaystyle=\frac{\int\alpha_{t}(\bm{x}_{t})d\bm{x}_{t}}{\int\alpha_{t-1}(\bm{x}_{t-1})d\bm{x}_{t-1}}=\frac{c_{t}}{c_{t-1}},$
(17)

where the last equality is due to the fact that $\alpha_{t}(\bm{x}_{t})$ is a
scaled Gaussian density function of $\bm{x}_{t}$ with coefficient $c_{t}$.
Note that before the change, the observations are i.i.d.. Therefore, it
suffices to consider the additive form for the post-change log likelihood
function. We have the following additive form for the log likelihood function
$\log p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})$:

$\displaystyle\log p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})$
$\displaystyle=\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}+\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-2})}+\cdots+\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\bm{y}_{t_{0}+1})}{p_{t_{0}}(\bm{y}_{t_{0}})}+\log
p_{t_{0}}(\bm{y}_{t_{0}})$
$\displaystyle=\sum_{i=t_{0}}^{t}\frac{c_{i}}{c_{i-1}},$ (18)

where $\frac{c_{t_{0}}}{c_{t_{0}-1}}=\log p_{t_{0}}(\bm{y}_{t_{0}})$. As shown
in (III-A), $\frac{c_{t}}{c_{t-1}}$ is a function of
$\bm{y}_{t},\bm{\mu}_{t-1},\bm{\Sigma}_{t-1}$. Therefore, the log likelihood
function can be written in an additive form using the parameters of the
forward variable.

It can be seen that at each time step, we only need to compute
$c_{t},\bm{\mu}_{t},\bm{\Sigma}_{t}$, and the forward variable $\alpha_{t}$
can be fully recovered. Moreover, $c_{t},\bm{\mu}_{t},\bm{\Sigma}_{t}$ can be
updated recursively, and thus the likelihood function
$p_{t_{0}}(\bm{y}_{1},\cdots,\bm{y}_{t})$ can be computed efficiently using
the forward variable $\alpha_{t}$.

###  III-B Universal Lower Bound on WADD

The general theory for QCD with non-i.i.d. samples relies on the assumption
that the normalized log likelihood ratio between the post- and pre-change
distributions converges to some finite and positive number [24, 49, 50, 51].
For the AR model, whether such assumption holds remains unexplored. Here, we
first show that as $t\rightarrow\infty$, the limit of
$\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
exists and is positive, denoted by $\mathcal{K}$, and further provide the
universal lower bound on WADD in the following theorem. The explicit
expression of $\mathcal{K}$ will be provided later after we introduce
necessary notations.

######  Theorem 1.

We have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}=\mathcal{K}$
(19)

almost surely under $\mathbb{P}_{t_{0}}$ with $\mathcal{K}>0$. Moreover, as
$\gamma\rightarrow\infty$,

$\displaystyle\inf_{\tau:\text{ARL}(\tau)\geq\gamma}\text{WADD}(\tau)\geq\frac{\log\gamma}{\mathcal{K}}(1+o(1)).$
(20)

Since $\bm{y}_{t_{0}},\bm{y}_{t_{0}+1},\cdots,\bm{y}_{t}$ are dependent, the
Law of Large Number used in the i.i.d. setting is not applicable, even if the
log likelihood function has an additive form in (III-A). We develop a novel
approach to show the convergence of
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
under $\mathbb{P}_{t_{0}}$.

Observe that $\frac{c_{t}}{c_{t-1}}$ is a function of
$\bm{y}_{t},\bm{\mu}_{t-1},\bm{\Sigma}_{t-1}$ and
$\bm{y}_{t},\bm{\mu}_{t},\bm{\Sigma}_{t}$ are dependent over time. This
motivates us to formulate a new Markov chain using
$\bm{y}_{t},\bm{\mu}_{t},\bm{\Sigma}_{t}$ and apply the ergodic theorem [52,
53] to characterize the convergence of
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
under $\mathbb{P}_{t_{0}}$.

From the update rule of $\bm{\Sigma}_{t}$ in (III-A), it can be seen that
$\bm{\Sigma}_{t}$ does not depend on $\bm{y}_{t}$ and thus is deterministic.
We then show that $\bm{\Sigma}_{t}$ converges as $t\rightarrow\infty$ and let
$\bm{\Sigma}^{*}\triangleq\lim_{t\rightarrow\infty}\bm{\Sigma}_{t}$. Let
$\bm{\mu}^{*}_{t}$ be the parameter of the forward variable when
$\bm{\Sigma}_{t_{0}}=\bm{\Sigma}^{*}$. We formulate an auxiliary Markov chain
using $\bm{y}_{t},\bm{\mu}^{*}_{t}$ as an intermediate step to prove the
convergence of the log likelihood ratio. Specifically, denote by
$p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})$ the likelihood function when
$\bm{\Sigma}_{t_{0}}=\bm{\Sigma}^{*}$. With the auxiliary Markov chain
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$, we show that
$\lim_{t\rightarrow\infty}\frac{1}{t}p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})$
converges under $\mathbb{P}_{t_{0}}$. The convergence of
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
under $\mathbb{P}_{t_{0}}$ is then proved by showing that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\big{(}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})-\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})\big{)}=0$ (21)

under $\mathbb{P}_{t_{0}}$ almost surely.

We first show that $\bm{\Sigma}_{t}$ converges in the following Lemma.

######  Lemma 2.

As $t\rightarrow\infty$, $\bm{\Sigma}_{t}$ converges to $\bm{\Sigma}^{*}$.

###### Proof.

The proof can be found in Appendix B. ∎

To construct our auxiliary Markov chain, let
$\bm{\Sigma}_{t_{0}}=\bm{\Sigma}^{*}$. Then
$\bm{\Sigma}_{t}=\bm{\Sigma}^{*},\forall t\geq t_{0}$. It follows that
$\bm{\mu}^{*}_{t}$ can be updated recursively as follows:

$\displaystyle\bm{\mu}^{*}_{t}$
$\displaystyle=(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\mu}_{t-1}^{*}+(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{y}_{t}$
$\displaystyle=(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\mu}_{t-1}^{*}+\bm{\Sigma}^{*}\bm{y}_{t}.$
(22)

In the following, we denote the stationary distribution of a Markov chain in
the post-change phase by $\pi$. For examples, $\pi(\bm{x}_{t})$ is the
stationary distribution of $\\{\bm{x}_{t}\\}_{t=t_{0}}^{\infty}$,
$\pi(\bm{x}_{t},\bm{y}_{t})$ is the stationary distribution of
$\\{\bm{x}_{t},\bm{y}_{t}\\}_{t=t_{0}}^{\infty}$. We first show that
$\\{\bm{x}_{t},\bm{y}_{t}\\}_{t=t_{0}}^{\infty}$ has a unique stationary
distribution in the following lemma, which will be used to formulate and
characterize the property of the auxiliary Markov chain.

######  Lemma 3.

The HMM $\\{\bm{x}_{t},\bm{y}_{t}\\}_{t=t_{0}}^{\infty}$ is $\pi$-irreducible.
Moreover, $\pi(\bm{x}_{t})g(\bm{y}_{t}|\bm{x}_{t})$ is the unique stationary
distribution of $\\{\bm{x}_{t},\bm{y}_{t}\\}_{t=t_{0}}^{\infty}$, where
$\pi(\bm{x}_{t})$ is a Gaussian distribution with mean $\bm{0}$ and covariance
matrix
$\bm{\Sigma}=\sum_{i=0}^{\infty}(\bm{A}^{\top})^{i}\bm{R}_{\omega}\bm{A}^{i}$.

###### Proof.

The detailed proof can be found in Appendix C. ∎

In the following lemma, we show that
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t=t_{0}}^{\infty}$ is a Markov chain and
the stationary distribution of
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t=t_{0}}^{\infty}$ exists and is unique.

######  Lemma 4.

$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t=t_{0}}^{\infty}$ is a Markov chain.
Moreover,
$\int\pi(\bm{\mu}_{t-1}^{*})p_{t_{0}}(\bm{y}_{t},\bm{\mu}^{*}_{t}|\bm{\mu}^{*}_{t-1})d\bm{\mu}^{*}_{t-1}$
is the unique stationary distribution of
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t=t_{0}}^{\infty}$, where
$\pi(\bm{\mu}^{*}_{t-1})$ is the stationary distribution of
$\bm{\mu}^{*}_{t-1}$ and is guaranteed to exist.

###### Proof.

The proof can be found in Appendix D. ∎

With the auxiliary Markov chain
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ and its stationary
distribution, we are ready to show the convergence of
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
and prove Theorem 1.

###### Proof sketch of Theorem 1.

We first consider the auxiliary Markov chain
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ and show that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]$
(23)

almost surely under $\mathbb{P}_{t_{0}}$ by applying the ergodic theorem [52,
53] to the Markov chain
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t=t_{0}}^{\infty}$, where
$h(\bm{\mu},\bm{y})$ is a quadratic function of $\bm{\mu},\bm{y}$ and
$(\bm{y},\bm{\mu})$ follows the stationary distribution of the auxiliary
Markov chain $\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$. The
explicit expression of $h(\bm{\mu},\bm{y})$ is as follows

$\displaystyle
h(\bm{\mu},\bm{y})=\log\Big{(}\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})}}\frac{1}{\sqrt{\det((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I})}}\Big{)}$
$\displaystyle\hskip
5.69046pt-\frac{1}{2}\Big{(}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}-\bm{\Sigma}^{*}\bm{y}))^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}-\bm{\Sigma}^{*}\bm{y}))$
$\displaystyle\hskip
5.69046pt+\bm{y}^{\top}\bm{y}-\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}-\bm{\Sigma}^{*}\bm{y}))+\bm{y}\big{)}^{\top}$
$\displaystyle\hskip
5.69046pt\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}-\bm{\Sigma}^{*}\bm{y}))+\bm{y}\big{)}\Big{)}.$
(24)

Moreover, if we let $\bm{y}=\bm{y}_{t}$ and $\bm{\mu}=\bm{\mu}^{*}_{t}$, we
have that
$h(\bm{\mu}^{*}_{t},\bm{y}_{t})=p^{*}_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})$.

We then show that $\lim_{t\rightarrow\infty}\frac{1}{t}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]$
under $\mathbb{P}_{t_{0}}$ almost surely by showing that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\big{(}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})-\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})\big{)}=0$ (25)

under $\mathbb{P}_{t_{0}}$ almost surely.

Since the observations are independent before the change point, we have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log
p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\log
p_{\infty}(\bm{y}_{i}).$

We then have that under $\mathbb{P}_{t_{0}}$,

$\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\log
p_{\infty}(\bm{y}_{i})=\mathbb{E}_{\pi}[\log p_{\infty}(\bm{y})]$

almost surely from the ergodic theorem of Markov chain [52, 53].

Let

$\displaystyle\mathcal{K}=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]-\mathbb{E}_{\pi}[\log
p_{\infty}(\bm{y})].$ (26)

It then follows that for any initial state $\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}$,

$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}=\mathcal{K}$

under $\mathbb{P}_{t_{0}}$ almost surely. Therefore, for any $\eta>0$,

$\displaystyle\lim_{t\rightarrow\infty}$ $\displaystyle\sup_{t_{0}\geq
1}\text{esssup}\mathbb{P}_{t_{0}}\Bigg{\\{}\max_{k\leq
t}\sum_{i=t_{0}}^{t_{0}+k-1}\log\frac{p_{t_{0}}(\bm{y}_{i}|\bm{y}_{t_{0}},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i})}\geq\mathcal{K}(1+\eta)t\big{|}\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}\Bigg{\\}}=0.$
(27)

Then (20) follows from [24, Theorem 1]. The full proof can be found in
Appendix E. ∎

######  Remark 1.

Results in [24] assumes
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
exists under $\mathbb{P}_{t_{0}}$ while in our results, we prove its existence
and characterize its value.

###  III-C $q$-th Order AR Models

In this section, we show that our results for first-order AR models can be
generalized to any $q$-th order AR models.

In the following theorem, we show that for any $q$-th order AR model, as in
Section III-B, as $t\rightarrow\infty$, the limit of
$\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
exists and is positive, denoted by $\widetilde{\mathcal{K}}$. The expression
of $\widetilde{\mathcal{K}}$ can be derived similarly as in (26). The
universal lower bound on WADD then follows from [24, Theorem 1].

######  Theorem 2.

For a $q$-th order AR model, we have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}=\widetilde{\mathcal{K}}$
(28)

almost surely under $\mathbb{P}_{t_{0}}$ where $\widetilde{\mathcal{K}}>0$.
Moreover, as $\gamma\rightarrow\infty$,

$\displaystyle\inf_{\tau:\text{ARL}(\tau)\geq\gamma}\text{WADD}(\tau)\geq\frac{\log\gamma}{\widetilde{\mathcal{K}}}(1+o(1)).$
(29)

The proof of Theorem 2 is similar to the proof of Theorem 1, and the idea is
to convert the $q$-th order AR model to a first-order AR model, and then apply
the proof of Theorem 1. Below, we provide a proof sketch.

###### Proof Sketch.

Note that for a $q$-th order AR model, it can be converted to a first-order AR
model (See Appendix A). Let $\lfloor x\rfloor$ denote the greatest integer
less than or equal to $x$. We have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\bigg{(}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1})}+\log\frac{p_{t_{0}}(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1})}{p_{\infty}(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1})}$
$\displaystyle\hskip
14.22636pt+\cdots+\log\frac{p_{t_{0}}(\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)},\cdots,\bm{y}_{t_{0}+q\lfloor\frac{t}{q}\rfloor-1}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)-1})}{p_{\infty}(\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)},\cdots,\bm{y}_{t_{0}+q\lfloor\frac{t}{q}\rfloor-1}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)-1})}\bigg{)}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}\Big{(}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1}),(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}),\cdots,(\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)},\cdots,\bm{y}_{t_{0}+q\lfloor\frac{t}{q}\rfloor-1})\Big{)}}{p_{\infty}\Big{(}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1}),(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}),\cdots,(\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)},\cdots,\bm{y}_{t_{0}+q\lfloor\frac{t}{q}\rfloor-1})\Big{)}}$
$\displaystyle=\lim_{t^{\prime}\rightarrow\infty}\frac{1}{qt^{\prime}}\log\frac{p_{t_{0}}\Big{(}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1}),(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}),\cdots,(\bm{y}_{t_{0}+q(t^{\prime}-1)},\cdots,\bm{y}_{t_{0}+qt^{\prime}-1})\Big{)}}{p_{\infty}\Big{(}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1}),(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}),\cdots,(\bm{y}_{t_{0}+q(t^{\prime}-1)},\cdots,\bm{y}_{t_{0}+qt^{\prime}-1})\Big{)}}$
$\displaystyle=\widetilde{\mathcal{K}}$ (30)

almost surely under $\mathbb{P}_{t_{0}}$, where the last equality is from (19)
and the fact that
$\\{(\bm{y}_{t_{0}+qi},\cdots,\bm{y}_{t_{0}+q(i+1)-1})\\}_{i=0}^{\infty}$
follows a first-order AR model. For the quickest change detection in $q$-th
order AR models, we then have that as $\gamma\rightarrow\infty$,

$\displaystyle\inf_{\tau:\text{ARL}(\tau)\geq\gamma}\text{WADD}(\tau)\geq\frac{\log\gamma}{\widetilde{\mathcal{K}}}(1+o(1)).$
(31)

∎

##  IV Asymptotically Optimal Stopping Time

In this section, we first present the algorithm and its optimality results for
the first-order AR model. We then show the generalization to the $q$-th order
AR model.

###  IV-A First-Order AR Model

The CuSum algorithm based on the generalized likelihood ratio (GLR) approach
has been widely used for QCD problems. For the AR model, the GLR statistic is
defined as follows:

$\displaystyle W_{t}=\max_{1\leq k\leq
t}\sum_{i=k}^{t}\log\frac{p_{k}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}.$
(32)

The CuSum algorithm [24] can then be designed:

$\displaystyle\tau_{c}=\inf\Big{\\{}t:W_{t}\geq c\Big{\\}}.$ (33)

In [24], a special non-i.i.d. case was studied where the post-change
distribution does not depend on the change point $t_{0}$, i.e.,
$p_{k}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})$ in (32) does not depend on
$k$. Under the assumption that
$\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\log\frac{p_{t_{0}}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}$
converges in probability under $\mathbb{P}_{t_{0}}$ to some positive constant
$\mathcal{K}$, the asymptotic optimality of the CuSum algorithm was proved.
The non-i.i.d. case under the Bayesian setting was studied in [56], where the
post-change distribution does not depend on the change point $t_{0}$ as in
[24]. The asymptotic optimality of Shiryaev procedure was established under
the same stability assumption.

If the post-change distribution does not depend on the change-point, then the
CuSum algorithm can be updated recursively and is computationally efficient.
However, for our AR model, under $\mathbb{P}_{t_{0}}$,
$p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})$ depends on $t_{0}$
for $t\geq t_{0}$. At each time $t$, we need to update
$p_{k}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})$ for every $1\leq k\leq t$,
the complexity of which scales with $t$, which is not practically feasible.

In this section, we propose a computationally efficient Ergodic CuSum
algorithm and further show that it is asymptotically optimal.

Motivated by the fact that the WADD is lower bounded by
$\frac{\log\gamma}{\mathcal{K}}(1+o(1))$, we aim to find a computationally
efficient test statistic that has a positive drift of $\mathcal{K}$ under the
post-change distribution and has a negative drift under the pre-change
distribution. Define the likelihood ratio of the first $t$ observations when
the change point $t_{0}=1$ by

$\displaystyle
L_{t}=\frac{p_{1}(\bm{y}_{1},\cdots,\bm{y}_{t})}{p_{\infty}(\bm{y}_{1},\cdots,\bm{y}_{t})},$
(34)

and let $L_{0}=1$. For any $t_{0}\geq 1$, we have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{1}(\bm{y}_{1},\cdots,\bm{y}_{t})}{p_{\infty}(\bm{y}_{1},\cdots,\bm{y}_{t})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\bigg{(}\log\frac{p_{1}(\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}+\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}\bigg{)}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}.$
(35)

Note that different sample trajectories $\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}$
lead to different values of $\bm{\mu}_{t_{0}}$. However, the ergodic theorem
of Markov chain [52, 53] implies that the convergence of
$\frac{1}{t}\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}$
does not depend on $\bm{\mu}_{t_{0}}$ and thus does not depend on the sample
trajectory $\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}$. Therefore, we have that for
any $\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}$, under $\mathbb{P}_{t_{0}}$ almost
surely

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}$
$\displaystyle=\mathcal{K},$ (36)

where the first equality is due to the fact that the convergence of the
ergodic Markov chain doesn’t depend on the initial state [52, 53]. Therefore,
under the post-change distribution, $L_{t}$ has a positive drift $\mathcal{K}$
on average. Under the pre-change distribution, we have that
$\mathbb{E}_{\infty}\big{[}\log\frac{p_{1}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})}{p_{\infty}(\bm{y}_{t})}\big{]}=-D\big{(}p_{\infty}(\bm{y}_{t})||p_{1}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})\big{)}\leq
0$ where the equality holds when
$p_{1}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})=p_{\infty}(\bm{y}_{t})$
almost surely. Therefore, under the pre-change distribution, $L_{t}$ has a
negative drift as long as
$p_{1}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})\neq p_{\infty}(\bm{y}_{t})$.

Motivated by these facts, we define the Ergodic CuSum statistic

$\displaystyle S_{t}$ $\displaystyle=\max_{0\leq i\leq t}(\log L_{t}-\log
L_{i})$ $\displaystyle=\max\big{(}0,S_{t-1}+\log L_{t}-\log L_{t-1}\big{)}.$
(37)

The Ergodic CuSum algorithm is then defined as

$\displaystyle\tau_{c}^{*}=\inf\\{t:S_{t}\geq c\\}.$ (38)

At each time $t$, we only need to compute $L_{t}$. It can be easily verified
using (III-A) that $L_{t}$ admits an additive form. Moreover, $L_{t}$ can be
recursively updated using the forward variable as shown in Lemma 1.
Specifically, let $\bm{\mu}_{0},\bm{\Sigma}_{0}$ be the parameters of the
initial distribution of $\bm{x}_{0}$. We have that
$\bm{\mu}_{t},\bm{\Sigma}_{t}$ can be updated recursively according to (III-A)
for each $t$. We then have that
$L_{t}=\sum_{i=1}^{t}\frac{c_{i}}{c_{i-1}}=L_{t-1}+\frac{c_{t}}{c_{t-1}}$,
where $\frac{c_{t}}{c_{t-1}}$ can be computed using
$\bm{\mu}_{t},\bm{\Sigma}_{t}$ as shown in (III-A).

We note that $S_{t}$ is not the actual generalized likelihood ratio, and is
different from $W_{t}$ in $\tau_{c}$. However, at each time $t$,
$\log\frac{p_{1}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})}{p_{\infty}(\bm{y}_{t})}$
can still be viewed as a log likelihood ratio of the current sample
$\bm{y}_{t}$. Therefore, the ARL lower bound can be derived following the
proof of the ARL lower bound in [24] for general non-i.i.d. case. Since the
convergence of
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}$
under $\mathbb{P}_{t_{0}}$ does not depend on the initial state
$\bm{\mu}_{t_{0}},\bm{y}_{t_{0}}$, it can be shown that $\tau_{c}^{*}$ is
asymptotically optimal for (7).

In the following theorem, we show 1) the ARL lower bound of $\tau_{c}^{*}$ and
2) the WADD upper bound of $\tau_{c}^{*}$.

######  Theorem 3.

1) Let $c=\log\gamma$ in (38), then
$\mathbb{E}_{\infty}[\tau_{c}^{*}]\geq\gamma$; and 2) as
$\gamma\rightarrow\infty$,
$\text{WADD}(\tau_{c}^{*})\leq\frac{\log\gamma}{\mathcal{K}}(1+o(1))$.

###### Proof.

The proof can be found in Appendix F. ∎

Based on Theorem 1 and Theorem 3, we establish the asymptotic optimality of
$\tau_{c}^{*}$ in the following theorem.

######  Theorem 4.

$\tau_{c}^{*}$ is asymptotically optimal.

###### Proof.

By Theorem 1 and Theorem 3, we establish the asymptotic optimality of
$\tau_{c}^{*}$. ∎

The Ergodic CuSum algorithm in (38) is computationally efficient and
asymptotically optimal for detecting the change in the AR model in Section II.

###  IV-B $q$-th Order AR Model

We first convert the post-change $q$-th order AR model equivalently into a
first-order AR model (see Appendix A). We then partition the sequence of the
observations into a sequence of non-overlapping blocks with size $q$.
Specially, define $\tilde{y}_{t}=(y_{(t-1)q+1},\ldots,y_{tq})$, for
$t=1,2,\ldots$. We then apply our Ergodic Cusum algorithm in (38) on the
sequence of $\\{\tilde{y}_{t}\\}_{t=1}^{\infty}$. We show that our Ergodic
Cusum algorithm is asymptotically optimal for the problem of QCD in $q$-th
order AR models.

For the ARL lower bound and WADD upper bound of $\tau^{*}_{c}$, we have the
following theorem.

######  Theorem 5.

For the QCD problem in $q$-th order AR models, consider $\tau^{*}_{c}$ applied
on $\\{\tilde{y}_{t}\\}_{t=1}^{\infty}$. 1) Let $c=\log\gamma$ in (38), then
$\mathbb{E}_{\infty}[\tau_{c}^{*}]\geq\gamma$. 2) As
$\gamma\rightarrow\infty$,
$\text{WADD}(\tau_{c}^{*})\leq\frac{\log\gamma}{\widetilde{\mathcal{K}}}(1+o(1))$.

The proof of Theorem 5 is similar to the proof of Theorem 3. Here, we only
provide a proof sketch.

###### Proof Sketch.

When the post-change disturbance signal follows a $q$-th order AR model, we
have that for any $\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}$, under
$\mathbb{P}_{t_{0}}$ almost surely

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{1}(\bm{y}_{1},\cdots,\bm{y}_{t})}{p_{\infty}(\bm{y}_{1},\cdots,\bm{y}_{t})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\bigg{(}\log\frac{p_{1}(\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}+\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}\bigg{)}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}\Big{(}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1}),(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}),\cdots,(\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)},\cdots,\bm{y}_{t_{0}+q\lfloor\frac{t}{q}\rfloor-1})\Big{)}}{p_{\infty}\Big{(}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+q-1}),(\bm{y}_{t_{0}+q},\cdots,\bm{y}_{t_{0}+2q-1}),\cdots,(\bm{y}_{t_{0}+q(\lfloor\frac{t}{q}\rfloor-1)},\cdots,\bm{y}_{t_{0}+q\lfloor\frac{t}{q}\rfloor-1})\Big{)}}$
$\displaystyle=\widetilde{\mathcal{K}}.$ (39)

where the third equality is due to the fact that the convergence of
$\frac{1}{t}\log\frac{p_{1}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1})}$
does not depend on the sample trajectory $\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}$
and the last equality is from (III-C) and the fact that
$\\{(\bm{y}_{t_{0}+qi},\cdots,\bm{y}_{t_{0}+q(i+1)-1})\\}_{i=0}^{\infty}$
follows a first-order AR model. Following the same steps as in Theorem 3, we
have that when $c=\log\gamma$, $\mathbb{E}_{\infty}[\tau_{c}^{*}]\geq\gamma$
and
$\text{WADD}(\tau_{c}^{*})\leq\frac{\log\gamma}{\widetilde{\mathcal{K}}}(1+o(1))$
as $\gamma\rightarrow\infty$. ∎

It can be seen that the WADD upper bound of $\tau_{c}^{*}$ matches with the
universal lower bound in Theorem 2 for the $q$-th order AR models. We then
have the following result of asymptotic optimality.

######  Theorem 6.

For the problem of QCD in $q$-th order AR models, $\tau_{c}^{*}$ (applied on
its first-order equivalence) is asymptotically optimal.

###### Proof.

By Theorem 2 and Theorem 5, we establish the asymptotic optimality of
$\tau_{c}^{*}$. ∎

##  V Data-Driven Setting: Online Gradient Ascent CuSum

In this section, we consider the practical data-driven setting where the post-
change parameters are unknown. This is motivated by practical applications
such as detecting cyber-attacks in dynamic systems, where there is no prior
knowledge about the disturbance signal $\bm{x}_{t}$. As shown in Sections III
and IV, the QCD problem in a $q$-th order AR model can be equivalently solved
using a first-order AR model. Therefore, in this section, we focus on the
first-order AR model. Specifically, the post-change parameters
$\bm{A},\bm{R}_{\omega}$ are unknown.

The generalized likelihood ratio test (GLRT) which replaces the unknown post-
change parameter with its maximum likelihood estimate (MLE) is widely used
when there are unknown parameters. Specifically, at each time $t$,
hypothesizing on the change point being $k$, we use samples
$\bm{y}_{k},\cdots,\bm{y}_{t}$ to compute the MLE of the unknown post-change
parameters, i.e.,
$\hat{\theta}_{k,t}=\arg\max_{\theta}p_{k,\theta}(\bm{y}_{k},\cdots,\bm{y}_{t})$,
where $p_{k,\theta}(\bm{y}_{k},\cdots,\bm{y}_{t})$ denotes the likelihood of
$\bm{y}_{k},\cdots,\bm{y}_{t}$ when the change point is $k$ and post-change
parameter is $\theta$. The MLE of the hypothesized change point is then
computed. Specifically, let $S_{\hat{\theta}_{k,t}}$ be the test statistic
with post-change parameters $\hat{\theta}_{k,t}$, and take the maximum of
$S_{\hat{\theta}_{k,t}}$ over all $k$, i.e,
$S_{t}=\max_{k}S_{\hat{\theta}_{k,t}}$. This approach needs to store all the
historical samples $\bm{y}_{1},\cdots,\bm{y}_{t}$ and recompute
$\hat{\theta}_{k,t}$ when there is a new sample, which is not efficient for
memory and computation. Though in practice, a sliding-window approach can be
used, but it is still hard in the non-i.i.d. setting [24].

In this section, we apply an online gradient ascent algorithm [54, 48] to
estimate the unknown parameters $\bm{A},\bm{R}_{\omega}$ and plug them back to
the Ergodic CuSum algorithm to design our online gradient ascent CuSum
algorithm (OGA-CuSum). The OGA-CuSum can be updated using only the most recent
sample and in a recursive way, and thus is much more memory and
computationally efficient.

Recall the definition of $h(\bm{\mu},\bm{y})$ in (III-B). Define
$h_{t}(\bm{\mu},\bm{y})$ by replacing $\bm{\Sigma}^{*}$ with $\bm{\Sigma}_{t}$
in $h(\bm{\mu},\bm{y})$. From Appendix E, letting $\bm{y}=\bm{y}_{t}$ and
$\bm{\mu}=\bm{\mu}_{t}$, we have that
$h_{t}(\bm{\mu}_{t},\bm{y}_{t})=p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})$.
Given the initial $\bm{\mu}_{0},\bm{\Sigma}_{0}$, it then follows that

$\displaystyle\log L_{t}=\sum_{i=1}^{t}h_{i}(\bm{\mu}_{i},\bm{y}_{i})-\log
p_{\infty}(\bm{y}_{i}).$ (40)

Since $\bm{A},\bm{R}_{\omega}$ are unknown, $h_{t}(\bm{\mu}_{t},\bm{y}_{t})$
cannot be computed. We then propose a one-step update rule to efficiently
estimate $\bm{A}$, $\bm{R}_{\omega}$ and further estimate
$h_{t}(\bm{\mu}_{t},\bm{y}_{t})$ at each time $t$.

Denote by $\widehat{\bm{A}}_{t},\widehat{\bm{R}}_{\omega,t}$ the estimate of
$\bm{A}$, $\bm{R}_{\omega}$ at time $t$. Denote by
$\widehat{\bm{\mu}}_{t+1},\widehat{\bm{\Sigma}}_{t+1}$ the parameters of the
forward variable. From (III-A), we have that
$\widehat{\bm{\mu}}_{t+1},\widehat{\bm{\Sigma}}_{t+1}$ can be updated
recursively using $\widehat{\bm{A}}_{t},\widehat{\bm{R}}_{\omega,t}$, i.e.,

$\displaystyle\widehat{\bm{\Sigma}}_{t+1}=(\widehat{\bm{A}}_{t}\widehat{\bm{\Sigma}}_{t}\widehat{\bm{A}}_{t}^{\top}+\widehat{\bm{R}}_{\omega,t})(\widehat{\bm{A}}_{t}\widehat{\bm{\Sigma}}_{t}\widehat{\bm{A}}_{t}^{\top}+\widehat{\bm{R}}_{\omega,t}+\bm{I})^{-1},$
$\displaystyle\widehat{\bm{\mu}}_{t+1}=(\widehat{\bm{A}}_{t}\widehat{\bm{\Sigma}}_{t}\widehat{\bm{A}}_{t}^{\top}+\widehat{\bm{R}}_{\omega,t}+\bm{I})^{-1}\widehat{\bm{A}}_{t}\widehat{\bm{\mu}}_{t}+(\widehat{\bm{A}}_{t}\widehat{\bm{\Sigma}}_{t}\widehat{\bm{A}}_{t}^{\top}+\widehat{\bm{R}}_{\omega,t})(\widehat{\bm{A}}_{t}\widehat{\bm{\Sigma}}_{t}\widehat{\bm{A}}_{t}^{\top}+\widehat{\bm{R}}_{\omega,t}+\bm{I})^{-1}\bm{y}_{t+1}.$
(41)

Define $\widehat{h}_{t+1}(\widehat{\bm{\mu}}_{t+1},\bm{y}_{t+1})$ by replacing
$\bm{A},\bm{R}_{\omega},\bm{\mu}_{t+1},\bm{\Sigma}_{t+1}$ in
$h_{t+1}(\bm{\mu}_{t+1},\bm{y}_{t+1})$ with
$\widehat{\bm{A}}_{t},\widehat{\bm{R}}_{\omega,t},\widehat{\bm{\mu}}_{t+1},\widehat{\bm{\Sigma}}_{t+1}$.
We iteratively update the estimate of the parameters when there is a new
observation based on the maximum likelihood principle using gradient ascent.
We note that $h_{t}(\bm{\mu}_{t},\bm{y}_{t})-\log p_{\infty}(\bm{y}_{t})$ is
the log-likelihood ratio of the observation $\bm{y}_{t}$ at time $t$, i.e.,
$h_{t}(\bm{\mu}_{t},\bm{y}_{t})-\log
p_{\infty}(\bm{y}_{t})=\frac{p_{1}(\bm{y}_{t}|\bm{y}_{1},\cdots,\bm{y}_{t-1})}{p_{\infty}(\bm{y}_{t})}$.
Therefore, based on the maximum likelihood principle, the online gradient
ascent estimator is updated as follows

$\displaystyle\widehat{\bm{A}}_{t}=\widehat{\bm{A}}_{t-1}+\beta\nabla_{\bm{A}}\widehat{h}_{t}(\widehat{\bm{\mu}}_{t},\bm{y}_{t}),$
$\displaystyle\widehat{\bm{R}}_{\omega,t}=\text{Proj}\big{(}\widehat{\bm{R}}_{\omega,t-1}+\beta\nabla_{\bm{R}_{\omega}}\widehat{h}_{t}(\widehat{\bm{\mu}}_{t},\bm{y}_{t})\big{)},$
(42)

where $\beta$ is a pre-specified step-size and Proj$(\bm{X})$ is an operator
that projects a matrix $\bm{X}$ to the set of positive definite matrices to
guarantee that $\widehat{\bm{R}}_{\omega,t}$ is positive definite. Define the
eigenvalue decomposition $\bm{X}=\sum_{i=1}^{n}\lambda_{i}v_{i}v_{i}^{\top}$,
where $\lambda_{i}$ is the eigenvalue of $\bm{X}$ and $v_{i}$ is its
corresponding eigenvector. We define
Proj$(\bm{X})=\sum_{i=1}^{n}\max\\{\lambda_{i},\epsilon\\}v_{i}v_{i}^{\top}$
[57], where $\epsilon>0$ is chosen to guarantee that the eigenvalue of
Proj$(\bm{X})$ is at least $\epsilon$. With
$\widehat{\bm{A}}_{t},\widehat{\bm{R}}_{\omega,t}$, we can update
$\widehat{\bm{\mu}}_{t+1},\widehat{\bm{\Sigma}}_{t+1}$ according to (V). We
further plug them back to (IV-A) and (40) to compute $\widehat{S}_{t+1}$,
which serves as the estimate of $S_{t+1}$. When $\widehat{S}_{t+1}<0$, we
claim that the change hasn’t occurred. Therefore, we ignore the previous
samples and reset our parameters
$\widehat{\bm{A}}_{t+1},\widehat{\bm{R}}_{\omega,t+1}$. The OGA-CuSum
algorithm is then defined as follows

$\displaystyle\tau_{\text{OGA}}=\inf\\{t:\widehat{S}_{t}\geq c\\}.$ (43)

We summarize the algorithm in Algorithm 1.

In general, it is hard to obtain theoretical optimality performance guarantees
for data-driven approaches. Nevertheless, we provide a lower bound on the ARL
for our OGA-CuSum algorithm, so that a threshold can be chosen analytically to
control the false alarm in practice. We also provide simulation results to
demonstrate the good performance of our algorithm.

In the following theorem, we present a lower bound on the ARL for our OGA-
CuSum algorithm.

######  Theorem 7.

Let $c=\log\gamma$ in Algorithm 1, then
$\mathbb{E}_{\infty}[\tau_{\text{OGA}}]\geq\gamma$.

###### Proof.

The proof can be found in Appendix G. ∎

Algorithm 1 Online gradient ascent CuSum algorithm

Initial parameters
$c,\beta,\widehat{\bm{\mu}}_{0},\widehat{\bm{\Sigma}}_{0},\widehat{\bm{A}}_{0},\widehat{\bm{R}}_{\omega,0}$,
$t\leftarrow 0,\widehat{L}_{0}\leftarrow 1,\widehat{S}_{0}\leftarrow 0$

while $\widehat{S}_{t}<c$ do

Take a new observation $\bm{y}_{t+1}$

Update $\widehat{\bm{\mu}}_{t+1},\widehat{\bm{\Sigma}}_{t+1}$ according to (V)

$\log\widehat{L}_{t+1}\leftarrow\log\widehat{L}_{t}+\widehat{h}_{t+1}(\widehat{\bm{\mu}}_{t+1},\bm{y}_{t+1})-\log
p_{\infty}(\bm{y}_{t+1})$

$\widehat{S}_{t+1}\leftarrow\widehat{S}_{t}+\log\widehat{L}_{t+1}-\log\widehat{L}_{t}$

if $\widehat{S}_{t+1}<0$ then

$\widehat{\bm{A}}_{t+1}\leftarrow\widehat{\bm{A}}_{0}$

$\widehat{\bm{R}}_{\omega,t+1}\leftarrow\widehat{\bm{R}}_{\omega,0}$

$\widehat{S}_{t+1}\leftarrow 0$

$\widehat{L}_{t+1}\leftarrow 1$

else

$\widehat{\bm{A}}_{t+1}=\widehat{\bm{A}}_{t}+\beta\nabla_{\bm{A}}\widehat{h}_{t+1}(\widehat{\bm{\mu}}_{t+1},\bm{y}_{t+1})$

$\widehat{\bm{R}}_{\omega,t+1}=\text{Proj}\big{(}\widehat{\bm{R}}_{\omega,t}+\beta\nabla_{\bm{R}_{\omega}}\widehat{h}_{t+1}(\widehat{\bm{\mu}}_{t+1},\bm{y}_{t+1})\big{)}$

end if

$t\leftarrow\ t+1$

end while

##  VI Simulation Results

In this section, we provide some numerical results to demonstrate the
performance of our Ergodic CuSum algorithm and OGA-CuSum algorithm.

###  VI-A Model-Based Setting

In this section, we consider the setting where the parameters of the
disturbance signal are known. We compare the Ergodic CuSum algorithm
$\tau_{c}^{*}$ with a stationary CuSum algorithm designed using the stationary
distribution of $\bm{y}$. The stationary CuSum algorithm uses the CuSum
algorithm that detects a change from $p_{\infty}$ to the stationary
distribution $\pi$:

$\displaystyle\tau_{s}=\inf\Big{\\{}t:\max_{1\leq k\leq
t}\sum_{i=k}^{t}\log\frac{\pi(\bm{y}_{i})}{p_{\infty}(\bm{y}_{i})}\geq
c\Big{\\}}.$ (44)

We consider three different cases. For the first case, we set
$\bm{A}=[[0.7,0.4],[0.2,0.6]]$ and $\bm{R}_{\omega}=[[1,0.5],[0.5,1]]$. For
the second case, we let $\bm{A}$ be a $10\times 10$ matrix and
$\bm{R}_{\omega}=\bm{I}$. For the third case, we consider a $q$-th order AR
model. We set
$\bm{A}_{1}=[[0.4,0.3],[0.2,0.1]],\bm{A}_{2}=[[0.3,0.2],[0.1,0.2]]$ and
$\bm{R}_{\omega}=\bm{I}$. Note that our Ergodic CuSum algorithm can be easily
generalized to the $q$-th order AR model by reformulating the $q$-th order AR
model to a first-order AR model. In Fig. 3, Fig. 3 and Fig. 3, we plot the
WADD as a function of ARL. It can be seen from Fig. 3, Fig. 3 and Fig. 3 that
with the same constraint on the ARL, our Ergodic CuSum algorithm has a lower
detection delay. Therefore, our Ergodic CuSum algorithm outperforms the
stationary CuSum algorithm under both the low-dimensional setting and the
high-dimensional setting. This further implies that for detecting changes in
the AR model, simply applying the stationary distribution and ignoring the
dependence among samples may not perform well. Moreover, the relationship
between WADD and log of the WARL is linear, which validates our theoretical
analysis.

Figure 1: Comparison of the two algorithms: case 1.

Figure 2: Comparison of the two algorithms: case 2.

Figure 3: Comparison of the two algorithms: case 3.

###  VI-B Data-Driven Setting

In this section, we focus on the data-driven setting. We first plot the
estimation errors of $\widehat{\bm{A}}_{t}$, $\widehat{\bm{R}}_{\omega,t}$ as
functions of number of samples under the pre- and post-change distribution,
where the estimation errors are defined as the Frobenius norm of
$\widehat{\bm{A}}_{t}-\bm{A}$ and
$\widehat{\bm{R}}_{\omega,t}-\bm{R}_{\omega}$. Under the post-change
distribution, we set $\bm{A}=[[0.7,0.4],[0.2,0.6]]$ and
$\bm{R}_{\omega}=[[1,0.5],[0.5,1]]$. We note that the pre-change distribution
is equivalent to the AR model with $\bm{A}=\bm{0}$ and
$\bm{R}_{\omega}=\bm{0}$. It can be seen from Fig. 5, Fig. 5, Fig. 7 and Fig.
7 that the estimated parameters $\widehat{\bm{A}}_{t}$,
$\widehat{\bm{R}}_{\omega,t}$ converge to the true parameters under the pre-
and post-change distribution, which demonstrates that our OGA-CuSum scheme
provides reliable estimates for the true parameters of the disturbance signal.

Figure 4: Convergence of $\widehat{\bm{A}}_{t}$ under the pre-change
distribution.

Figure 5: Convergence of $\widehat{\bm{A}}_{t}$ under the post-change
distribution.

Figure 6: Convergence of $\widehat{\bm{R}}_{\omega,t}$ under the pre-change
distribution.

Figure 7: Convergence of $\widehat{\bm{R}}_{\omega,t}$ under the post-change
distribution

We then compare our OGA-CuSum algorithm with an existing approach in [15]. In
[15], a GLRT based data-driven algorithm was proposed to detect disturbance
signals in the AR model. We also plot the curve of the Ergodic CuSum algorithm
under the model-based setting, which serves as a lower bound for the data-
driven setting. We use the same parameters as in Section VI-A and plot the
WADD as a function of ARL. From Fig. 10, Fig. 10 and Fig. 10, it can be seen
that our OGA-CuSum algorithm outperforms the GLRT based algorithm in [15].
Moreover, the performance of the OGA-CuSum algorithm is close to the
performance of the Ergodic CuSum algorithm. Therefore, our OGA-CuSum algorithm
has a small performance loss compared with the asymptotically optimal
algorithm under the model-based setting.

Figure 8: Comparison of the three algorithms: case 1.

Figure 9: Comparison of the three algorithms: case 2.

Figure 10: Comparison of the three algorithms: case 3.

##  VII Conclusion

In this paper, we studied the QCD problem in AR models under the model-based
setting and the data-driven setting. For the model-based setting, we proposed
a novel Ergodic CuSum algorithm. By introducing the forward variable for
general state HMMs, our Ergodic CuSum algorithm can be implemented
efficiently. We further showed that our algorithm is asymptotically optimal
under the Lorden’s criterion [41] based on the ergodic theorem [52, 53]. For
the data-driven setting, we proposed an OGA-CuSum algorithm, which replaces
the unknown parameters in the Ergodic CuSum algorithm with their estimates
based on the online convex optimization. We provided simulation results to
demonstrate the performance of our algorithm. Our approaches provide useful
insights for general detection problems for AR models and general state HMMs.

##  Appendix A Convert $q$-th order AR Model to First-order AR Model

Consider the case with $q=2$. Let
$\tilde{\bm{x}}_{t}^{(2)}=[\bm{x}_{2t-1}^{\top},\bm{x}_{2t}^{\top}]^{\top}$.
From (3), we have that

$\displaystyle\tilde{\bm{x}}_{t}^{(2)}=\tilde{\bm{A}}\tilde{\bm{x}}_{t-1}^{(2)}+\tilde{\bm{\omega}}_{t}^{(2)},$
(45)

where $\tilde{\bm{A}}=\begin{pmatrix}\bm{A}_{2}&\bm{A}_{1}\\\
\bm{A}_{1}\bm{A}_{2}&\bm{A}_{1}\bm{A}_{1}+\bm{A}_{2}\end{pmatrix}$ and
$\tilde{\bm{\omega}}_{t}^{(2)}=\begin{pmatrix}\bm{\omega}_{2t-1}\\\
\bm{A}_{1}\bm{\omega}_{2t-1}+\bm{\omega}_{2t}\end{pmatrix}$ is Gaussian
distributed. Therefore, $\tilde{\bm{x}}_{t}^{(2)}$ follows a first-order AR
model. Let
$\tilde{\bm{y}}_{t}^{(2)}=[\bm{y}_{2t-1}^{\top},\bm{y}_{2t}^{\top}]^{\top}$.
We have that
$\tilde{\bm{y}}_{t}^{(2)}=\tilde{\bm{x}}_{t}^{(2)}+\tilde{\bm{\nu}}_{t}^{(2)}$,
where $\tilde{\bm{\nu}}_{t}^{(2)}=\begin{pmatrix}\bm{\nu}_{2t-1}\\\
\bm{\nu}_{2t}\end{pmatrix}$.

Following the same steps, for a $q$-th order AR model, let
$\tilde{\bm{x}}_{t}^{(q)}=[\bm{x}_{q(t-1)+1}^{\top},\bm{x}_{q(t-1)+2}^{\top},\cdots,\bm{x}_{qt}^{\top}]^{\top}$
and
$\tilde{\bm{y}}_{t}^{(q)}=[\bm{y}_{q(t-1)+1}^{\top},\bm{y}_{q(t-1)+2}^{\top},\cdots,\bm{y}_{qt}^{\top}]^{\top}$.
We have that

$\displaystyle\tilde{\bm{x}}_{t}^{(q)}$
$\displaystyle=\tilde{\bm{A}}\tilde{\bm{x}}_{(t-1)}^{(q)}+\tilde{\bm{\omega}}_{t}^{(q)},$
$\displaystyle\tilde{\bm{y}}_{t}^{(q)}$
$\displaystyle=\tilde{\bm{x}}_{t}^{(q)}+\tilde{\bm{\nu}}_{t}^{(q)},$ (46)

where $\tilde{\bm{\omega}}_{t}^{(q)}$ and $\tilde{\bm{\nu}}_{t}^{(q)}$ are the
innovative noise and measurement noise respectively for the $q$-th order AR
model. We denote the covariance of $\tilde{\bm{\omega}}_{t}^{(q)}$ by
${\widetilde{\bm{R}}}_{\bm{\omega}}$. Here, we omit the expression of
$\widetilde{\bm{A}}$ since it is cumbersome but it can be computed directly.
It can be easily verified that the covariance matrix of
$\tilde{\bm{\nu}}_{t}^{(q)}$ is $\bm{I}$. Therefore,
$\\{\tilde{\bm{x}}_{t}^{(q)},\tilde{\bm{y}}_{t}^{(q)}\\}_{t=1}^{\infty}$ is a
first-order AR model.

##  Appendix B Proof of Lemma 2

###### Proof.

To show that $\bm{\Sigma}_{t}$ converges, it suffices to show that for any
$\epsilon>0$, there exists an integer $T$ such that for any $s,t>T$,
$\|\bm{\Sigma}_{s}-\bm{\Sigma}_{t}\|<\epsilon$. Without loss of generality, we
assume that $s>t$.

We first note that

$\displaystyle\bm{\Sigma}_{t}$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1},$
$\displaystyle\bm{\Sigma}_{s}$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}.$
(47)

We then have that

$\displaystyle\bm{\Sigma}_{s}-\bm{\Sigma}_{t}$
$\displaystyle=\bm{\Sigma}_{t}(\bm{\Sigma}_{t}^{-1}-\bm{\Sigma}_{s}^{-1})\bm{\Sigma}_{s}$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}-(\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{)}$
$\displaystyle\hskip
14.22636pt\big{(}(\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}(\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})-(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})\big{)}$
$\displaystyle\hskip
14.22636pt(\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}(\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}$
$\displaystyle=(\bm{I}+\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}(\bm{\Sigma}_{s-1}-\bm{\Sigma}_{t-1})\bm{A}^{\top}(\bm{I}+\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}.$
(48)

From the update rule of $\bm{\Sigma}_{t}$ in (B), $\bm{\Sigma}_{t}$ is
positive definite for every $t\geq t_{0}$. Therefore,
$\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega}$ is positive definite.
It then allows the eigendecomposition and can be factorized as
$\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega}=\bm{Q}\bm{\Lambda}\bm{Q}^{-1}$,
where $\bm{Q}$ is the square $K\times K$ matrix whose $i$th column is the
eigenvector $q_{i}$ of $\bm{A}$ and $\bm{\Lambda}$ is the diagonal matrix
whose $i$th diagonal elements are the corresponding eigenvalues $\lambda_{i}$.
We then have that

$\displaystyle(\bm{I}+\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}$
$\displaystyle=\big{(}\bm{Q}(\bm{\Lambda}+\bm{I})\bm{Q}^{-1}\big{)}^{-1}=\bm{Q}(\bm{\Lambda}+\bm{I})^{-1}\bm{Q}^{-1}.$
(49)

Since $\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega}$ is positive
definite for any $t>0$, we have that $\lambda_{i}>0$ for $i=1,\cdots,K$.
Therefore, all the diagonal elements of $\bm{\Lambda}+\bm{I}$ are strictly
larger than 1. Therefore, the diagonal elements of
$(\bm{\Lambda}+\bm{I})^{-1}$ are strictly less than 1. From the definition of
operator norm, we have that there exists a $0<\delta<1$ such that
$\|(\bm{I}+\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\|<(1-\delta)$
for any $t>0$.

It then follows that

$\displaystyle\|\bm{\Sigma}_{s}-\bm{\Sigma}_{t}\|$
$\displaystyle=\big{\|}(\bm{I}+\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}(\bm{\Sigma}_{s-1}-\bm{\Sigma}_{t-1})\bm{A}^{\top}(\bm{I}+\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{\|}$
$\displaystyle\leq\|(\bm{I}+\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\|\|\bm{A}\|\|\bm{\Sigma}_{s-1}-\bm{\Sigma}_{t-1}\|\|\bm{A}^{\top}\|\|(\bm{I}+\bm{A}\bm{\Sigma}_{s-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\|$
$\displaystyle<(1-\delta)^{2}\|\bm{\Sigma}_{s-1}-\bm{\Sigma}_{t-1}\|,$ (50)

where the first inequality is from the submultiplicative of operator norm and
the second inequality is due to the facts that
$\|(\bm{I}+\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\|<1-\delta,\|(\bm{I}+\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\|<1-\delta$
and the assumption $\|\bm{A}\|<1$. By applying (B) recursively, we have that

$\displaystyle\|\bm{\Sigma}_{s}-\bm{\Sigma}_{t}\|<(1-\delta)^{2(t-t_{0})}\|\bm{\Sigma}_{s-t_{0}+t}-\bm{\Sigma}_{t_{0}}\|.$
(51)

Moreover, from the update rule of $\bm{\Sigma}_{t}$, it can be easily verified
that $\|\bm{\Sigma}_{s-t_{0}+t}-\bm{\Sigma}_{1}\|$ is bounded. Therefore, for
any $\epsilon>0$, there exists an integer $T$ such that for any $s,t>T$,
$\|\bm{\Sigma}_{s}-\bm{\Sigma}_{t}\|<\epsilon$. This completes the proof. ∎

##  Appendix C Proof of Lemma 3

###### Proof.

Let $\pi(\bm{x}_{t})$ be a Gaussian distribution with mean $\bm{\mu}$ and
covariance matrix $\bm{\Sigma}$. We will show that $\pi(\bm{x}_{t})$ is a
stationary distribution of $\\{\bm{x}_{t}\\}_{t=t_{0}}^{\infty}$ if
$\bm{\mu}=\bm{0}$ and
$\bm{\Sigma}=\sum_{i=0}^{\infty}(\bm{A}^{\top})^{i}\bm{R}_{\omega}\bm{A}^{i}$.
It can be easily proved that $\bm{\Sigma}$ exists since $\|\bm{A}\|<1$, where
$\|\cdot\|$ denotes the operator norm of a matrix.

From the definition of stationary distribution, $\pi(\bm{x}_{t})$ should
satisfy that

$\displaystyle\int\pi(\bm{x}_{t-1})f(\bm{x}_{t}|\bm{x}_{t-1})d\bm{x}_{t-1}=\pi(\bm{x}_{t}).$
(52)

We have that

$\displaystyle\int\pi(\bm{x}_{t-1})f(\bm{x}_{t}|\bm{x}_{t-1})d\bm{x}_{t-1}$
$\displaystyle=\int\pi(\bm{x}_{t-1})\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{R}_{\omega})}}\exp\Big{(}-\frac{1}{2}(\bm{x}_{t}-\bm{A}\bm{x}_{t-1})^{\top}\bm{R}_{\omega}^{-1}(\bm{x}_{t}-\bm{A}\bm{x}_{t-1})\Big{)}d\bm{x}_{t-1}$
$\displaystyle=\int\pi(\bm{x}_{t-1})\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{R}_{\omega})}}\exp\Big{(}-\frac{1}{2}(\bm{A}^{-1}\bm{x}_{t}-\bm{x}_{t-1})^{\top}\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A}(\bm{A}^{-1}\bm{x}_{t}-\bm{x}_{t-1})\Big{)}d\bm{x}_{t-1}$
$\displaystyle=\int\pi(\bm{x}_{t-1})\frac{\sqrt{\det((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1})}}{\sqrt{\det(\bm{R}_{\omega})}\sqrt{(2\pi)^{K}\det((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A}))^{-1}}}\exp$
$\displaystyle\hskip
5.69046pt\Big{(}-\frac{1}{2}(\bm{A}^{-1}\bm{x}_{t}-\bm{x}_{t-1})^{\top}\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A}(\bm{A}^{-1}\bm{x}_{t}-\bm{x}_{t-1})\Big{)}d\bm{x}_{t-1}$
$\displaystyle\overset{(a)}{=}\frac{\sqrt{\det((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1})}}{\sqrt{\det(\bm{R}_{\omega})}\sqrt{(2\pi)^{K}\det((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1}+\bm{\Sigma})}}\exp$
$\displaystyle\hskip
5.69046pt\Big{(}-\frac{1}{2}(\bm{A}^{-1}\bm{x}_{t}-\bm{\mu})^{\top}\big{(}(\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1}+\bm{\Sigma}\big{)}^{-1}(\bm{A}^{-1}\bm{x}_{t}-\bm{\mu})\Big{)}$
$\displaystyle=\frac{\sqrt{\det((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1})}\sqrt{\det\big{(}\bm{A}((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1}+\bm{\Sigma})\bm{A}^{\top}\big{)}}}{\sqrt{\det(\bm{R}_{\omega})}\sqrt{\det\big{(}(\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1}+\bm{\Sigma}\big{)}}}$
$\displaystyle\hskip
5.69046pt\frac{1}{\sqrt{(2\pi)^{K}\det\big{(}\bm{A}((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1}+\bm{\Sigma})\bm{A}^{\top}\big{)}}}\exp\Big{(}-\frac{1}{2}$
$\displaystyle\hskip
5.69046pt(\bm{x}_{t}-\bm{A}\bm{\mu})^{\top}\big{(}\bm{A}((\bm{A}^{\top}\bm{R}_{\omega}^{-1}\bm{A})^{-1}+\bm{\Sigma})\bm{A}^{\top}\big{)}^{-1}(\bm{x}_{t}-\bm{A}\bm{\mu})\Big{)},$
(53)

where $(a)$ is from the convolution of two Gaussian functions.

Let
$\bm{\Sigma}=\sum_{i=0}^{\infty}(\bm{A}^{\top})^{i}\bm{R}_{\omega}\bm{A}^{i}$.
We have that $\bm{\Sigma}=\bm{R}_{\omega}+\bm{A}\bm{\Sigma}\bm{A}^{\top}$. If
we choose $\bm{\mu}=0$ and
$\bm{\Sigma}=\sum_{i=0}^{\infty}(\bm{A}^{\top})^{i}\bm{R}_{\omega}\bm{A}^{i}$,
it can then be easily verified that (52) holds. Therefore, $\pi(\bm{x}_{t})$
is a stationary distribution of $\\{\bm{x}_{t}\\}_{t=t_{0}}^{\infty}$. It then
follows that $\\{\bm{x}_{t},\bm{y}_{t}\\}_{t=t_{0}}^{\infty}$ has a stationary
distribution $\pi(\bm{x}_{t})g(\bm{y}_{t}|\bm{x}_{t})$.

Let $P^{1}\big{(}\\{\bm{x},\bm{y}\\},E\big{)}$ denote the probability of
reaching a measurable set $E$ from state $\\{\bm{x},\bm{y}\\}$ in one step. We
have that for any $\\{\bm{x},\bm{y}\\}$ and $E\in\mathbb{R}^{2K}$ such that
$\pi(E)>0$,

$\displaystyle P^{1}\big{(}\\{\bm{x},\bm{y}\\},E\big{)}$
$\displaystyle=\int_{E}f\big{(}\bm{x}^{\prime}|\bm{x}\big{)}g(\bm{y}^{\prime}|\bm{x}^{\prime})d\bm{x}^{\prime}\bm{y}^{\prime}>0,$
(54)

where the inequality is due to the fact that Gaussian density functions are
positive. From the definition of irreducible [53] Markov chain, we have that
$\\{\bm{x}_{t},\bm{y}_{t}\\}_{t=t_{0}}^{\infty}$ is $\pi$-irreducible. ∎

##  Appendix D Proof of Lemma 4

###### Proof.

To prove Lemma 4, we will first show that
$\\{\bm{x}_{t},\bm{y}_{t},\bm{\mu}_{t}^{*}\\}_{t=t_{0}}^{\infty}$ is a Markov
chain. We will then show that the stationary distribution for
$\bm{\mu}_{t}^{*}$ exists by finding the limiting distribution of
$\bm{\mu}_{t}^{*}$. We then construct a new Markov chain using
$\bm{y}_{t},\bm{\mu}_{t}^{*}$ and show that its stationary distribution exists
and is unique.

Step 1. We first show that
$\\{\bm{x}_{t},\bm{y}_{t},\bm{\mu}_{t}^{*}\\}_{t=t_{0}}^{\infty}$ is a Markov
chain. From the update rule of $\bm{\mu}^{*}_{t}$, we have that for $t\geq
t_{0}$,

$\displaystyle
p_{t_{0}}(\bm{\mu}_{t}^{*}|\bm{x}_{t_{0}},\cdots,\bm{x}_{t},\bm{y}_{t_{0}},\cdots,\bm{y}_{t},\bm{\mu}_{t_{0}}^{*},\cdots,\bm{\mu}_{t-1}^{*})=p_{t_{0}}(\bm{\mu}_{t}^{*}|\bm{\mu}_{t-1}^{*},\bm{y}_{t}).$
(55)

We then have that

$\displaystyle
p_{t_{0}}(\bm{x}_{t},\bm{y}_{t},\bm{\mu}_{t}^{*}|\bm{x}_{t_{0}},\cdots,\bm{x}_{t-1},\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}_{t_{0}}^{*},\cdots,\bm{\mu}_{t-1}^{*})$
$\displaystyle=p_{t_{0}}(\bm{x}_{t}|\bm{x}_{t-1})p_{t_{0}}(\bm{y}_{t}|\bm{x}_{t})p_{t_{0}}(\bm{\mu}_{t}^{*}|\bm{\mu}_{t-1}^{*},\bm{y}_{t})$
$\displaystyle=p_{t_{0}}(\bm{x}_{t},\bm{y}_{t},\bm{\mu}_{t}^{*}|\bm{x}_{t-1},\bm{y}_{t-1},\bm{\mu}_{t-1}^{*}).$
(56)

Therefore, $\\{\bm{x}_{t},\bm{y}_{t},\bm{\mu}_{t}^{*}\\}_{t={t_{0}}}^{\infty}$
is a Markov chain.

Step 2. We then show that the stationary distribution for $\bm{\mu}_{t}^{*}$
exists by finding the limiting distribution of $\bm{\mu}_{t}^{*}$. Let the
initial state be $\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$.
From (III-B), we have that

$\displaystyle\bm{\mu}^{*}_{t}$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{t-{t_{0}}}\bm{\mu}_{t_{0}}^{*}+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{t-{t_{0}}-1}\bm{\Sigma}^{*}\bm{y}_{t_{0}+1}+\cdots$
$\displaystyle\hskip
14.22636pt+(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\Sigma}^{*}\bm{y}_{t-1}+\bm{\Sigma}^{*}\bm{y}_{t}.$
(57)

Note that $\bm{\mu}^{*}_{t}$ is the sum of Gaussian random variables. If
$\bm{y}_{t_{0}},\cdots,\bm{y}_{t}$ are jointly Gaussian distributed, then
$\bm{\mu}^{*}_{t}$ is a Gaussian random variable. Since the Gaussian density
function is continuous in its mean and covariance, it then suffices to show
that the limiting mean and limiting covariance matrix of $\bm{\mu}^{*}_{t}$
exist and are independent of the initial state.

Given the initial state
$\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$, we have that

$\displaystyle
p_{t_{0}}(\bm{y}_{t_{0}+1}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*})$
$\displaystyle=\int
p_{t_{0}}(\bm{x}_{t_{0}+1}|\bm{x}_{t_{0}})p_{t_{0}}(\bm{y}_{t_{0}+1}|\bm{x}_{t_{0}+1})d\bm{x}_{t_{0}+1}$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{R}_{\omega})}}\exp\Big{(}-\frac{1}{2}(\bm{x}_{t_{0}+1}-\bm{A}\bm{x}_{t_{0}})^{\top}\bm{R}_{\omega}^{-1}(\bm{x}_{t_{0}+1}-\bm{A}\bm{x}_{t_{0}})\Big{)}$
$\displaystyle\hskip
5.69046pt\frac{1}{\sqrt{(2\pi)^{K}}}\exp(-\frac{1}{2}(\bm{y}_{t_{0}+1}-\bm{x}_{t_{0}+1})^{\top}(\bm{y}_{t_{0}+1}-\bm{x}_{t_{0}+1}))d\bm{x}_{t_{0}+1}$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{R}_{\omega}+\bm{I})}}\exp\Big{(}-\frac{1}{2}(\bm{y}_{t_{0}+1}-\bm{A}\bm{x}_{t_{0}})^{\top}(\bm{I}+\bm{R}_{\omega})^{-1}(\bm{y}_{t_{0}+1}-\bm{A}\bm{x}_{t_{0}})\Big{)}.$
(58)

Therefore, conditioning on
$\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$, $\bm{y}_{t_{0}+1}$
is Gaussian distributed with mean $\bm{A}\bm{x}_{t_{0}}$ and covariance matrix
$\bm{I}+\bm{R}_{\omega}$. Similarly, we can show that for any $t>t_{0}$,
conditioning on $\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$,
$\bm{y}_{t}\sim\mathcal{N}(\bm{A}^{(t-{t_{0}})}\bm{x}_{t_{0}},\bm{I}+\bm{R}_{\omega}+\bm{A}^{\top}\bm{R}_{\omega}\bm{A}+\cdots+(\bm{A}^{\top})^{(t-t_{0}-1)}\bm{R}_{\omega}\bm{A}^{t-t_{0}-1})$.
Moreover, since $\bm{x}_{t}$ is Gaussian distributed and
$\bm{x}_{t+1}|\bm{x}_{t}\sim\mathcal{N}(\bm{A}\bm{x}_{t},\bm{R}_{\omega})$, we
have that $\bm{x}_{t},\bm{x}_{t+1}$ are jointly Gaussian distributed.
Similarly, since $\bm{x}_{t},\bm{x}_{t+1}$ are jointly Gaussian and
$\bm{y}_{t}|(\bm{x}_{t},\bm{x}_{t+1})\sim\mathcal{N}(\bm{x}_{t},\bm{I})$, we
have that $\bm{x}_{t},\bm{y}_{t},\bm{x}_{t+1}$ are jointly Gaussian
distributed. Following the same idea, we can show that
$\bm{x}_{t_{0}},\bm{x}_{t_{0}+1},\cdots,\bm{x}_{t},\bm{y}_{t_{0}},\bm{y}_{t_{0}+1},\cdots,\bm{y}_{t}$
are jointly Gaussian distributed and thus
$\bm{y}_{t_{0}},\bm{y}_{t_{0}+1},\cdots,\bm{y}_{t}$ are jointly Gaussian
distributed. Therefore, $\bm{\mu}_{t}^{*}$ follows a Gaussian distribution.

To find the limiting distribution of $\bm{\mu}_{t}^{*}$, we first consider the
limiting expectation of $\bm{\mu}^{*}_{t}$. We have that

$\displaystyle\mathbb{E}_{t_{0}}[\bm{\mu}^{*}_{t}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{t-{t_{0}}}\bm{\mu}_{t_{0}}^{*}+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{t-{t_{0}}-1}\bm{\Sigma}^{*}\bm{A}\bm{x}_{t_{0}}+\cdots$
$\displaystyle\hskip
5.69046pt+(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\Sigma}^{*}\bm{A}^{t-{t_{0}}-1}\bm{x}_{t_{0}}+\bm{\Sigma}^{*}\bm{A}^{t-{t_{0}}}\bm{x}_{t_{0}}.$

Since there exists $\delta>0$ such that
$\|(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\|<1-\delta$
and $\|\bm{A}\|<1-\delta$, we have that

$\displaystyle\big{\|}\mathbb{E}_{t_{0}}[\bm{\mu}^{*}_{t}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{\|}$
$\displaystyle\overset{(a)}{\leq}\big{\|}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{t-{t_{0}}}\bm{\mu}_{t_{0}}^{*}\big{\|}+\big{\|}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{t-{t_{0}}-1}\bm{\Sigma}^{*}\bm{A}\bm{x}_{t_{0}}\big{\|}+\cdots$
$\displaystyle\hskip
5.69046pt+\big{\|}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\Sigma}^{*}\bm{A}^{t-{t_{0}}-1}\bm{x}_{t_{0}}\big{\|}+\big{\|}\bm{\Sigma}^{*}\bm{A}^{t-{t_{0}}}\bm{x}_{t_{0}}\big{\|}$
$\displaystyle\overset{(b)}{\leq}\big{\|}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{\|}^{t-{t_{0}}}\big{\|}\bm{\mu}_{t_{0}}^{*}\big{\|}+\big{\|}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{\|}^{t-{t_{0}}-1}\big{\|}\bm{\Sigma}^{*}\big{\|}\big{\|}\bm{A}\big{\|}\big{\|}\bm{x}_{t_{0}}\big{\|}+\cdots$
$\displaystyle\hskip
5.69046pt+\big{\|}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{\|}\big{\|}\bm{\Sigma}^{*}\big{\|}\big{\|}\bm{A}\big{\|}^{t-{t_{0}}-1}\big{\|}\bm{x}_{t_{0}}\big{\|}+\big{\|}\bm{\Sigma}^{*}\big{\|}\big{\|}\bm{A}\big{\|}^{t-{t_{0}}}\big{\|}\bm{x}_{t_{0}}\big{\|}$
$\displaystyle\leq(1-\delta)^{t-{t_{0}}}\big{\|}\bm{\mu}_{t_{0}}^{*}\big{\|}+(t-{t_{0}})(1-\delta)^{t-{t_{0}}}\|\bm{\Sigma}^{*}\|\|\bm{x}_{t_{0}}\|,$

where $(a)$ is from the triangle inequality of norms and $(b)$ is due to the
submultiplicative of the operator norm. We then have that
$\lim_{t\rightarrow\infty}\big{\|}\mathbb{E}_{t_{0}}[\bm{\mu}^{*}_{t}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{\|}=\bm{0}$.
Therefore, we have that
$\lim_{t\rightarrow\infty}\mathbb{E}_{t_{0}}[\bm{\mu}^{*}_{t}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]$
exists and is independent of the initial state
$\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$.

We then consider the covariance matrix of $\bm{\mu}^{*}_{t}$. Given the
initial state $\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$, let

$\displaystyle\bm{\Sigma}_{\bm{\mu}_{t}^{*}}=\mathbb{E}_{t_{0}}\Big{[}\big{(}\bm{\mu}_{t}^{*}-\mathbb{E}_{t_{0}}[\bm{\mu}_{t}^{*}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{)}\big{(}\bm{\mu}_{t}^{*}-\mathbb{E}_{t_{0}}[\bm{\mu}_{t}^{*}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{)}^{\top}\big{|}\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\Big{]},$
$\displaystyle\bm{\Sigma}_{ij}=\mathbb{E}_{t_{0}}\Big{[}\big{(}\bm{y}_{i}-\mathbb{E}_{t_{0}}[\bm{y}_{i}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{)}\big{(}\bm{y}_{j}-\mathbb{E}_{t_{0}}[\bm{y}_{j}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{)}^{\top}\big{|}\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\Big{]}.$
(59)

From (D) and (D), we have that

$\displaystyle\bm{\Sigma}_{\bm{\mu}_{t}^{*}}=\sum_{i={t_{0}}}^{t-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{i-t_{0}}\bm{\Sigma}^{*}\bm{\Sigma}_{(t-i+t_{0})(t-i+t_{0})}\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{i-t_{0}}\bm{\Sigma}^{*}\Big{)}^{\top}$
$\displaystyle+2\sum_{i={t_{0}}}^{t-1}\sum_{j={t_{0}}}^{i-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{i-t_{0}}\bm{\Sigma}^{*}\bm{\Sigma}_{(t-i+t_{0})(t-j+t_{0})}\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{j-t_{0}}\bm{\Sigma}^{*}\Big{)}^{\top}.$
(60)

We first show that $\|\bm{\Sigma}_{ij}\|$ is bounded for any $i,j>t_{0}$. When
$i=j$, we have that

$\displaystyle\|\bm{\Sigma}_{ii}\|$
$\displaystyle=\|\bm{I}+\bm{R}_{\omega}+\bm{A}^{\top}\bm{R}_{\omega}\bm{A}+\cdots+(\bm{A}^{\top})^{(i-{t_{0}}-1)}\bm{R}_{\omega}\bm{A}^{i-{t_{0}}-1}\|$
$\displaystyle\leq\|\bm{I}\|+\|\bm{R}_{\omega}\|+\|\bm{A}^{\top}\|\|\bm{R}_{\omega}\|\|\bm{A}\|+\cdots+\|\bm{A}^{\top}\|^{i-{t_{0}}-1}\|\bm{R}_{\omega}\|\|\bm{A}\|^{i-{t_{0}}-1}$
$\displaystyle=\|\bm{I}\|+\big{(}1+\|\bm{A}\|^{2}+\cdots+\|\bm{A}\|^{2(i-{t_{0}})}\big{)}\|\bm{R}_{\omega}\|$
$\displaystyle=\|\bm{I}\|+\frac{1-\|\bm{A}\|^{2(i-{t_{0}})}}{1-\|\bm{A}\|^{2}}\|\bm{R}_{\omega}\|$
$\displaystyle\leq\|\bm{I}\|+\frac{1}{1-\|\bm{A}\|^{2}}\|\bm{R}_{\omega}\|,$
(61)

where the last inequality is due to the fact that there exists $0<\delta<1$
such that $\|\bm{A}\|<1-\delta$. When $i\neq j$, we have that

$\displaystyle\|\bm{\Sigma}_{ij}\|$
$\displaystyle=\Big{\|}\mathbb{E}_{t_{0}}\Big{[}\big{(}\bm{y}_{i}-\mathbb{E}_{t_{0}}[\bm{y}_{i}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{)}\big{(}\bm{y}_{j}-\mathbb{E}_{t_{0}}[\bm{y}_{j}|\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}]\big{)}^{\top}\big{|}\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\Big{]}\Big{\|}$
$\displaystyle\leq\frac{1}{2}\big{(}\|\bm{\Sigma}_{ii}\|+\|\bm{\Sigma}_{jj}\|\big{)}$
$\displaystyle\leq\|\bm{I}\|+\frac{1}{1-\|\bm{A}\|^{2}}\|\bm{R}_{\omega}\|,$
(62)

where the first inequality is due to the fact that for any two vectors
$\bm{u},\bm{v}\in\mathbb{R}^{K}$, $(\bm{u}-\bm{v})(\bm{u}-\bm{v})^{\top}$ is
positive semi-definite, and thus
$\bm{u}\bm{u}^{\top}+\bm{v}\bm{v}^{\top}-2\bm{u}\bm{v}^{\top}\succeq 0$. Since
the operator norm equals to the largest eigenvalue for the positive semi-
definite matrix, we have that
$\|\bm{u}\bm{v}^{\top}\|\leq\frac{1}{2}\|\bm{u}\bm{u}^{\top}+\bm{v}\bm{v}^{\top}\|\leq\frac{1}{2}(\|\bm{u}\bm{u}^{\top}\|+\|\bm{v}\bm{v}^{\top}\|)$.

To show that $\bm{\Sigma}_{\bm{\mu}_{t}^{*}}$ converges, it suffices to show
that for any $\epsilon>0$, there exists an integer $T\geq t_{0}$ such that for
any $s,t>T$,
$\|\bm{\Sigma}_{\bm{\mu}_{s}^{*}}-\bm{\Sigma}_{\bm{\mu}_{t}^{*}}\|<\epsilon$.
Without loss of generality, we assume that $s>t$. We then have that

$\displaystyle\|\bm{\Sigma}_{\bm{\mu}_{s}^{*}}-\bm{\Sigma}_{\bm{\mu}_{t}^{*}}\|$
$\displaystyle\leq\Big{\|}\sum_{i=t}^{s-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{i-t_{0}}\bm{\Sigma}^{*}\bm{\Sigma}_{(s-i+t_{0})(s-i+t_{0})}\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{i-t_{0}}\bm{\Sigma}^{*}\Big{)}^{\top}\Big{\|}$
$\displaystyle\hskip
2.84544pt+2\Big{\|}\sum_{i=t}^{s-1}\sum_{j=t_{0}}^{i-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{i-t_{0}}\bm{\Sigma}^{*}$
$\displaystyle\hskip
14.22636pt\bm{\Sigma}_{(s-i+t_{0})(s-j+t_{0})}\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\big{)}^{j-t_{0}}\bm{\Sigma}^{*}\Big{)}^{\top}\Big{\|},$

where the inequality is from the triangle inequality of operator norm. To
simplify the notation, let
$b=\|(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\|<1-\delta,\|\widehat{\bm{\Sigma}}\|=\|\bm{I}\|+\frac{1}{1-\|\bm{A}\|^{2}}\|\bm{R}_{\omega}\|$.
It then follows that

$\displaystyle\|\bm{\Sigma}_{\bm{\mu}_{s}^{*}}-\bm{\Sigma}_{\bm{\mu}_{t}^{*}}\|$
$\displaystyle\leq\sum_{i=T+1}^{s-1}b^{2(i-t_{0})}\|\bm{\Sigma}^{*}\|^{2}\|\bm{\Sigma}_{(s-i+t_{0})(s-i+t_{0})}\|+2\sum_{i=T+1}^{s-1}\sum_{j=t_{0}}^{i-1}b^{(i+j-2t_{0})}\|\bm{\Sigma}^{*}\|^{2}\|\bm{\Sigma}_{(s-i+t_{0})(s-j+t_{0})}\|$
$\displaystyle\leq\sum_{i=T+1}^{s-1}b^{2(i-t_{0})}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|+2\sum_{i=T+1}^{s-1}\sum_{j=t_{0}}^{i-1}b^{(i+j-2t_{0})}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|$
$\displaystyle=\frac{b^{2(T+1-t_{0})}(1-b^{2(s-T-1)})}{1-b^{2}}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|+2\sum_{i=T+1}^{s-1}\frac{b^{(i-t_{0})}(1-b^{(i-t_{0})})}{1-b}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|$
$\displaystyle\leq\frac{b^{2(T+1-t_{0})}}{1-b^{2}}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|+2\frac{b^{(T+1-t_{0})}(1-b^{(s-T-1)})}{(1-b)^{2}}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|$
$\displaystyle\leq\frac{b^{2(T+1-t_{0})}}{1-b^{2}}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|+2\frac{b^{(T+1-t_{0})}}{(1-b)^{2}}\|\bm{\Sigma}^{*}\|^{2}\|\widehat{\bm{\Sigma}}\|.$
(63)

Since there exists $0<\delta<1$ such that $b<1-\delta$, for any $\epsilon>0$,
there exists an integer $T$ such that for any $s,t>T$,
$\|\bm{\Sigma}_{\bm{\mu}_{s}^{*}}-\bm{\Sigma}_{\bm{\mu}_{t}^{*}}\|<\epsilon$.
Therefore, $\bm{\Sigma}_{\mu^{*}_{t}}$ converges as $t\rightarrow\infty$.
Therefore, given the initial state
$\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$, the limiting
covariance matrix of $\bm{\mu}^{*}_{t}$ exists. Moreover, from the definition
of $\bm{\Sigma}_{\bm{\mu}^{*}_{t}}$ in (D), $\bm{\Sigma}_{\bm{\mu}^{*}_{t}}$
is independent of the initial state
$\\{\bm{x}_{t_{0}},\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}^{*}\\}$. Since the limiting
distribution of $\bm{\mu}^{*}_{t}$ exists and is independent of the initial
state, from the definition of stationary distribution of Markov chain [53],
$\bm{\mu}^{*}_{t}$ has a stationary distribution, which is a Gaussian
distribution.

Step 3. We then show that
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ is a Markov chain and
its stationary distribution exists and is unique. We have that

$\displaystyle
p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}_{t_{0}}^{*},\cdots,\bm{\mu}^{*}_{t-1})$
$\displaystyle=\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t},\bm{\mu}^{*}_{t_{0}})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}^{*}_{t_{0}})}\cdot\frac{p_{t_{0}}(\bm{\mu}^{*}_{t_{0}+1},\cdots,\bm{\mu}^{*}_{t-1}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t},\bm{\mu}^{*}_{t_{0}})}{p_{t_{0}}(\bm{\mu}^{*}_{t_{0}+1},\cdots,\bm{\mu}^{*}_{t-1}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}^{*}_{t_{0}})}$
$\displaystyle\overset{(a)}{=}\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t},\bm{\mu}^{*}_{t_{0}})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}^{*}_{t_{0}})}$
$\displaystyle=p_{t_{0}}(\bm{y}_{t}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}^{*}_{t_{0}})$
$\displaystyle\overset{(b)}{=}\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})}}\frac{1}{\sqrt{\det((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I})}}$
$\displaystyle\hskip
14.22636pt\exp\bigg{(}-\frac{1}{2}\Big{(}(\bm{A}\bm{\mu}^{*}_{t-1})^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}^{*}_{t-1})$
$\displaystyle\hskip
14.22636pt+\bm{y}_{t}^{\top}\bm{y}_{t}-\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}^{*}_{t-1})+\bm{y}_{t}\big{)}^{\top}$
$\displaystyle\hskip
14.22636pt\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}^{*}_{t-1})+\bm{y}_{t}\big{)}\Big{)}\bigg{)},$
(64)

which depends only on $\bm{\mu}^{*}_{t-1}$, and $(a)$ is due to the fact that
conditioning on $(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})$,
$(\bm{\mu}^{*}_{t_{0}+1},\cdots,\bm{\mu}^{*}_{t-1})$ is independent of
$\bm{y}_{t}$ and $(b)$ is from the update rule of the forward variable in
Lemma 1. We then have that

$\displaystyle
p_{t_{0}}(\bm{y}_{t},\bm{\mu}_{t}^{*}|\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1},\bm{\mu}_{t_{0}}^{*},\cdots,\bm{\mu}_{t-1}^{*})$
$\displaystyle=p_{t_{0}}(\bm{y}_{t}|\bm{\mu}_{t-1}^{*})p_{t_{0}}(\bm{\mu}_{t}^{*}|\bm{\mu}_{t-1}^{*},\bm{y}_{t})$
$\displaystyle=p_{t_{0}}(\bm{y}_{t},\bm{\mu}_{t}^{*}|\bm{y}_{t-1},\bm{\mu}_{t-1}^{*}).$
(65)

Therefore, $\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ is a
Markov chain.

Since the stationary distribution of $\bm{\mu}^{*}_{t}$ exists, we have that

$\displaystyle\int\big{(}\pi(\bm{\mu}_{t-1}^{*})p_{t_{0}}(\bm{y}_{t},\bm{\mu}^{*}_{t}|\bm{\mu}^{*}_{t-1})d\bm{\mu}^{*}_{t-1}\big{)}p_{t_{0}}(\bm{y}_{t+1}|\bm{\mu}_{t}^{*})p_{t_{0}}(\bm{\mu}_{t+1}^{*}|\bm{\mu}_{t},\bm{y}_{t+1})d\bm{\mu}^{*}_{t}\bm{y}_{t}$
$\displaystyle=\int\pi(\bm{\mu}^{*}_{t})p_{t_{0}}(\bm{y}_{t+1}|\bm{\mu}_{t}^{*})p_{t_{0}}(\bm{\mu}_{t+1}^{*}|\bm{\mu}_{t},\bm{y}_{t+1})d\bm{\mu}^{*}_{t}$
$\displaystyle=\int\pi(\bm{\mu}^{*}_{t})p_{t_{0}}(\bm{y}_{t+1},\bm{\mu}^{*}_{t+1}|\bm{\mu}_{t}^{*})d\bm{\mu}^{*}_{t}.$
(66)

Therefore,
$\int\pi(\bm{\mu}_{t-1}^{*})p_{t_{0}}(\bm{y}_{t},\bm{\mu}^{*}_{t}|\bm{\mu}^{*}_{t-1})d\bm{\mu}^{*}_{t-1}$
is a stationary distribution of
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$. Following the same
techniques as in the proof of Lemma 3, we have that
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ is $\pi$-irreducible.
Therefore, the stationary distribution of
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ is unique. ∎

##  Appendix E Proof of Theorem 1

###### Proof.

Step 1. We first consider a special case where
$\bm{\Sigma}_{t_{0}}=\bm{\Sigma}^{*}$. We have that

$\displaystyle\log p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})$
$\displaystyle=\log\frac{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}+\log\frac{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-2})}+\cdots+\log\frac{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\bm{y}_{t_{0}+1})}{p_{t_{0}}^{*}(\bm{y}_{t_{0}})}+\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}})$ (67)

and

$\displaystyle\log\frac{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}$
$\displaystyle=\log\frac{\int\alpha_{t}(\bm{x}_{t})d\bm{x}_{t}}{\int\alpha_{t-1}(\bm{x}_{t-1})d\bm{x}_{t-1}}$
$\displaystyle=\log\Big{(}\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})}}\frac{1}{\sqrt{\det((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I})}}\Big{)}$
$\displaystyle\hskip
5.69046pt-\frac{1}{2}\Big{(}(\bm{A}\bm{\mu}_{t-1}^{*})^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1}^{*})$
$\displaystyle\hskip
5.69046pt+\bm{y}_{t}^{\top}\bm{y}_{t}-\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1}^{*})+\bm{y}_{t}\big{)}^{\top}$
$\displaystyle\hskip
5.69046pt\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1}^{*})+\bm{y}_{t}\big{)}\Big{)},$
(68)

where the second equality is from the updated rule of
$\alpha_{t}(\bm{x}_{t})$.

Since
$\bm{\mu}_{t}^{*}=(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}\bm{A}\bm{\mu}_{t-1}^{*}+\bm{\Sigma}^{*}\bm{y}_{t}$,
we replace $\bm{A}\bm{\mu}_{t-1}^{*}$ by
$(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}_{t}^{*}-\bm{\Sigma}^{*}\bm{y}_{t})$.
Therefore,
$\log\frac{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}$
can be equivalently written as a function of only
$\bm{\mu}^{*}_{t},\bm{y}_{t}$. Let
$h(\bm{\mu}^{*}_{t},\bm{y}_{t})=\log\frac{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}$.
We have the following explicit expression for
$h(\bm{\mu}^{*}_{t},\bm{y}_{t})$,

$\displaystyle h(\bm{\mu}^{*}_{t},\bm{y}_{t})$
$\displaystyle=\log\Big{(}\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})}}\frac{1}{\sqrt{\det((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I})}}\Big{)}$
$\displaystyle\hskip
5.69046pt-\frac{1}{2}\Big{(}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}_{t}^{*}-\bm{\Sigma}^{*}\bm{y}_{t}))^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}_{t}^{*}-\bm{\Sigma}^{*}\bm{y}_{t}))$
$\displaystyle\hskip
5.69046pt+\bm{y}_{t}^{\top}\bm{y}_{t}-\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}_{t}^{*}-\bm{\Sigma}^{*}\bm{y}_{t}))+\bm{y}_{t}\big{)}^{\top}$
$\displaystyle\hskip
5.69046pt\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}((\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})(\bm{\mu}_{t}^{*}-\bm{\Sigma}^{*}\bm{y}_{t}))+\bm{y}_{t}\big{)}\Big{)}$
(69)

Since $\pi(\bm{\mu}^{*}_{t}),\pi(\bm{y}_{t})$ are Gaussian distributions and
$h(\bm{\mu},\bm{y})$ is a quadratic function of $\bm{\mu},\bm{y}$, it can be
easily verified that $\mathbb{E}_{\pi}[|h(\bm{\mu},\bm{y})|]<\infty$. We
further note that $\mathbb{R}^{2K}$ is a state space with countably generated
$\sigma$-algebra [52] and
$\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$ is $\pi$-irreducible
and aperiodic. Therefore, from the ergodic theorem of Markov chain [52, 53],
we have that under $\mathbb{P}_{t_{0}}$,

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]$
(70)

almost surely. Moreover, the convergence result does not depend on the initial
state of $\\{\bm{y}_{t},\bm{\mu}^{*}_{t}\\}_{t={t_{0}}}^{\infty}$.

Step 2. We then show that $\lim_{t\rightarrow\infty}\frac{1}{t}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]$
almost surely under $\mathbb{P}_{t_{0}}$ by showing that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\big{(}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})-\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})\big{)}=0$ (71)

almost surely under $\mathbb{P}_{t_{0}}$.

We have that

$\displaystyle\log p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})$
$\displaystyle=\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}+\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-2})}+\cdots+\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\bm{y}_{t_{0}+1})}{p_{t_{0}}(\bm{y}_{t_{0}})}+p_{t_{0}}(\bm{y}_{t_{0}})$
(72)

and

$\displaystyle\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t})}{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t-1})}$
$\displaystyle=\log\Big{(}\frac{1}{\sqrt{(2\pi)^{K}\det(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})}}\frac{1}{\sqrt{\det((\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I})}}\Big{)}$
$\displaystyle\hskip
5.69046pt-\frac{1}{2}\Big{(}(\bm{A}\bm{\mu}_{t-1})^{\top}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1})$
$\displaystyle\hskip
5.69046pt+\bm{y}_{t}^{\top}\bm{y}_{t}-\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1})+\bm{y}_{t}\big{)}^{\top}$
$\displaystyle\hskip
5.69046pt\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{I}\big{)}^{-1}\big{(}(\bm{A}\bm{\Sigma}_{t-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{t-1})+\bm{y}_{t}\big{)}\Big{)}.$
(73)

To show that $\lim_{t\rightarrow\infty}\frac{1}{t}\big{(}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})-\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})\big{)}=0$. It suffices
to show the limit sum of each term in (E) converges. Here we provide the proof
of
$\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\Big{(}(\bm{A}\bm{\mu}_{i})^{\top}(\bm{A}\bm{\Sigma}_{i-1}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i})-(\bm{A}\bm{\mu}_{i}^{*})^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})\Big{)}=0$.
The rest of the terms can be proved using the same techniques. Let

$\displaystyle(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}=(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t}^{1},$
$\displaystyle(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}=(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t}^{2},$
$\displaystyle(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}=(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{\Lambda}_{t}^{3}.$

Since $\bm{\Sigma}_{t}$ converges to $\bm{\Sigma}^{*}$ as $t\rightarrow\infty$
and
$(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1},(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1},(\bm{A}\bm{\Sigma}_{t}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}$
are continuous in $\bm{\Sigma}_{t}$, we have that
$\lim_{t\rightarrow\infty}\|\bm{\Lambda}_{t}^{1}\|=0$,
$\lim_{t\rightarrow\infty}\|\bm{\Lambda}_{t}^{2}\|=0$ and
$\lim_{t\rightarrow\infty}\|\bm{\Lambda}_{t}^{3}\|=0$. Therefore,
$\|\bm{\Lambda}_{t}^{1}\|,\|\bm{\Lambda}_{t}^{2}\|,\|\bm{\Lambda}_{t}^{3}\|$
are uniformly upper bounded for any $t$. Denote by
$\|\bm{\Lambda}^{1}\|,\|\bm{\Lambda}^{2}\|,\|\bm{\Lambda}^{3}\|$ the upper
bound of
$\|\bm{\Lambda}_{t}^{1}\|,\|\bm{\Lambda}_{t}^{2}\|,\|\bm{\Lambda}_{t}^{3}\|$
respectively. Given the initial $\bm{\mu}_{t_{0}}$, from (III-A), we have the
following representation for $\bm{\mu}_{t}$:

$\displaystyle\bm{\mu}_{t_{0}+1}=\bm{\Lambda}_{t_{0}+1}^{1}\bm{A}\bm{\mu}_{t_{0}}+\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1}+\bm{\mu}_{t_{0}+1}^{*},$
$\displaystyle\bm{\mu}_{t_{0}+2}=\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+1}^{1}\bm{A}\bm{\mu}_{t_{0}}+\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1}+\bm{\mu}_{t_{0}+1}^{*})$
$\displaystyle\hskip
5.69046pt+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{2}\big{)}\bm{y}_{t_{0}+2}$
$\displaystyle=\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+1}^{1}\bm{A}\bm{\mu}_{t_{0}}+\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})+\bm{\Lambda}_{t_{0}+2}^{2}\bm{y}_{t_{0}+2}+\bm{\mu}_{t_{0}+2}^{*}$
$\displaystyle\hskip 85.35826pt\vdots\hskip 85.35826pt\vdots\hskip
85.35826pt\vdots$
$\displaystyle\bm{\mu}_{t_{0}+t-1}=\bm{\Lambda}_{t_{0}+t-1}^{2}\bm{y}_{t_{0}+t-1}+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+t-1}^{1}\big{)}\bm{A}\bm{\Lambda}_{t_{0}+t-2}^{2}\bm{y}_{t_{0}+t-2}$
$\displaystyle\hskip
5.69046pt+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+t-1}^{1}\big{)}\bm{A}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+t-2}^{1}\big{)}\bm{A}$
$\displaystyle\hskip
5.69046pt\bm{\Lambda}_{t_{0}+t-3}^{2}\bm{y}_{t_{0}+t-3}+\cdots+\prod_{i=t_{0}+2}^{t_{0}+t-1}\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{i}^{1}\big{)}\bm{A}\Big{)}(\bm{\Lambda}_{t_{0}+1}^{1}\bm{A}\bm{\mu}_{t_{0}}+\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})$
$\displaystyle\hskip 5.69046pt+\bm{\mu}_{t_{0}+t-1}^{*}.$ (74)

Let $\bm{Y}_{t}=(\bm{y}_{t_{0}},\bm{y}_{t_{0}+1},\cdots,\bm{y}_{t_{0}+t-1})$
and denote $\bm{\mu}_{t}$ as $F_{t}(\bm{Y}_{t})+\bm{\mu}^{*}_{t}$. We then
have that

$\displaystyle\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\Big{(}(\bm{A}\bm{\mu}_{i})^{\top}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{\Lambda}_{i}^{3}\big{)}(\bm{A}\bm{\mu}_{i})-(\bm{A}\bm{\mu}_{i}^{*})^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})\Big{)}$
$\displaystyle=\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\Big{(}\big{(}\bm{A}F_{i}(\bm{Y}_{i})+\bm{A}\bm{\mu}^{*}_{i}\big{)}^{\top}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{\Lambda}_{i}^{3}\big{)}$
$\displaystyle\hskip
5.69046pt\big{(}\bm{A}F_{i}(\bm{Y}_{i})+\bm{A}\bm{\mu}^{*}_{i}\big{)}-(\bm{A}\bm{\mu}_{i}^{*})^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})\Big{)}$
$\displaystyle=\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\Big{(}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}$
$\displaystyle\hskip
5.69046pt+\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}\bm{\Lambda}_{i}^{3}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}+\big{(}\bm{A}\bm{\mu}^{*}_{i}\big{)}^{\top}\bm{\Lambda}_{i}^{3}\big{(}\bm{A}\bm{\mu}^{*}_{i}\big{)}$
$\displaystyle\hskip
5.69046pt+2\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})+2\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}\bm{\Lambda}_{i}^{3}(\bm{A}\bm{\mu}_{i}^{*})\Big{)}.$
(75)

In the following, we will show that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}=0$

and

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})=0.$

The convergence of the rest of terms can be proved using the same techniques.
To simplify the notation, we set $\bm{\mu}_{t_{0}}=\bm{0}$. The proof for any
arbitrary $\bm{\mu}_{t_{0}}$ can also be derived similarly. For
$\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}$,
$i=t_{0},\cdots,t$, we have that

$\displaystyle\big{(}\bm{A}F_{t_{0}+1}(\bm{Y}_{t_{0}+1})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{t_{0}+1}(\bm{Y}_{t_{0}+1})\big{)}$
$\displaystyle=(\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})^{\top}\bm{A}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1},$
$\displaystyle\big{(}\bm{A}F_{t_{0}+2}(\bm{Y}_{t_{0}+2})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{t_{0}+2}(\bm{Y}_{t_{0}+2})\big{)}$
$\displaystyle=\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})+\bm{\Lambda}_{t_{0}+2}^{2}\bm{y}_{t_{0}+2}\Big{)}^{\top}\bm{A}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}$
$\displaystyle\hskip
5.69046pt\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})+\bm{\Lambda}_{t_{0}+2}^{2}\bm{y}_{t_{0}+2}\Big{)},$
$\displaystyle\big{(}\bm{A}F_{t_{0}+3}(\bm{Y}_{t_{0}+3})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{t_{0}+3}(\bm{Y}_{t_{0}+3})\big{)}$
$\displaystyle=\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+3}^{1}\big{)}\bm{A}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})$
$\displaystyle\hskip
5.69046pt+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+3}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+2}^{2}\bm{y}_{t_{0}+2})+\bm{\Lambda}_{t_{0}+3}^{2}\bm{y}_{t_{0}+3}\Big{)}^{\top}\bm{A}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}$
$\displaystyle\hskip
5.69046pt\Big{(}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+3}^{1}\big{)}\bm{A}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+2}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+1}^{2}\bm{y}_{t_{0}+1})$
$\displaystyle\hskip
5.69046pt+\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{t_{0}+3}^{1}\big{)}\bm{A}(\bm{\Lambda}_{t_{0}+2}^{2}\bm{y}_{t_{0}+2})+\bm{\Lambda}_{t_{0}+3}^{2}\bm{y}_{t_{0}+3}\Big{)}$
$\displaystyle\hskip 85.35826pt\vdots\hskip 85.35826pt\vdots\hskip
85.35826pt\vdots$ (76)

We note that there exists a $0<\delta<1$ such that for any $i\geq t_{0}$,
$\big{\|}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I})^{-1}+\bm{\Lambda}_{i}^{1}\big{)}\bm{A}\big{\|}=\big{\|}\big{(}\bm{A}\bm{\Sigma}_{i}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I}\big{)}^{-1}\bm{A}\big{\|}\leq\big{\|}\big{(}\bm{A}\bm{\Sigma}_{i}\bm{A}^{\top}+\bm{R}_{\omega}+\bm{I}\big{)}^{-1}\big{\|}\|\bm{A}\|\leq(1-\delta)$.
Let $B=1-\delta$ and
$C=\big{\|}\bm{A}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}\big{\|}$.
We then have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\Big{\|}\Big{(}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}\Big{\|}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{\|}F_{i}(\bm{Y}_{i})\big{\|}^{2}\big{\|}\bm{A}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}\big{\|}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}C\big{\|}F_{i}(\bm{Y}_{i})\big{\|}^{2}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}C\bigg{(}\sum_{j=t_{0}}^{i}\
B^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|\bigg{)}^{2}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}\sum_{k=t_{0}}^{i}CB^{2i-j-k}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|\|\bm{\Lambda}^{2}_{k}\|\|\bm{y}_{k}\|$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}\sum_{k=t_{0}}^{i}CB^{2i-j-k}\big{(}\|\bm{\Lambda}^{2}_{j}\|^{2}\|\bm{y}_{j}\|^{2}+\|\bm{\Lambda}^{2}_{k}\|^{2}\|\bm{y}_{k}\|^{2}\big{)},$
(77)

where the first inequality is due to the submultiplicative of operator norm,
and the third inequality is from the explicit expression of
$F_{i}(\bm{Y}_{i})$ and the triangle inequality and the submultiplicative of
operator norm. We then consider the coefficient of the sum of all terms
containing $\|\bm{\Lambda}^{2}_{t_{0}}\|^{2}\|\bm{y}_{t_{0}}\|^{2}$, denoted
by $Co(\|\bm{\Lambda}^{2}_{t_{0}}\|^{2}\|\bm{y}_{t_{0}}\|^{2})$. We have that

$\displaystyle Co(\|\bm{\Lambda}^{2}_{t_{0}}\|^{2}\|\bm{y}_{t_{0}}\|^{2})$
$\displaystyle=\lim_{t\rightarrow\infty}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}2CB^{2i-t_{0}-j}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\sum_{i=t_{0}}^{t_{0}+t-1}\frac{2CB^{i-t_{0}}}{1-B}$
$\displaystyle\leq\frac{2C}{(1-B)^{2}}.$ (78)

Similarly, we can show that the coefficient of
$\|\bm{\Lambda}^{2}_{t_{0}+1}\|^{2}\|\bm{y}_{t_{0}+1}\|^{2},\cdots,\|\bm{\Lambda}^{2}_{t_{0}+t-1}\|^{2}\|\bm{y}_{t_{0}+t-1}\|^{2}$
are not larger than $\frac{2C}{(1-B)^{2}}$. Since
$\lim_{t\rightarrow\infty}\|\bm{\Lambda}_{t}^{2}\|^{2}=0$, for any
$\epsilon>0$, there exists an integer $T$ such that for any $t>T$,
$\|\bm{\Lambda}_{t}^{2}\|^{2}<\epsilon$. We then have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}2C\frac{\|\bm{\Lambda}^{2}_{t_{0}}\|^{2}\|\bm{y}_{t_{0}}\|^{2}+\cdots+\|\bm{\Lambda}^{2}_{t_{0}+t-1}\|^{2}\|\bm{y}_{t_{0}+t-1}\|^{2}}{(1-B)^{2}}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}2C\bigg{(}\frac{\|\bm{\Lambda}^{2}_{t_{0}}\|^{2}\|\bm{y}_{t_{0}}\|^{2}+\cdots+\|\bm{\Lambda}^{2}_{T}\|^{2}\|\bm{y}_{T}\|^{2})}{(1-B)^{2}}+\frac{\epsilon\big{(}\|\bm{y}_{T+1}\|^{2}+\cdots+\|\bm{y}_{t_{0}+t-1}\|^{2}\big{)}}{(1-B)^{2}}\bigg{)}$
$\displaystyle\leq\frac{2\epsilon
C\mathbb{E}_{\pi}[\|\bm{y}\|^{2}]}{(1-B)^{2}},$ (79)

where the last inequality is from the ergodic theorem of Markov chain [52,
53].

Since $\mathbb{E}_{\pi}[\|\bm{y}\|^{2}]$ is bounded, by letting
$\epsilon\rightarrow 0$, we have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{\|}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}\big{\|}=0.$
(80)

Therefore, we have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}=\bm{0}.$
(81)

We then consider
$\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})$.
We have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{\|}\big{(}\bm{A}F_{i}(\bm{Y}_{i})\big{)}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})\big{\|}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\big{\|}F_{i}(\bm{Y}_{i})\big{\|}\big{\|}\bm{A}^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}\bm{A}\big{\|}\big{\|}\bm{\mu}^{*}_{i}\big{\|}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}C\big{\|}F_{i}(\bm{Y}_{i})\big{\|}\big{\|}\bm{\mu}^{*}_{i}\big{\|}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|\|\bm{\mu}^{*}_{i}\|$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\big{(}\|\bm{y}_{j}\|^{2}+\|\bm{\mu}^{*}_{i}\|^{2}\big{)}.$
(82)

It suffices to show that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|^{2}=0$
(83)

and

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{\mu}^{*}_{i}\|^{2}=0.$
(84)

We have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|^{2}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=t_{0}}^{t_{0}+t-1}\sum_{i=j}^{t_{0}+t-1}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|^{2}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=t_{0}}^{t_{0}+t-1}\frac{C}{1-B}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|^{2}.$
(85)

For any $\epsilon>0$, there exists an integer $T$ such that for any $t>T$,
$\|\bm{\Lambda}^{2}_{t}\|<\epsilon$. It then follows that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=t_{0}}^{t_{0}+t-1}\frac{C}{1-B}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|^{2}$
$\displaystyle\leq\lim_{t\rightarrow\infty}\frac{1}{t}C\Big{(}\frac{\|\bm{\Lambda}^{2}_{t_{0}}\|\|\bm{y}_{t_{0}}\|^{2}+\cdots+\|\bm{\Lambda}^{2}_{T}\|\|\bm{y}_{T}\|^{2}}{1-B}+\frac{\epsilon\|\bm{y}_{T+1}\|^{2}+\cdots+\epsilon\|\bm{y}_{t_{0}+t-1}\|^{2}}{1-B}\Big{)}$
$\displaystyle\leq\frac{\epsilon C\mathbb{E}_{\pi}[\|\bm{y}\|^{2}]}{1-B}.$
(86)

Since $\mathbb{E}_{\pi}[\|\bm{y}\|^{2}]$ is bounded, by letting
$\epsilon\rightarrow 0$, we have that
$\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=t_{0}}^{t_{0}+t-1}\frac{C}{1-B}\|\bm{\Lambda}^{2}_{j}\|\|\bm{y}_{j}\|^{2}=0$.

Consider
$\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{\mu}^{*}_{i}\|^{2}$.
We have that $\lim_{t\rightarrow\infty}B^{t}=0$ and
$\lim_{t\rightarrow\infty}\|\bm{\Lambda}_{t}^{2}\|=0$. Therefore, for any
$\epsilon>0$, there exists $S$ such that for $s>S$,
$\|\bm{\Lambda}_{s-1}^{2}\|<\epsilon$ and $B^{s-1}<\epsilon$. Let $T=2S$. For
any $t>t_{0}+T$, we have that

$\displaystyle\sum_{j=t_{0}}^{t}CB^{t-j}\|\bm{\Lambda}^{2}_{j}\|$
$\displaystyle=C\big{(}\|\bm{\Lambda}^{2}_{t}\|+B\|\bm{\Lambda}^{2}_{t-1}\|+\cdots+B^{t-t_{0}}\|\bm{\Lambda}_{t_{0}}^{2}\|\big{)}$
$\displaystyle\leq
C\big{(}\epsilon+B\epsilon+\cdots+B^{\lfloor\frac{t-t_{0}}{2}\rfloor}\epsilon+B^{\lfloor\frac{t-t_{0}}{2}\rfloor+1}\|\bm{\Lambda}_{t-\lfloor\frac{t-t_{0}}{2}\rfloor-1}^{2}\|+\cdots+B^{t-t_{0}}\|\bm{\Lambda}_{t_{0}}^{2}\|\big{)}$
$\displaystyle\leq
C\Big{(}\frac{\epsilon}{1-B}+\frac{B^{\lfloor\frac{t-t_{0}}{2}\rfloor+1}}{1-B}\|\bm{\Lambda}^{2}\|\Big{)}$
$\displaystyle\leq\frac{C(1+\|\bm{\Lambda}^{2}\|)}{1-B}\epsilon,$ (87)

where for the second inequality, we compute the sum of the first half and the
second half respectively and use the fact that
$\|\bm{\Lambda}_{t}^{2}\|\leq\|\bm{\Lambda}^{2}\|$ for any $t$. Since
$\frac{C(1+\|\bm{\Lambda}^{2}\|)}{1-B}$ is bounded, the coefficient of
$\|\bm{\mu}^{*}_{t}\|$ will converge to zero as $t\rightarrow\infty$.
Therefore, following the same steps as in (E), we have that
$\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\sum_{j=t_{0}}^{i}CB^{i-j}\|\bm{\Lambda}^{2}_{j}\|\|\bm{\mu}^{*}_{i}\|^{2}=0$.
Therefore, we have that

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=t_{0}}^{t_{0}+t-1}\Big{(}(\bm{A}\bm{\mu}_{i})^{\top}\big{(}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}+\bm{\Lambda}_{i}^{3}\big{)}(\bm{A}\bm{\mu}_{i})-(\bm{A}\bm{\mu}_{i}^{*})^{\top}(\bm{A}\bm{\Sigma}^{*}\bm{A}^{\top}+\bm{R}_{\omega})^{-1}(\bm{A}\bm{\mu}_{i}^{*})\Big{)}=0,$
(88)

and thus

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\big{(}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})-\log
p_{t_{0}}^{*}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})\big{)}=0.$ (89)

Step 3. Since the observations are independent before the change point, we
have that $\lim_{t\rightarrow\infty}\frac{1}{t}\log
p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i=1}^{t_{0}+t-1}\log
p_{\infty}(\bm{y}_{i})$. Since $\log p_{\infty}(\bm{y})$ is a quadratic
function of $\bm{y}$, we have that $\mathbb{E}_{\pi}[|\log
p_{\infty}(\bm{y})|]<\infty$. From the ergodic theorem of Markov chain [52,
53], we have that under $\mathbb{P}_{t_{0}}$,

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{i={t_{0}}}^{t_{0}+t-1}\log
p_{\infty}(\bm{y}_{i})=\mathbb{E}_{\pi}[\log p_{\infty}(\bm{y})]$ (90)

almost surely.

It then follows that under $\mathbb{P}_{t_{0}}$,

$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}$
$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{t}\Big{(}\log
p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})-\sum_{i=t_{0}}^{t_{0}+t-1}\log
p_{\infty}(\bm{y}_{i})\Big{)}$
$\displaystyle=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]-\mathbb{E}_{\pi}[\log
p_{\infty}(\bm{y})]$ (91)

almost surely.

Let

$\displaystyle\mathcal{K}=\mathbb{E}_{\pi}[h(\bm{\mu},\bm{y})]-\mathbb{E}_{\pi}[\log
p_{\infty}(\bm{y})].$ (92)

We then have that for any initial state $\bm{y}_{t_{0}},\bm{\mu}_{t_{0}}$,
$\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{p_{t_{0}}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}{p_{\infty}(\bm{y}_{t_{0}},\cdots,\bm{y}_{t_{0}+t-1})}=\mathcal{K}$
under $\mathbb{P}_{t_{0}}$ almost surely.

From (E), we have that given $\bm{\mu}^{*}_{t}$, the distribution of
$\bm{y}_{t}$ under $\mathbb{P}_{t_{0}}$ can be fully specified, denoted by
$p_{t_{0}}^{*}(\bm{y}_{t}|\bm{\mu}^{*}_{t})$. Therefore, we have that
$\pi(\bm{\mu},\bm{y})=\pi(\bm{\mu})p_{t_{0}}^{*}(\bm{y}|\bm{\mu})$. Note that
$h(\bm{\mu},\bm{y})=\log p_{t_{0}}^{*}(\bm{y}|\bm{\mu})$. We then have that

$\displaystyle\mathcal{K}$
$\displaystyle=\mathbb{E}_{\pi}\Big{[}\log\frac{p_{t_{0}}^{*}(\bm{y}|\bm{\mu})}{p_{\infty}(\bm{y})}\Big{]}$
$\displaystyle=\mathbb{E}_{\pi}\Big{[}\log\frac{\pi(\bm{\mu})p_{t_{0}}^{*}(\bm{y}|\bm{\mu})}{\pi(\bm{\mu})p_{\infty}(\bm{y})}\Big{]}$
$\displaystyle=D\big{(}\pi(\bm{\mu})p_{t_{0}}^{*}(\bm{y}|\bm{\mu})\|\pi(\bm{\mu})p_{\infty}(\bm{y})\big{)}$
$\displaystyle\geq 0,$ (93)

where $D(\cdot\|\cdot)$ denotes the KL-Divergence between two distributions
and the equality holds if and only if
$\pi(\bm{\mu})p_{t_{0}}^{*}(\bm{y}|\bm{\mu})\neq\pi(\bm{\mu})p_{\infty}(\bm{y})$.
From (E), we have that
$\pi(\bm{\mu})p_{t_{0}}^{*}(\bm{y}|\bm{\mu})\neq\pi(\bm{\mu})p_{\infty}(\bm{y})$.
Therefore, we have that $\mathcal{K}>0$.

Therefore, for any $\eta>0$,

$\displaystyle\lim_{t\rightarrow\infty}$ $\displaystyle\sup_{t_{0}\geq
1}\text{esssup}\mathbb{P}_{t_{0}}\Big{\\{}\max_{k\leq
t}\sum_{i=t_{0}}^{t_{0}+k-1}\log\frac{p_{t_{0}}(\bm{y}_{i}|\bm{y}_{t_{0}},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i})}\mathcal{K}(1+\eta)t\big{|}\bm{y}_{1},\cdots,\bm{y}_{t_{0}-1}\Big{\\}}=0.$
(94)

Then (20) follows from [24, Theorem 1]. ∎

##  Appendix F Proof of Theorem 3

###### Proof.

Let $\sigma_{0}=0$ and define the stopping times

$\displaystyle\sigma_{m+1}=\inf\Big{\\{}t>\sigma_{m}:\sum^{t}_{i=\sigma_{m}+1}\log\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}\leq
0\Big{\\}},\ \text{for}\ m\geq 0.$ (95)

We have that

$\displaystyle\mathbb{E}_{\infty}\bigg{[}\prod_{i=\sigma_{m}+1}^{t+1}\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}\
\text{for some}\ t>\sigma_{m}\Big{|}\mathcal{F}_{t}\bigg{]}$
$\displaystyle=\prod_{i=\sigma_{m}+1}^{t}\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}\times\mathbb{E}_{\infty}\bigg{[}\frac{p_{1}(\bm{y}_{t+1}|\bm{y}_{1},\cdots,\bm{y}_{t})}{p_{\infty}(\bm{y}_{t+1}|\bm{y}_{1},\cdots,\bm{y}_{t})}\Big{|}\mathcal{F}_{t}\bigg{]}$
$\displaystyle=\prod_{i=\sigma_{m}+1}^{t}\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}\times\int
p_{1}(\bm{y}_{t+1}|\bm{y}_{1},\cdots,\bm{y}_{t})d\bm{y}_{t+1}$
$\displaystyle=\prod_{i=\sigma_{m}+1}^{t}\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}.$
(96)

Therefore,
$\\{\prod_{i=\sigma_{m}+1}^{t}\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})},\mathcal{F}_{t},t>\sigma_{m}\\}$
is a martingale under the pre-change distribution with mean 1. Therefore, from
Doob’s submartingale inequality[58] and the optional sampling theorem[58], we
have that

$\displaystyle\mathbb{P}_{\infty}\Big{\\{}\sum_{i=\sigma_{m}+1}^{t}\log\frac{p_{1}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}{p_{\infty}(\bm{y}_{i}|\bm{y}_{1},\cdots,\bm{y}_{i-1})}\geq