\phdthesis\department

Mathematics \degreeyear2024 \advisorYuliy Baryshnikov \committeeProfessor Vadim Zharnitsky, Chair
Professor Yuliy Baryshnikov, Director of Research
Professor Richard Sowers
Professor Lee DeVille

Cyclicity Analysis of the Ornstein-Uhlenbeck Process

Vivek Kaushik

Abstract

In this thesis, we consider an $N$ -dimensional Ornstein-Uhlenbeck (OU) process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ satisfying the linear stochastic differential equation $d\mathbf{x}(t)=-\mathbf{B}\ \mathbf{x}(t)\ dt+\boldsymbol{\Sigma}\ d\mathbf{w}% (t).$ Here, $\mathbf{B}$ is a fixed $N\times N$ circulant friction matrix whose eigenvalues have positive real parts, $\boldsymbol{\Sigma}$ is a fixed $N\times M$ matrix for some $M\in\mathbb{N},$ and $\left\{\mathbf{w}(t)\right\}_{t\geq 0}$ is the standard $M$ -dimensional Wiener process.

We consider a signal propagation model governed by this OU process. In this model, an underlying signal propagates throughout a network consisting of $N$ linked sensors located in space. For each $t\geq 0,$ we interpret $x_{n}(t),$ the $n$ -th component of the OU process at time $t,$ as the measurement of the propagating effect made by the $n$ -th sensor. The matrix $\mathbf{B}$ represents the sensor network structure: if $\mathbf{B}$ has first row $(b_{1}\ ,\ \dots\ ,\ b_{N}),$ where $b_{1}>0$ and $b_{2}\ ,\ \dots\ ,\ b_{N}\leq 0,$ then the magnitude of $b_{p}$ quantifies how receptive the $n$ -th sensor is to activity within the $(n+p-1)$ -th sensor, where $n+p-1$ is indexed mod $N.$ Finally, the $(m,n)$ -th entry of the matrix $\mathbf{D}=\frac{\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\text{T}}}{2}$ is the covariance of the component noises injected into the $m$ -th and $n$ -th sensors.

For different choices of $\mathbf{B}$ and $\boldsymbol{\Sigma},$ we investigate whether Cyclicity Analysis enables us to recover the structure of network. Roughly speaking, Cyclicity Analysis studies the lead-lag dynamics pertaining to the components of a multivariate signal. We specifically consider an $N\times N$ skew-symmetric matrix $\mathbf{Q},$ known as the lead matrix, in which the sign of its $(m,n)$ -th entry captures the lead-lag relationship between the $m$ -th and $n$ -th component OU processes. We investigate whether the structure of the leading eigenvector of $\mathbf{Q},$ the eigenvector corresponding to the largest eigenvalue of $\mathbf{Q}$ in modulus, reflects the network structure induced by $\mathbf{B}.$

Acknowledgements.

Many people have supported me throughout my Ph.D. journey. I would like to acknowledge them here. Firstly, I would like to thank my advisor Professor Baryshnikov for his guidance. Professor Baryshnikov introduced me to the topic of Cyclicity Analysis and strongly encouraged me to explore its applications on real-world data. Later on, we formulated the topic for this thesis. Throughout my Ph.D. journey, Professor Baryshnikov instilled the meaning of a Ph.D. scholar: a scholar who is able to conduct independent research. Via his meticulous comments and feedback on my dissertation drafts, he taught me that being a good researcher does not just involve doing research but also involves articulating my findings in a clear and concise manner. Next, I would like to thank all the other thesis committee members. During my undergraduate career, Professor Sowers supervised my first data science project, which initiated my curiosity for the field. He also taught me Stochastic Calculus, a course that was essential to my thesis research. Professor DeVille strongly encouraged me to pursue a Ph.D. at the University of Illinois and taught me Dynamical Systems, another course that was essential to my thesis research. Although I did not take courses or do research with him, Professor Zharnitsky and I would communicate frequently especially in Russian. I thank him for kindly lending me a book on the Russian language. Furthermore, I would like to thank Professor Reznick and Professor D’Angelo from the University of Illinois and Professor Ritelli from the University of Bologna. All professors played a significant role during my early research career. I enrolled as an undergraduate freshman in Professor Reznick’s senior level research seminar, during which I learned how to conduct and share research. The work done in that seminar would later turn into joint work with Professor Ritelli during my undergraduate career. Professor Ritelli guided me in publishing my first undergraduate paper in the AMS Quarterly of Applied Mathematics. Lastly, Professor D’Angelo has always been supportive of my research and regularly encouraged me to share my research findings with him. He also shared with me his insights on the math department and on issues related to the broader mathematical community. Moreover, I thank all of my colleagues at eBay. During my Ph.D. journey, I had the opportunity to work at eBay as a data scientist and collaborate with many experts in the field. eBay gave me a unique platform to showcase my skills in solving business related problems. My colleagues were very supportive of me and helped me grow as a data scientist. The work I did at eBay helped shaped my decision to pursue a Ph.D. Finally, I would like to thank my loving family. Admittedly, I faced various challenges pertaining to the completion of the thesis. But my family has always been extremely supportive of me nonetheless.

Chapter 1 Introduction

Throughout this chapter, we fix $N\in\mathbb{N}$ and let $T$ be an interval in $\mathbb{R},$ which we regard as an interval of time. We let $\mathbf{x}:T\rightarrow\mathbb{R}^{N}$ be a continuous function. Note that we can write $\mathbf{x}=(x_{1}\ ,\ \dots\ ,\ x_{N}),$ in which $x_{n}:T\rightarrow\mathbb{R}$ is a real-valued continuous function. We refer to $x_{n}$ as the $n$ -th component of $\mathbf{x}.$

In many practical applications, $\mathbf{x}$ is a signal (wave/pulse) representing some observed phenomenon, in which we interpret $\mathbf{x}(t)$ as a measurement of the observed phenomenon recorded at the fixed time $t\in T.$ For example, if we are observing the behavior of the entire stock market, then we could consider an $N$ -dimensional signal $\mathbf{x}:T\rightarrow\mathbb{R}^{N},$ where $N$ is the total number of companies in the market; here, the component $x_{n}(t)$ could be the stock price for the $n$ -th company at the time $t\in T.$

1.1 Leader-Follower Dynamics Questions

In this section, we pose two very general questions pertaining to the leader-follower dynamics amongst the $N$ components of $\mathbf{x}$ relative to their evolution throughout time. These questions are central to the thesis as a whole.

First, fix two specific indices $1\leq m,n\leq N$ and assume that the corresponding component signals $x_{m}$ and $x_{n}$ evolve very similarly up to a time shift. Does $x_{n}$ lead (precede) or follow (lag) $x_{m}$ throughout time ? To illustrate the meaning of this question, we consider the signals $f,g:\mathbb{R}\rightarrow\mathbb{R}$ defined by $f(t)=\exp(-\pi t^{2})$ and $g(t)=f(t-1).$ We plot these signals in Figure 1.1. Note that $g$ is a horizontal translation of $f$ to the right by $1$ unit of time, and the number $1$ is the value of the time shift. As a result, throughout time, we can describe the evolution of $f$ and $g$ in relation to one another: when $f$ increases, $g$ later increases, and when $f$ decreases, $g$ later decreases. Thus, we say that $f$ leads $g,$ or equivalently, $g$ follows $f.$

Refer to caption — Figure 1.1: An illustration of two signals exhibiting a leader follower relationship. Here, we plot $f(t)=e^{-\pi t^{2}}$ and $g(t)=f(t-1).$ We interpret $f$ to lead $g.$

Next, assume all of the $N$ component signals of $\mathbf{x}$ trace some underlying real-valued signal up to scaling constants and time shifts. What is the order in which the component signals evolve throughout time ? More specifically, what is the order of the time shifts ? To illustrate the meaning of this question, let $T$ be the time interval $[0,1],$ and consider the sinusoidal signal $\mathbf{x}:[0,1]\rightarrow\mathbb{R}^{N},$ whose $n$ -th component signal is defined by

\displaystyle x_{n}(t)

\displaystyle=\sin\left(2\pi\left(t-\frac{n-1}{N}\right)\right)

(1.1)

for each $1\leq n\leq N.$ In Figure 1.2, we fix $N=20$ and plot some of these component sinusoidal signals. Each $x_{n}$ is a time-shifted copy of the signal $\phi(t)=\sin(2\pi t)$ with $\frac{n-1}{N}$ being the value of the time shift. Since $\frac{n-1}{N}<\frac{n}{N}$ for each $1\leq n\leq N,$ we interpret the $n$ -th component signal as leading the $(n+1)$ -th component signal. Thus, the order in which the $N$ component signals evolve throughout time is $x_{1}\ ,\ \dots\ ,\ x_{N}.$

1.2 Analysis of Periodic Signals

In this section, we focus on analyzing periodic signals. Recall that a signal $\mathbf{x}:\mathbb{R}\rightarrow\mathbb{R}^{N}$ is periodic if there exists a minimal constant $P>0$ such that $\mathbf{x}(t+P)=\mathbf{x}(t)$ for all $t\in\mathbb{R}.$ The constant $P$ is known as the period of $\mathbf{x}.$ We specifically say $\mathbf{x}$ is $P$ -periodic in order to emphasize $\mathbf{x}$ has period $P$ . Equivalently, we can view a $P$ -periodic signal as a well-defined map between $\mathbb{R}/(P\mathbb{Z})$ and $\mathbb{R}^{N},$ where $\mathbb{R}/(P\mathbb{Z})$ is the additive group $[0,P)$ equipped with addition modulo $P.$

To analyze periodic signals in general, one uses tools from harmonic analysis, which include Fourier series [1, Chapter 2] and the Fourier transform [1, Chapter 5]. Another common tool that is suitable for analyzing periodic signals is the cross-correlation function [2, Chapter 12]. Recall if $f,g:\mathbb{R}\rightarrow\mathbb{R}$ are two $P$ -periodic signals, then the (periodic) cross-correlation of $f$ and $g$ is the $P$ -periodic signal $f\star g:\mathbb{R}\rightarrow\mathbb{R}$ defined by

\displaystyle(f\star g)(\tau)

\displaystyle=\int_{0}^{P}f(t-\tau)\ g(t)\ dt.

(1.2)

The cross-correlation function measures the similarity between the input signals $f$ and $g$ as a function of the displacement variable $\tau.$ The displacement at which the cross-correlation attains its absolute maximum reveals information about $f$ and $g$ . Up to an integer multiple of the period, the maximum displacement coincides with the time shift value such that shifting $f$ by this value would result in the signals $f$ and $g$ being most aligned.

For example, let $f,g:\mathbb{R}\rightarrow\mathbb{R}$ be the periodic signals $f(t)=\sin(2\pi t)$ and $g(t)=\cos(2\pi t).$ Explicitly computing the cross-correlation of $f$ and $g$ using (1.2), we obtain

	$\displaystyle(f\star g)(\tau)$	$\displaystyle=\int_{0}^{1}\sin(2\pi(t-\tau))\cos(2\pi t)\ dt$
		$\displaystyle=\int_{0}^{1}\left(\sin(2\pi t)\cos(2\pi\tau)-\sin(2\pi\tau)\cos(% 2\pi t)\right)\cos(2\pi t)\ dt$
		$\displaystyle=\int_{0}^{1}\sin(2\pi t)\cos(2\pi t)\cos(2\pi\tau)-\sin(2\pi\tau% )\cos^{2}(2\pi t)\ dt$
		$\displaystyle=\int_{0}^{1}\frac{\sin(4\pi t)}{2}\cos(2\pi\tau)-\sin(2\pi\tau)% \left(\frac{1-\cos(4\pi t)}{2}\right)\ dt$
		$\displaystyle=-\frac{\sin(2\pi\tau)}{2}.$

We plot this cross-correlation function in Figure 1.3. On one hand, note that $g$ is a time-shifted copy of $f$ to the right by $\frac{3}{4}$ units of time. On the other hand, we observe the cross-correlation function $f\star g$ attains its maximum at $\tau=\frac{3}{4}+k$ for any $k\in\mathbb{Z}.$ This means the displacement at which the cross-correlation function attains its maximum coincides with the time shift between $f$ and $g$ up to an integer.

1.2.1 The Chain of Offsets Model

We answer the leader follower dynamics questions that we posed in Section 1.1 for periodic signals. We assume that $\mathbf{x}:\mathbb{R}\rightarrow\mathbb{R}^{N}$ obeys a Chain of Offsets model (COOM), which was first introduced in [3]. This model stipulates the existence of a $P$ -periodic signal $\phi:\mathbb{R}\rightarrow\mathbb{R},$ scaling constants $c_{1}\ ,\ \dots\ ,\ c_{N}>0,$ and offsets $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{N}\in\mathbb{R}/(P\mathbb{Z})$ such that

\displaystyle x_{n}(t)

\displaystyle=c_{n}\phi(t-\alpha_{n})

(1.3)

for each $1\leq n\leq N.$ In other words, under COOM, all component signals of $\mathbf{x}$ trace an underlying real-valued periodic signal up to scaling constants and offsets.

Answering our posed leader follower dynamics questions amounts to determining the order of the offsets $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{N}$ under COOM. However, we emphasize this order of offsets is a cyclic order; it corresponds to the order in which the points $e^{2\pi i\alpha_{1}}\ ,\ \dots\ ,\ e^{2\pi i\alpha_{N}}$ are traversed on the unit circle in the counterclockwise direction. We represent the cyclic order as a permutation $\sigma\in S_{N}$ of the indices $\left\{1\ ,\ \dots\ ,\ N\right\}$ satisfying the condition that $e^{2\pi i\alpha_{\sigma(n)}}$ comes before $e^{2\pi i\alpha_{\sigma(n+1)}}$ in the counterclockwise direction. In terms of the signals themselves, we interpret $x_{\sigma(n)}$ as (cyclically) preceding $x_{\sigma(n+1)}$ for each $1\leq n<N.$ Furthermore, $\sigma$ is unique up to cyclic shifts of indices. This means if $\widetilde{\sigma}\in S_{N}$ is another permutation representing the cyclic order of the offsets, then there exists an integer $d$ such that

\displaystyle\widetilde{\sigma}(n)=\sigma(d+n-1)

for all $1\leq n\leq N,$ in which $d+n-1$ is indexed mod $N.$

We give an example of how a permutation representing the cyclic order is unique up to cyclic shifts. Reconsider the sinsoidal signal $\mathbf{x}$ with $n$ -th component signal defined in (1.1) for each $1\leq n\leq N.$ Note that $\mathbf{x}$ satisfies COOM in (1.3) if we set $\phi(t)=\sin(2\pi t)$ and $c_{n}=1$ and $\alpha_{n}=\frac{n-1}{N}$ for each $1\leq n\leq N.$ On one hand, note that $0=\alpha_{1}<\ \dots\ <\alpha_{N}<1,$ which means the points $e^{2\pi i\alpha_{1}}\ ,\ \dots\ ,\ e^{2\pi i\alpha_{N}}$ are traversed in that order on the unit circle in the counterclockwise direction beginning from the positive real axis. Hence, the identity permutation $\text{Id}\in S_{N}$ is one permutation that represents the cyclic order of the offsets. In fact, any permutation $\widetilde{\sigma}\in S_{N}$ defined by $\widetilde{\sigma}(n)=n+d,$ where $d\in\mathbb{N}$ is a fixed integer and $n+d$ is indexed mod $N,$ is a valid permutation representing the cyclic order of the offsets. Nevertheless, observe that

\displaystyle\text{Id}(n)

\displaystyle=\widetilde{\sigma}(n+N-d),

(1.4)

which means the identity permutation is equal to $\widetilde{\sigma}$ up to $N-d$ cyclic shifts.

1.2.2 Recovery of the Cyclic Order under COOM

Given the $P$ -periodic signal $\mathbf{x},$ we outline a procedure to recover the cyclic order of the offsets $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{N}$ under COOM. This procedure involves the utilization of the cross-correlation function, which we defined previously in (1.7).

First, we write the underlying $P$ -periodic signal $\phi$ in COOM as a Fourier series:

\displaystyle\phi(t)

\displaystyle=\sum_{k\in\mathbb{Z}}\hat{\phi}_{k}\ e^{\frac{2\pi ikt}{P}},

(1.5)

where $\hat{\phi}_{k}$ is the $k$ -th Fourier coefficient of $\phi$ for each $k\in\mathbb{Z}.$ Then, for each $1\leq m,n\leq N,$ the cross-correlation of the component signals $x_{m}$ and $x_{n}$ is explicitly

	$\displaystyle(x_{m}\star x_{n})(\tau)$	$\displaystyle=\int_{0}^{P}x_{m}(t-\tau)\ x_{n}(t)\ dt$
		$\displaystyle=\int_{0}^{P}\sum_{k,\ell\in\mathbb{Z}}c_{m}c_{n}\ \hat{\phi}_{k}% \ \hat{\phi_{\ell}}\ e^{\frac{2\pi ik(t-\tau-\alpha_{m})+2\pi i\ell(t-\alpha_{% n})}{P}}\ dt$
		$\displaystyle=c_{m}c_{n}\sum_{k,\ell\in\mathbb{Z}}\hat{\phi}_{k}\ \hat{\phi_{% \ell}}\ e^{\frac{-2\pi ik(\tau+\alpha_{m})-2\pi i\ell\alpha_{n}}{P}}\int_{0}^{% P}e^{\frac{2\pi i(k+\ell)t}{P}}dt$
		$\displaystyle=c_{m}c_{n}\sum_{k,\ell\in\mathbb{Z}}\hat{\phi}_{k}\ \hat{\phi_{% \ell}}\ e^{\frac{-2\pi ik(\tau+\alpha_{m})-2\pi i\ell\alpha_{n}}{P}}\ P\ % \delta_{k,-\ell}$
		$\displaystyle=Pc_{m}c_{n}\sum_{k\in\mathbb{Z}}\left\|\hat{\phi}_{k}\right\|^{2}% \ e^{\frac{2\pi ik(\alpha_{n}-\alpha_{m}-\tau)}{P}}$
		$\displaystyle=Pc_{m}c_{n}\left\|\hat{\phi}_{0}\right\|^{2}+2Pc_{m}c_{n}\sum_{k% \in\mathbb{N}}\left\|\hat{\phi}_{k}\right\|^{2}\cos\left(\frac{2\pi k(\alpha_{n}% -\alpha_{m}-\tau)}{P}\right),$

where $\delta_{\cdot,\cdot}$ is the Kronecker delta function. Differentiating the cross-correlation function with respect to $\tau$ , we obtain

\displaystyle(x_{m}\star x_{n})^{\prime}(\tau)

\displaystyle=-4\pi c_{m}c_{n}\sum_{k\in\mathbb{N}}k\left|\hat{\phi}_{k}\right% |^{2}\sin\left(\frac{2\pi k(\alpha_{n}-\alpha_{m}-\tau)}{P}\right).

Equating the derivative to $0,$ we see that any critical point of $x_{m}\star x_{n}$ is of the form

\displaystyle\tau

\displaystyle=\alpha_{n}-\alpha_{m}+kP,

where $k\in\mathbb{Z}.$ This means that the critical points of the $x_{m}\star x_{n}$ are equal to the difference between the offsets $\alpha_{n}$ and $\alpha_{m}$ up to integer multiples of the period $P.$

Now, we are ready to construct a permutation of the indices $\left\{1\ ,\ \dots\ ,\ N\right\}$ representing the cyclic order of the offsets under COOM. For each $1\leq m,n\leq N,$ let

\displaystyle\tau_{m,n}

\displaystyle=\alpha_{n}-\alpha_{m},

in which $\alpha_{n}-\alpha_{m}$ is specifically reduced modulo $P.$ More explicitly, we choose an integer $k\in\mathbb{Z}$ such that $\alpha_{n}-\alpha_{m}+kP\in[0,P)$ and set $\tau_{m,n}$ equal to $\alpha_{n}-\alpha_{m}+kP$ for this choice of $k.$ Define the permutation $\sigma\in S_{N}$ be a permutation as follows: let $\sigma(1)$ be arbitrary, and for each $1<n\leq N,$ let $\sigma(n)$ be the smallest index not in $\left\{\sigma(1)\ ,\ \dots\ ,\ \sigma(n-1)\right\}$ such that $\tau_{\sigma(n-1)\ ,\ \sigma(n)}$ is minimal i.e. we have $\tau_{\sigma(n-1)\ ,\ \sigma(n)}\leq\tau_{\sigma(n-1)\ ,\ m}$ for all $m\notin\left\{\sigma(1)\ ,\ \dots\ ,\ \sigma(n-1)\right\}.$ Observe $\sigma$ ensures that $e^{2\pi i\alpha_{\sigma(n)}}$ is the closest point on the unit circle to the point $e^{2\pi i\alpha_{\sigma(n-1)}}$ in the counterclockwise direction. Therefore, $\sigma$ is a candidate permutation representing the cyclic order.

It remains to show that our constructed permutation $\sigma$ from the previous paragraph is unique up to cyclic shifts. To this end, let $\widetilde{\sigma}\in S_{N}$ be another permutation constructed according to the procedure in the previous paragraph. Note that there exists a unique index $1\leq d\leq N$ such that $\widetilde{\sigma}(1)=\sigma(d),$ which follows from the fact that the composition $\widetilde{\sigma}\circ\sigma^{-1}$ is itself a permutation in $S_{N}.$ We claim that for this chosen index $d,$ we have

\displaystyle\widetilde{\sigma}(n)

\displaystyle=\sigma(d+n-1)

(1.6)

for all $1\leq n\leq N,$ in which $d+n-1$ is indexed mod $N.$ We prove this by induction. The statement (1.6) already holds for $n=1$ by the property of $d.$ Now, suppose the statement (1.6) holds for $n=k$ for some $1<k<N.$ We show it must also hold for $n=k+1.$ By the construction of $\widetilde{\sigma},$ we have that $\widetilde{\sigma}(k+1)$ is the smallest index not in $\left\{\widetilde{\sigma}(1)\ ,\ \dots\ ,\ \widetilde{\sigma}(k)\right\}$ such that $\tau_{\widetilde{\sigma}(k)\ ,\ \widetilde{\sigma}(k+1)}$ is minimal. On the other hand, by the induction hypothesis, $\widetilde{\sigma}(k+1)$ is the smallest index not in $\left\{\sigma(1)\ ,\ \dots\ ,\ \sigma(d+k-1)\right\}$ such that $\tau_{\sigma(d+k-1)\ ,\ \widetilde{\sigma}(k+1)}$ is minimal. By the construction of $\sigma,$ we conclude $\widetilde{\sigma}(k+1)=\sigma(d+k).$ Therefore, the statement (1.6) holds for $n=k+1.$ Thus, we have shown our candidate permutation $\sigma$ is unique up to cyclic shifts.

1.3 Analysis of Non-Periodic Signals

In the previous section, we utilized the cross-correlation function to answer our posed leader follower dynamics questions in Section 1.1 for periodic signals. In this section, we investigate whether cross-correlation can be used to answer our posed leader follower dynamics questions for signals that are not necessarily periodic. The cross-correlation of two generic signals $f,g:\mathbb{R}\rightarrow\mathbb{R}$ is defined as the map $f\star g:\mathbb{R}\rightarrow\mathbb{R},$ where

\displaystyle(f\star g)(\tau)

\displaystyle=\int_{\mathbb{R}}f(t-\tau)g(t)\ dt,

(1.7)

provided the integral on the right hand side exists.

1.3.1 Time Shifted Model

Suppose $f:\mathbb{R}\rightarrow\mathbb{R}$ is a signal that is compactly supported in a closed interval i.e. the closure of the set of all points $t\in\mathbb{R}$ with $f(t)\neq 0$ is a closed interval. Let $g:\mathbb{R}\rightarrow\mathbb{R}$ be the map defined by

\displaystyle g(t)

\displaystyle=f(t-\tau_{0}).

(1.8)

for some time shift $\tau_{0}\in\mathbb{R}.$ We refer to the model (1.8) as a time shifted model. Determining the leader follower relationship between $f$ and $g$ under the time shifted model amounts to determining the sign of the time shift $\tau_{0}.$ If $\tau_{0}>0,$ then we interpret $f$ to lead $g.$ If $\tau_{0}<0,$ then we interpret $f$ to follow $g.$

Given the signals $f$ and $g,$ we determine the value of the time shift $\tau_{0}$ under the time shifted model (1.8). In order to do this, we utilize the cross-correlation of $f$ and $g$ defined in (1.7). For each $\tau\in\mathbb{R},$ we have

$\displaystyle(f\star g)(\tau)$	$\displaystyle=\int_{\mathbb{R}}f(t-\tau)\ g(t)\ dt$
	$\displaystyle=\int_{\mathbb{R}}\frac{(f(t-\tau))^{2}+(g(t))^{2}-(f(t-\tau)-g(t% ))^{2}}{2}\ dt$
	$\displaystyle=\int_{\mathbb{R}}\frac{(f(t-\tau))^{2}+(f(t-\tau_{0}))^{2}}{2}\ % dt-\int_{\mathbb{R}}\frac{\left(f(t-\tau)-f(t-\tau_{0})\right)^{2}}{2}\ dt$
	$\displaystyle=\int_{\mathbb{R}}\frac{(f(t-\tau))^{2}}{2}\ dt+\int_{\mathbb{R}}% \frac{(f(t-\tau_{0}))^{2}}{2}\ dt-\int_{\mathbb{R}}\frac{\left(f(t-\tau)-f(t-% \tau_{0})\right)^{2}}{2}\ dt$	(1.9)
	$\displaystyle=(f\star f)(0)-\int_{\mathbb{R}}\frac{\left(f(t-\tau)-f(t-\tau_{0% })\right)^{2}}{2}\ dt,$	(1.10)

in which the first term of (1.10) follows from observing the first two integral terms in (1.9) are both equal to $\frac{(f\star f)(0)}{2}$ via the respective changes of variables $t\mapsto t-\tau$ and $t\mapsto t-\tau_{0}.$

Note that the integral after the minus sign in (1.9) is a non-negative quantity because it is the integral of a non-negative function. This implies the cross-correlation function has an upper bound

\displaystyle(f\star g)(\tau)\leq(f\star f)(0).

(1.11)

for all $\tau\in\mathbb{R}.$ Equality in (1.11) holds if and only if the integral term in (1.10) is $0$ for some $\tau\in\mathbb{R}.$ Because $f$ is continuous and compactly supported in a closed interval, equality in (1.11) holds if and only if the integrand of the integral term in (1.10) is $0$ i.e.

\displaystyle f(t-\tau)

\displaystyle=f(t-\tau_{0}).

for all $t\in\mathbb{R}.$ But since $f$ is not periodic, we have $(f\star g)(\tau)=(f\star f)(0)$ if and only if $\tau=\tau_{0}.$ Therefore, the unique displacement at which the cross-correlation $f\star g$ attains its maximum is precisely the time shift $\tau_{0}$ under the model (1.8). Having recovered the time shift $\tau_{0},$ we can determine the leader follower relationship between the signals $f$ and $g$ upon inspecting the sign of $\tau_{0}.$

1.3.2 Time Reparameterization Invariance

Based on our analysis so far, it seems that the cross-correlation function is a useful signal processing tool that enables us to determine the leader follower relationship between any two input signals, provided these signals have nice enough properties. We wonder whether there are drawbacks to using the cross-correlation.

We focus on one drawback central to the thesis. The cross-correlation function suffers under time reparameterization. A time reparameterization of the time interval $T$ is a strictly increasing homeomorphism between $T$ and itself. This means that if we use the cross-correlation as a tool to analyze signals, then the results given by cross-correlation may not be the same across all time reparameterizations. In other words, the results we obtain via cross-correlation may vary according to how we observe the signals throughout time.

To illustrate with an example, let $T=[-3,2],$ and consider the signals $f:[-3,2]\rightarrow\mathbb{R},$ where

\displaystyle f(t)

\displaystyle=\begin{cases}\sin^{2}(2\pi t)&t\in[-1,1]\\ 0&\text{else}\end{cases}

(1.12)

and $g(t)=f\left(t-\frac{1}{4}\right).$ Consider the map $r:[-3,2]\rightarrow[-3,2]$ defined by

\displaystyle r(t)=\frac{t^{3}+6}{7}.

(1.13)

Note that $r$ is a valid time reparameterization of $[-3,2].$ Indeed, we have $r(-3)=-3$ and $r(2)=2,$ and since $t\mapsto t^{3}$ is strictly increasing on $\mathbb{R}$ , the map $t\mapsto r(t)$ is also strictly increasing on $\mathbb{R}$ and thus strictly increasing in $[-3,2].$ Therefore, $r([-3,2])=[-3,2].$ Finally, the inverse map $r^{-1}(y)=(7y-6)^{\frac{1}{3}}$ is continuous on all of $\mathbb{R}$ and is thus continuous in $[-3,2].$ Therefore, $r$ is a strictly increasing homeomorphism, which means it is a valid time reparameterization. However, in Figure 1.4, we plot the correlation functions $f\star g$ and $(f\circ r)\star(g\circ r).$ We see that the displacement maximum for the original cross-correlation $f\star g$ and the displacement maximum for the reparameterized cross-correlation $(f\circ r)\star(g\circ r)$ differ in sign. Hence, the cross-correlation function does not allow us to conclude whether $f$ leads $g$ or vice versa.

Thus, in order to perform any meaningful analysis of signals, we need to consider using tools that are time reparameterization invariant.

1.4 Cyclicity Analysis

In this section, we consider cyclic signals. A signal $\mathbf{x}:\mathbb{R}\rightarrow\mathbb{R}^{N}$ is cyclic if $\mathbf{x}\circ r$ is periodic for some time reparameterization $r:\mathbb{R}\rightarrow\mathbb{R}.$ Cyclic signals generalize periodic signals. Roughly speaking, cyclic signals have repeating temporal patterns, but such patterns do not repeat in a predictable fashion, which makes them different from periodic signals. Many real world phenomena are cyclic in nature. For example, the economy undergoes expansion and recession, but not predictably due to various external factors that control the economy. Other examples include cardiac cycles and fMRI signals [4].

We illustrate the difference between a periodic and a cyclic signal with a theoretical example. Consider the signals $f(t)=\sin(2\pi t)$ and $g(t)=f(t^{3})$ on $\mathbb{R}.$ In Figure 1.5, we plot $f$ and $g$ and the resulting signals we observe after shifting the time axis by $1$ unit. Since $f$ is $1$ -periodic, shifting the time axis by $1$ unit results in us observing the exact same signal that we started with. On the other hand, $g$ is not $1$ -periodic; shifting the time axis by $1$ unit results in us observing a completely different signal than the one we started with. In fact, $g$ is not $P$ -periodic for any constant $P>0.$ However, $g$ is cyclic: if $r:\mathbb{R}\rightarrow\mathbb{R}$ is defined by $r(t)=t^{\frac{1}{3}},$ then note that $r$ is a valid time reparameterization satisfying $f=g\circ r.$

In this section, we describe the procedure of Cyclicity Analysis, first introduced in [3], which is the study of the leader-follower dynamics pertaining to cyclic signals. Cyclicity Analysis directly answers both of our posed leader follower questions in Section 1.1 in a time reparameterization invariant fashion.

1.4.1 Iterated Path Integrals

First, we assume $\mathbf{x}:T\rightarrow\mathbb{R}^{N}$ is a generic signal in which $T=[a,b]$ is a closed interval. We recall the first question we posed in Section 1.1, which is to determine the pairwise leader follower relationship between the $m$ -th component and $n$ -th component signals of $\mathbf{x}$ for fixed indices $1\leq m,n\leq N,$ assuming such component signals evolve similarly throughout time.

In order to answer this question, Cyclicity Analysis invokes the theory of iterated path integrals [5, 6, 7, 8]. For each $K\in\mathbb{N}$ , a $K$ -th order iterated path integral of $\mathbf{x}$ is a functional $\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}$ of the form

\displaystyle\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{x})

\displaystyle=\frac{1}{K!}\int_{a\leq t_{1}\leq\ \dots\ \leq t_{K}\leq b}x^{% \prime}_{i_{1}}(t_{1})\ \dots\ x^{\prime}_{i_{K}}(t_{K})\ dt_{K}\ \dots\ dt_{1},

(1.14)

where $1\leq i_{1}\leq\ \dots\ \leq i_{K}\leq N.$ For example, a first order iterated path integral of $\mathbf{x}$ is explicitly of the form

	$\displaystyle\text{IP}_{n}(\mathbf{x})$	$\displaystyle=\int_{a}^{b}x^{\prime}_{n}(t)\ dt$
		$\displaystyle=x_{n}(b)-x_{n}(a),$

which is the increment of the trajectory parameterized by $x_{n}:T\rightarrow\mathbb{R}.$

Iterated path integrals date back to the times of Riemann and Picard. They were initially used as tools to construct and analyze solutions to ordinary differential equations. Later on, K.T. Chen [5, 6, 7, 8] discovered the relevance of iterated path integrals to algebraic topology, while Terry Lyons [9, 10] discovered the relevance to the analysis of rough paths, which are Holder trajectories.

Iterated path integrals of the signal $\mathbf{x}$ obey several invariant properties [3], which make them a suitable tool for data analysis. Firstly, they are indeed time reparameterization invariant: for any time reparameterization $r:T\rightarrow T,$ we have

\displaystyle\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{x})

\displaystyle=\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{x}\circ r),

which follows directly from performing a change of variables $t\mapsto r(t)$ in (1.14). Next, the iterated path integrals of $\mathbf{x}$ are invariant under parallel translation: if $\mathbf{a}\in\mathbb{R}^{N}$ is fixed and if $\mathbf{y}:T\rightarrow\mathbb{R}^{N}$ is the signal defined by $\mathbf{y}(t)=\mathbf{x}(t)+\mathbf{a},$ then

\displaystyle\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{x})

\displaystyle=\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{y}).

Moreover, the iterated path integrals are invariant under detours. If $C_{1}\ ,\ \dots\ ,\ C_{L}$ are oriented curves in $\mathbb{R}^{N}$ such that the ending point of $C_{\ell}$ is the starting point of $C_{\ell+1},$ their concatenation possesses a detour if for some $1\leq\ell<L,$ the curves $C_{\ell+1}$ and $C_{\ell}$ are the same but with opposite orientations. In other words, if we traverse $C_{1}\ ,\ \dots\ ,\ C_{L}$ in that order, then our traversal will involve backtracking. If $\mathbf{x}$ and $\mathbf{y}$ parameterize the same curve, in which $\mathbf{x}$ traverses it with detours but $\mathbf{y}$ traverses it without detours, then

\displaystyle\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{x})

\displaystyle=\text{IP}_{i_{1}\ ,\ \dots\ ,\ i_{K}}(\mathbf{y}).

Finally, K.T. Chen [5, 6, 7, 8] proved a stronger statement: two signals $\mathbf{x},\mathbf{y}:T\rightarrow\mathbb{R}^{N}$ that parameterize their respective curves without detours are the same up to time reparameterization and parallel translation if and only if all of their corresponding iterated path integrals agree.

1.4.2 Oriented Areas

Now, assume $\mathbf{x}:T\rightarrow\mathbb{R}^{N}$ parameterizes a closed curve i.e. $\mathbf{x}(a)=\mathbf{x}(b).$ We consider the second order iterated path integrals of $\mathbf{x}.$ For each $1\leq m,n\leq N,$ we substitute $K=2$ and $i_{1}=m$ and $i_{2}=n$ into (1.14) to obtain the explicit form of the second order iterated path integral:

$\displaystyle\text{IP}_{m,n}(\mathbf{x})$	$\displaystyle=\frac{1}{2}\int_{a\leq s\leq t\leq b}x_{m}^{\prime}(s)\ x_{n}^{% \prime}(t)\ ds\ dt$
	$\displaystyle=\frac{1}{2}\int_{a}^{b}\int_{a}^{t}x_{m}^{\prime}(s)\ x_{n}^{% \prime}(t)\ ds\ dt$
	$\displaystyle=\frac{1}{2}\int_{a}^{b}(x_{m}(t)-x_{m}(a))\ x_{n}^{\prime}(t)\ dt$
	$\displaystyle=\frac{1}{2}\int_{a}^{b}x_{m}(t)\ x_{n}^{\prime}(t)\ dt-\frac{% \left(x_{n}(b)-x_{n}(a)\right)\ x_{m}(a)}{2}$	(1.15)
	$\displaystyle=\frac{1}{2}\int_{a}^{b}x_{m}(t)\ x_{n}^{\prime}(t)\ dt,$

in which the second term of (1.15) vanishes due to the assumption $\mathbf{x}$ parameterizes a closed curve.

For all $1\leq m,n\leq N,$ we define

	$\displaystyle A_{m,n}(\mathbf{x})$	$\displaystyle=\text{IP}_{m,n}(\mathbf{x})-\text{IP}_{n,m}(\mathbf{x})$
		$\displaystyle=\frac{1}{2}\int_{a}^{b}x_{m}(t)\ x_{n}^{\prime}(t)-x_{n}(t)\ x_{% m}^{\prime}(t)\ dt.$		(1.16)

By Green’s Theorem, the integral (1.16) is a line integral that is equal to the area of the region enclosed by the curve parameterized by the map $(x_{m},x_{n}):T\rightarrow\mathbb{R}^{2}.$ We note, however, this area is actually an oriented (signed) area. In particular, $A_{m,n}(\mathbf{x})$ is positive if $(x_{m},x_{n}):T\rightarrow\mathbb{R}^{2}$ traverses the aforementioned curve in the counterclockwise direction and is negative otherwise. We heuristically interpret the oriented area as an indicator of the lead-lag relationship between the signals $x_{m}$ and $x_{n}$ . We say that $x_{m}$ leads (precedes) $x_{n}$ if $A_{m,n}(\mathbf{x})>0.$ We say that $x_{m}$ follows (lags) $x_{n}$ if $A_{m,n}(\mathbf{x})<0.$ Otherwise, we say $x_{m}$ and $x_{n}$ are in sync.

We illustrate the oriented area with an example. Reconsider the signals $f,g:[-\infty,\infty]\rightarrow\mathbb{R}$ with $f(t)=e^{-\pi t^{2}}$ and $g(t)=f(t-1)$ with $f(\pm\infty)=g(\pm\infty)=0.$ The curve $(f,g):[-\infty,\infty]\rightarrow\mathbb{R}^{2}$ is a closed curve oriented in the counterclockwise direction, which we plot in Figure 1.6. The oriented area of the region enclosed by this curve is explicitly

	$\displaystyle\frac{1}{2}\int_{\mathbb{R}}f(t)\ g^{\prime}(t)-f^{\prime}(t)\ g(% t)\ dt$	$\displaystyle=\frac{1}{2}\int_{\mathbb{R}}-2\pi(t-1)e^{-\pi t^{2}}e^{-\pi(t-1)% ^{2}}+2\pi te^{-\pi t^{2}}e^{-\pi(t-1)^{2}}\ dt$
		$\displaystyle=\int_{\mathbb{R}}\pi e^{-\pi t^{2}-\pi(t-1)^{2}}\ dt=\frac{\pi e% ^{-\frac{\pi}{2}}}{\sqrt{2}},$

which is positive. Therefore, we interpret $f$ to lead $g$ based on the oriented area.

Recall our assumption $\mathbf{x}$ parameterizes a closed curve. We can modify the definition of the oriented area in the case $\mathbf{x}$ does not parameterize a closed curve. In this situation, we may consider any map $\widetilde{\mathbf{x}}:T\rightarrow\mathbb{R}^{N}$ that parameterizes a closed curve formed by concatenating the curve parameterized by the original signal $\mathbf{x}$ and any oriented curve in $\mathbb{R}^{N}$ with starting point $\mathbf{x}(b)$ and ending point $\mathbf{x}(a).$ Different choices of $\widetilde{\mathbf{x}}$ yield different values of $A_{m,n}(\widetilde{\mathbf{x}}).$ For our purposes, however, we will choose a signal $\widetilde{\mathbf{x}}:T\rightarrow\mathbb{R}^{N}$ such that

\displaystyle A_{m,n}(\widetilde{\mathbf{x}})

\displaystyle=\frac{1}{2}\int_{a}^{b}x_{m}(t)\ x_{n}^{\prime}(t)-x_{n}(t)\ x_{% m}^{\prime}(t)\ dt.

(1.17)

i.e. the $(m,n)$ -th oriented area of $\widetilde{\mathbf{x}}$ is the same as the $(m,n)$ -th line integral of $\frac{x_{m}x_{n}^{\prime}-x_{n}x_{m}^{\prime}}{2}$ over the curve parameterized by $(x_{m},x_{n}):T\rightarrow\mathbb{R}^{2}.$

Consider the signal $\widetilde{\mathbf{x}}:T\rightarrow\mathbb{R}^{N}$ defined by

\displaystyle\widetilde{\mathbf{x}}(t)

\displaystyle=\begin{cases}\mathbf{x}(2t-a)&t\in\left[a,\frac{a+b}{2}\right]\\ \mathbf{x}(b)+\left(\mathbf{x}(a)-\mathbf{x}(b)\right)\left(\frac{a+b-2t}{a-b}% \right)&t\in\left[\frac{a+b}{2},b\right].\end{cases}

(1.18)

This signal parameterizes the curve formed by concatenating the curve parameterized by the original signal $\mathbf{x}$ and the oriented line segment with starting point $\mathbf{x}(b)$ and ending point $\mathbf{x}(a).$ We show this signal satisfies our requirement (1.17). Explicitly computing $A_{m,n}(\widetilde{\mathbf{x}}),$ we have

$\displaystyle A_{m,n}(\widetilde{\mathbf{x}})$	$\displaystyle=\frac{1}{2}\int_{a}^{b}\widetilde{x}_{m}(t)\ \widetilde{x}_{n}^{% \prime}(t)-\widetilde{x}_{n}(t)\ \widetilde{x}_{m}^{\prime}(t)\ dt$
	$\displaystyle=\frac{1}{2}\int_{a}^{\frac{a+b}{2}}\widetilde{x}_{m}(t)\ % \widetilde{x}_{n}^{\prime}(t)-\widetilde{x}_{n}(t)\ \widetilde{x}_{m}^{\prime}% (t)\ dt$
	$\displaystyle\qquad\qquad+\frac{1}{2}\int_{\frac{a+b}{2}}^{b}\widetilde{x}_{m}% (t)\ \widetilde{x}_{n}^{\prime}(t)-\widetilde{x}_{n}(t)\ \widetilde{x}_{m}^{% \prime}(t)\ dt$
	$\displaystyle=\int_{a}^{\frac{a+b}{2}}x_{m}(2t-a)\ x_{n}^{\prime}(2t-a)-x_{n}(% 2t-a)\ x_{m}^{\prime}(2t-a)\ dt$
	$\displaystyle\qquad-\frac{x_{n}(a)-x_{n}(b)}{b-a}\int_{\frac{a+b}{2}}^{b}x_{m}% (b)+\left(x_{m}(a)-x_{m}(b)\right)\left(\frac{a+b-2t}{a-b}\right)\ dt$
	$\displaystyle\qquad\qquad+\frac{x_{m}(a)-x_{m}(b)}{b-a}\int_{\frac{a+b}{2}}^{b% }x_{n}(b)+\left(x_{n}(a)-x_{n}(b)\right)\left(\frac{a+b-2t}{a-b}\right)\ dt$
	$\displaystyle=\int_{a}^{\frac{a+b}{2}}x_{m}(2t-a)\ x_{n}^{\prime}(2t-a)-x_{n}(% 2t-a)\ x_{m}^{\prime}(2t-a)\ dt$
	$\displaystyle\qquad-\frac{x_{m}(b)(x_{n}(a)-x_{n}(b))}{2}-\frac{(x_{m}(a)-x_{m% }(b))(x_{n}(a)-x_{n}(b))}{4}$
	$\displaystyle\qquad\quad+\frac{x_{n}(b)(x_{m}(a)-x_{m}(b))}{2}+\frac{(x_{n}(a)% -x_{n}(b))(x_{m}(a)-x_{m}(b))}{4}$
	$\displaystyle=\int_{a}^{\frac{a+b}{2}}x_{m}(2t-a)\ x_{n}^{\prime}(2t-a)-x_{n}(% 2t-a)\ x_{m}^{\prime}(2t-a)\ dt$	(1.19)
	$\displaystyle=\frac{1}{2}\int_{a}^{b}x_{m}(t)\ x_{n}^{\prime}(t)-x_{n}(t)\ x_{% m}^{\prime}(t)\ dt,$	(1.20)

in which (1.20) follows from performing the change of variables $t\mapsto 2t-a$ in (1.19). Thus, our choice of $\widetilde{\mathbf{x}}$ satisfies our specific requirement in (1.17).

So from henceforth, if the signal $\mathbf{x}$ does not parameterize a closed curve, then we define $A_{m,n}(\mathbf{x})=A_{m,n}(\widetilde{\mathbf{x}}),$ where $\widetilde{\mathbf{x}}$ is the signal defined in (1.18).

1.4.3 The Lead Matrix

If $\mathbf{x}:T\rightarrow\mathbb{R}^{N}$ is a signal with $T=[a,b],$ then we will write $A_{m,n}$ as the oriented area in place of the original notation $A_{m,n}(\mathbf{x}).$ In isolation, the oriented area is simply a number that heuristically indicates the pairwise leader follower relationship between any two specific components of $\mathbf{x}.$ For future purposes, however, we would like to consider all of these $N^{2}$ total oriented areas together. So we will arrange all the oriented areas into an $N\times N$ matrix. Let $\mathbf{A}$ be the $N\times N$ matrix whose $(m,n)$ -th entry is the oriented area $A_{m,n}$ for each $1\leq m,n\leq N.$ We refer to $\mathbf{A}$ as the lead(-lag) matrix corresponding to $\mathbf{x}$ . More succinctly, we can rewrite $\mathbf{A}$ as an integral of a matrix-valued function:

\displaystyle\mathbf{A}

\displaystyle=\frac{1}{2}\int_{a}^{b}\mathbf{x}(t)\ \dot{\mathbf{x}}^{\text{T}% }(t)-\dot{\mathbf{x}}(t)\ \mathbf{x}^{\text{T}}(t)\ dt,

(1.21)

where $\dot{\mathbf{x}}(t)=\left(x_{1}^{\prime}(t)\ ,\ \dots\ ,\ x_{N}^{\prime}(t)\right)$ is the derivative of $\mathbf{x}$ with respect to $t.$

We state several properties pertaining to the lead matrix. Firstly, the matrix $\mathbf{A}$ is skew-symmetric i.e. $\mathbf{A}^{\text{T}}=-\mathbf{A}.$ Equivalently, $A_{m,n}=-A_{n,m}$ for all $1\leq m,n\leq N,$ which follows immediately upon swapping the indices $m$ and $n$ in the equation (1.16). Next, each eigenvalue of $\mathbf{A}$ is either $0$ or purely imaginary [11, Exercise 7B.4]. These purely imaginary eigenvalues of $\mathbf{A}$ come in complex conjugate pairs: if $i\lambda$ is an eigenvalue of $\mathbf{A}$ for some nonzero $\lambda\in\mathbb{R},$ then $-i\lambda$ is also an eigenvalue of $\mathbf{A}.$ Lastly, $\mathbf{A}$ is diagonalizable [11, Exercise 7B.4]. In particular, there is an orthonormal basis of $\mathbb{C}^{N}$ consisting of the eigenvectors of $\mathbf{A}.$

For each $1\leq n\leq N,$ we let $\lambda_{n}(\mathbf{A})$ be the $n$ -th largest eigenvalue of $\mathbf{A}$ in modulus and consider the corresponding orthonormal basis of eigenvectors $\left\{\mathbf{v}_{n}(\mathbf{A})\right\}_{n=1}^{N},$ in which $\mathbf{v}_{n}(\mathbf{A})$ is the associated unit eigenvector with $\lambda_{n}(\mathbf{A})$ for each $1\leq n\leq N.$ We refer to the eigenvector $\mathbf{v}_{1}(\mathbf{A})$ as the leading (dominant) eigenvector of $\mathbf{A}.$ If $\lambda_{1}(\mathbf{A})$ is nonzero, then our notation for the largest eigenvalue and leading eigenvector is actually ambiguous because $\lambda_{2}(\mathbf{A})$ would be the complex conjugate of $\lambda_{1}(\mathbf{A})$ and thus satisfy $|\lambda_{1}(\mathbf{A})|=|\lambda_{2}(\mathbf{A})|.$ For our purposes, however, this ambiguity does not matter. By [11, Exercise 7B.4], we have $\mathbf{v}_{2}(\mathbf{A})=\overline{\mathbf{v}_{1}(\mathbf{A})},$ in which $\overline{\mathbf{v}_{1}(\mathbf{A})}$ is the vector whose $n$ -th component is equal to the complex conjugate of the $n$ -th component of $\mathbf{v}_{1}(\mathbf{A}).$ Therefore, regardless of whether $\lambda_{1}(\mathbf{A})$ is assigned to have a positive or negative imaginary part, the leading eigenvector $\mathbf{v}_{1}(\mathbf{A})$ is unique up to the complex conjugation of its components.

Finally, we mention how to construct a lead matrix so that it is suitable for practical applications. In the real world, one does not directly observe a continuous signal. Rather, one observes a signal at finitely many different times. Given a signal $\mathbf{x}:T\rightarrow\mathbb{R}^{N}$ with $T=[a,b],$ a time series is a finite collection of the form $\left\{\mathbf{x}_{t_{k}}\right\}_{k=1}^{K},$ where $K\in\mathbb{N}$ is fixed, $\left\{t_{k}\right\}_{k=1}^{K}$ is an increasing sequence in $T$ with $t_{1}=a$ and $t_{K}=b,$ and $\mathbf{x}_{t_{k}}=\mathbf{x}(t_{k}).$ Typically, $K$ is a very large integer. We define the lead matrix for the time series $\left\{\mathbf{x}_{t_{k}}\right\}_{k=1}^{K}$ to be the lead matrix $\mathbf{A}$ for the signal $\widetilde{\mathbf{x}}:T\rightarrow\mathbb{R}^{N}$ parameterizing the oriented curve formed by concatenating the oriented line segments $L_{1}\ ,\ \dots\ ,\ L_{K},$ in which $L_{k}$ has starting point $\mathbf{x}_{t_{k}}$ and ending point $\mathbf{x}_{t_{k+1}}$ for each $1\leq k<K$ and $L_{K}$ has starting point $\mathbf{x}_{t_{K}}$ and ending point $\mathbf{x}_{t_{1}}.$ Explicitly,

\displaystyle\mathbf{A}

\displaystyle=\frac{1}{2}\sum_{k=1}^{K-1}\mathbf{x}_{t_{k}}\mathbf{x}_{t_{k+1}% }^{\text{T}}-\mathbf{x}_{t_{k+1}}\mathbf{x}_{t_{k}}^{\text{T}}.

(1.22)

If we write out $\mathbf{x}_{t_{k}}=(x_{t_{k},1}\ ,\ \dots\ ,\ x_{t_{k},N}),$ then the $(m,n)$ -th entry of $\mathbf{A}$ is the oriented area of the convex polygon with vertices $\left(x_{t_{1},m}\ ,\ x_{t_{1},n}\right)\ ,\ \dots\ ,\ \left(x_{t_{K},m}\ ,\ x% _{t_{K},n}\right)$ computed via the Shoelace Formula [12].

1.4.4 Chain of Offsets Model

Now, we assume that $\mathbf{x}:\mathbb{R}\rightarrow\mathbb{R}^{N}$ is a cyclic signal. We recall our second leader-follower dynamics question that we posed in Section 1.1, which is to determine the order in which the $N$ component signals of $\mathbf{x}$ evolve throughout time assuming all these signals trace some underlying signal up to scaling constants and time shifts. Formally, we assume that the cyclic signal $\mathbf{x}$ satisfies the following Chain of Offsets Model (COOM) [3], which stipulates the existence of an underlying real-valued $P$ -periodic signal $\phi,$ scaling constants $c_{1}\ ,\ \dots\ ,\ c_{N}>0,$ offsets $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{N}\in\mathbb{R}/(P\mathbb{Z}),$ and a time reparameterization $r:\mathbb{R}\rightarrow\mathbb{R}$ such that

\displaystyle x_{n}(t)

\displaystyle=c_{n}\phi(r(t)-\alpha_{n}).

(1.23)

for each $1\leq n\leq N$ and each $t\in\mathbb{R}.$

Given the cyclic signal $\mathbf{x}:\mathbb{R}\rightarrow\mathbb{R}^{N},$ we determine the cyclic order of the offsets under COOM. But this time, we will utilize the lead matrix to do so. Without loss of generality, we may assume that $r$ is the identity map, which reduces the COOM in (1.23) to the COOM defined earlier in (1.3). If $r$ is not the identity map to begin with, then we can the consider the COOM in (1.3) for the $N$ -dimensional signal $\widetilde{\mathbf{x}}=\mathbf{x}\circ r^{-1}.$

As before, we begin by writing the underlying signal $\phi$ as the Fourier series in (1.5). Let $\mathbf{A}$ be the lead matrix for $\mathbf{x}$ over one period under COOM, that is, the lead matrix for the restricted signal $\mathbf{x}|_{[0,P]}.$ Then, for each $1\leq m,n\leq N,$ the $(m,n)$ -th entry of $\mathbf{A}$ is explicitly

$\displaystyle A_{m,n}$	$\displaystyle=\frac{1}{2}\int_{0}^{P}x_{m}(t)\ x_{n}^{\prime}(t)-x_{n}(t)\ x_{% m}^{\prime}(t)\ dt$
	$\displaystyle=\frac{1}{2}\int_{0}^{P}\sum_{k,\ell\in\mathbb{Z}}c_{n}c_{m}\left% (\frac{2\pi i\ell}{P}\right)\hat{\phi}_{k}\ \hat{\phi}_{\ell}\ e^{\frac{2\pi ik% (t-\alpha_{m})}{P}}\ e^{\frac{2\pi i\ell(t-\alpha_{n})}{P}}\ dt$
	$\displaystyle\qquad-\frac{1}{2}\int_{0}^{P}\sum_{k,\ell\in\mathbb{Z}}c_{n}c_{m% }\left(\frac{2\pi ik}{P}\right)\hat{\phi}_{k}\ \hat{\phi}_{\ell}\ e^{\frac{2% \pi ik(t-\alpha_{m})}{P}}\ e^{\frac{2\pi i\ell(t-\alpha_{n})}{P}}\ dt$
	$\displaystyle=\pi i\ c_{n}c_{m}\ \sum_{k,\ell\in\mathbb{Z}}\left(\ell-k\right)% \ \delta_{k,-\ell}\ \hat{\phi}_{k}\ \hat{\phi}_{\ell}\ e^{-\frac{2\pi ik\alpha% _{m}+2\pi i\ell\alpha_{n}}{P}}$
	$\displaystyle=-2\pi i\ c_{n}c_{m}\ \sum_{k\in\mathbb{Z}}k\left\|\hat{\phi}_{k}% \right\|^{2}e^{\frac{2\pi ik\left(\alpha_{n}-\alpha_{m}\right)}{P}}$
	$\displaystyle=2\pi c_{m}c_{n}\sum_{k\in\mathbb{N}}k\left\|\hat{\phi}_{k}\right\|% ^{2}\sin\left(\frac{2\pi k(\alpha_{m}-\alpha_{n})}{P}\right).$	(1.24)

1.4.5 One Harmonic

Assume that the underlying $P$ -periodic signal $\phi$ in (1.23) has only one harmonic. This means the Fourier Coefficient $\hat{\phi}_{k}$ is nonzero for only one index $k\in\mathbb{N}.$ We show that we can recover the cyclic order of the offsets via the leading eigenvector of the lead matrix $\mathbf{A}$ defined in (1.24). Note for each $1\leq m,n\leq N,$ the $(m,n)$ -th entry of the lead matrix $\mathbf{A}$ is explicitly

	$\displaystyle A_{m,n}$	$\displaystyle=2\pi c_{m}c_{n}k\left\|\hat{\phi}_{k}\right\|^{2}\sin\left(\frac{2% \pi k(\alpha_{m}-\alpha_{n})}{P}\right)$
		$\displaystyle=2\pi c_{m}c_{n}k\left\|\hat{\phi}_{k}\right\|^{2}\left(\sin\left(% \frac{2\pi k\alpha_{m}}{P}\right)\ \cos\left(\frac{2\pi k\alpha_{n}}{P}\right)% -\cos\left(\frac{2\pi k\alpha_{m}}{P}\right)\ \sin\left(\frac{2\pi k\alpha_{n}% }{P}\right)\right).$		(1.25)

We show that if the lead matrix $\mathbf{A}$ is nonzero, then it must have rank two. Recall if $\mathbf{a},\mathbf{b}\in\mathbb{C}^{N},$ then the outer product of $\mathbf{a}$ and $\mathbf{b}$ is the complex $N\times N$ matrix $\mathbf{a}\otimes\mathbf{b}=\mathbf{a}\mathbf{b}^{*},$ where $\mathbf{b}^{*}=\overline{\mathbf{b}^{\text{T}}}$ is the conjugate transpose of $\mathbf{b}.$ A statement made in [3] is that a complex skew-symmetric matrix has rank two if and only if it is of the form $\mathbf{a}\otimes\mathbf{b}-\mathbf{b}\otimes\mathbf{a}$ for some linearly independent $\mathbf{a},\mathbf{b}\in\mathbb{C}^{N}.$ So we verify that our lead matrix $\mathbf{A}$ indeed satisfies this condition stated in [3]. Define the vectors

	$\displaystyle\mathbf{a}$	$\displaystyle=\sqrt{2\pi k}\|\hat{\phi}_{k}\|\left(c_{1}\cos\left(\frac{2\pi k% \alpha_{1}}{P}\right)\ ,\ \dots\ ,\ c_{N}\cos\left(\frac{2\pi k\alpha_{N}}{P}% \right)\right)$
	$\displaystyle\mathbf{b}$	$\displaystyle=\sqrt{2\pi k}\|\hat{\phi}_{k}\|\left(c_{1}\sin\left(\frac{2\pi k% \alpha_{1}}{P}\right)\ ,\ \dots\ ,\ c_{N}\sin\left(\frac{2\pi k\alpha_{N}}{P}% \right)\right).$

Then, by (1.25) immediately implies

\displaystyle\mathbf{A}

\displaystyle=\mathbf{a}\otimes\mathbf{b}-\mathbf{b}\otimes\mathbf{a}.

(1.26)

Now, we need to show that $\mathbf{a}$ and $\mathbf{b}$ are linearly independent in $\mathbb{C}^{N}.$ Consider the $N\times 2$ matrix $\mathbf{X}$ whose first row is $\mathbf{a}$ and second row is $\mathbf{b}.$ The condition that $\mathbf{a}$ and $\mathbf{b}$ are linearly independent is equivalent to the condition $\mathbf{X}$ has rank $2.$ The latter condition is equivalent to the condition there exists a $2\times 2$ submatrix of $\mathbf{X}$ with a nonzero determinant. Recall our assumption $\mathbf{A}$ is nonzero. This means $A_{m,n}$ has a nonzero entry for some two indices $1\leq m<n\leq N.$ On the other hand, consider the submatrix $\widetilde{\mathbf{X}}$ formed by extracting the $m$ -th and $n$ -th rows of $\mathbf{X}$ for the same choices of $m$ and $n,$ namely the $2\times 2$ matrix

\displaystyle\widetilde{\mathbf{X}}_{m,n}

\displaystyle=\begin{bmatrix}\sqrt{2\pi k}|\widehat{\phi}_{k}|c_{m}\cos\left(% \frac{2\pi k\alpha_{m}}{P}\right)&\sqrt{2\pi k}|\widehat{\phi}_{k}|c_{m}\sin% \left(\frac{2\pi k\alpha_{m}}{P}\right)\\ \sqrt{2\pi k}|\widehat{\phi}_{k}|c_{n}\cos\left(\frac{2\pi k\alpha_{n}}{P}% \right)&\sqrt{2\pi k}|\widehat{\phi}_{k}|c_{n}\sin\left(\frac{2\pi k\alpha_{n}% }{P}\right)\end{bmatrix}.

By (1.25), the determinant of $\widetilde{\mathbf{X}}_{m,n}$ is precisely $A_{m,n},$ which shows $\mathbf{A}$ has rank two.

We now determine the leading eigenvector of $\mathbf{A}.$ Since $\mathbf{A}$ has rank two, it has exactly two nonzero eigenvalues coming in a purely imaginary pair, namely $\lambda_{1}(\mathbf{A})$ and $\lambda_{2}(\mathbf{A}),$ in which we recall $\lambda_{n}(\mathbf{A})$ is the $n$ -th largest eigenvalue of $\mathbf{A}$ in modulus. By [3],

\displaystyle\lambda_{1}(\mathbf{A})

\displaystyle=i\sin(\theta)\|\mathbf{a}\|\|\mathbf{b}\|

is the largest eigenvalue of $\mathbf{A}$ with associated eigenvector

\displaystyle\mathbf{v}_{1}(\mathbf{A})

\displaystyle=-e^{-i\theta}\frac{\mathbf{a}}{\|\mathbf{a}\|}+\frac{\mathbf{b}}% {\|\mathbf{b}\|},

where $\theta$ is the angle between the vectors $\mathbf{a}$ and $\mathbf{b},$ and $\|\mathbf{a}\|$ is the Euclidean $2$ -norm of $\mathbf{a}.$ In particular, if $v_{1,n}(\mathbf{A})$ is the $n$ -th component of $\mathbf{v}_{1}(\mathbf{A}),$ then upon extracting real and imaginary parts of $v_{1,n}(\mathbf{A}),$ we obtain

	$\displaystyle\Re(v_{1,n}(\mathbf{A}))$	$\displaystyle=-\frac{a_{n}\cos(\theta)}{\\|\mathbf{a}\\|}+\frac{b_{n}}{\\|\mathbf% {b}\\|}$		(1.27)
	$\displaystyle\Im(v_{1,n}(\mathbf{A}))$	$\displaystyle=\frac{a_{n}\sin(\theta)}{\\|\mathbf{a}\\|},$		(1.28)

where

	$\displaystyle a_{n}$	$\displaystyle=\sqrt{2\pi k}\|\hat{\phi}_{k}\|c_{n}\cos\left(\frac{2\pi k\alpha_{% n}}{P}\right)$
	$\displaystyle b_{n}$	$\displaystyle=\sqrt{2\pi k}\|\hat{\phi}_{k}\|c_{n}\sin\left(\frac{2\pi k\alpha_{% n}}{P}\right)$

are the $n$ -th components of $\mathbf{a}$ and $\mathbf{b},$ respectively.

Through the equations (1.27) and (1.28), we see $\left(\Re(v_{1,n}(\mathbf{A}))\ ,\ \Im(v_{1,n}(\mathbf{A}))\right)$ is the result of a certain linear transformation evaluated at the point $\left(a_{n},b_{n}\right).$ Since the points $(a_{1},b_{1})\ ,\ \dots\ ,\ (a_{N},b_{N})$ lie on a circle, the components of the leading eigenvector $\mathbf{v}_{1}(\mathbf{A})$ lie on an ellipse. The cyclic order in which the points $\left\{(a_{n},b_{n})\right\}_{n=1}^{N}$ are traversed on the circle is the cyclic order in which the points $\left\{\left(\Re(v_{1,n}(\mathbf{A}))\ ,\ \Im(v_{1,n}(\mathbf{A}))\right)% \right\}_{n=1}^{N}$ are traversed on the ellipse. Thus, the cyclic order of $\left\{\text{Arg}\left(v_{1,n}(\mathbf{A})\right)\right\}_{n=1}^{N},$ where $\text{Arg}(z)$ is the principal argument of the complex number $z,$ matches the cyclic order of the offsets in COOM. We refer to the principal arguments of the components of the leading eigenvector as phases.

We show an example of how the cyclic order of the leading eigenvector component phases and the cyclic order of the offsets match when the underlying signal $\phi$ has only one harmonic. Reconsider the sinusoidal signal (1.1) for each $1\leq n\leq N.$ On one hand, we already established the identity permutation in $S_{N}$ represents the cyclic order of the offsets. On the other hand, in Figure 1.7, for $N=10,$ we plot the heat map of the lead matrix $\mathbf{A}$ over one period, as well as the complex components, component moduli, and phases of the leading eigenvector $\mathbf{v}_{1}(\mathbf{A}).$ Since the components of $\mathbf{v}_{1}(\mathbf{A})$ have equal moduli, they lie on a circle. Observe the permutation $\sigma\in S_{10}$ that sorts the leading eigenvector component phases in ascending order is given by $\sigma(n)=n+1,$ in which $n+1$ is indexed mod $10.$ This is the identity permutation up to a cyclic shift.

1.4.6 Multiple Harmonics

Suppose the underlying $P$ -periodic signal $\phi$ in (1.23) has multiple harmonics, meaning $\phi$ has multiple nonzero Fourier coefficients. Assume, however, one of these Fourier coefficients dominates: there exists exactly one $k\in\mathbb{N}$ such that

\displaystyle|\hat{\phi}_{k}|^{2}\gg\sum_{\ell\neq k}|\hat{\phi}_{\ell}|^{2}.

With these assumptions, it is still possible to approximately recover the cyclic order of the offsets $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{N}$ under COOM via the leading eigenvector of the lead matrix over one period.

Let $\widetilde{\mathbf{A}}_{k}$ be the rank two skew-symmetric $N\times N$ matrix, whose $(m,n)$ -th entry is equal to (1.25) for all $1\leq m,n\leq N$ . Then, the $(m,n)$ -th entry of $\widetilde{\mathbf{A}}_{k}$ will dominate the series expansion in (1.24). This means we well approximate the original lead matrix $\mathbf{A}$ with the lower rank matrix $\widetilde{\mathbf{A}}_{k}.$

On the other hand, a rank two approximation of $\mathbf{A}$ is the matrix $\mathbf{P}\mathbf{A}\mathbf{P},$ in which $\mathbf{P}$ is the $N\times N$ matrix corresponding to the projection operator onto the subspace of $\mathbb{C}^{N}$ spanned by the vectors $\mathbf{v}_{1}(\mathbf{A})$ and $\mathbf{v}_{2}(\mathbf{A}).$ The smaller the Frobenius norm $\|\mathbf{P}\mathbf{A}\mathbf{P}-\mathbf{A}\|_{F}$ is, the better this low rank approximation is of $\mathbf{A}$ . Moreover, since $\mathbf{P}$ sends $\mathbf{v}_{1}(\mathbf{A})$ to itself, we have

	$\displaystyle(\mathbf{P}\mathbf{A}\mathbf{P})\ \mathbf{v}_{1}(\mathbf{A})$	$\displaystyle=(\mathbf{P}\mathbf{A})\ \mathbf{v}_{1}(\mathbf{A})$
		$\displaystyle=\lambda_{1}(\mathbf{A})\ \mathbf{P}\ \mathbf{v}_{1}(\mathbf{A})$
		$\displaystyle=\lambda_{1}(\mathbf{A})\ \mathbf{v}_{1}(\mathbf{A}),$

which means the leading eigenvector of $\mathbf{P}\mathbf{A}\mathbf{P}$ is $\mathbf{v}_{1}(\mathbf{A}).$ Therefore, assuming the low rank matrix $\mathbf{P}\mathbf{A}\mathbf{P}$ approximates $\mathbf{A}$ well enough, we can still approximately recover the cyclic order of the offsets $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{N}$ in COOM via the cyclic order of the component phases of $\mathbf{v}_{1}(\mathbf{A}).$

For practical purposes, one heuristic indicator [3] that measures the dominance of the low rank $2$ approximation of the lead matrix $\mathbf{A}$ is the magnitude of $\left|\frac{\lambda_{1}(\mathbf{A})}{\lambda_{3}(\mathbf{A})}\right|,$ the ratio of the first to third largest eigenvalue of $\mathbf{A}.$ This is because in the rank $2$ situation, this ratio would be infinite.

1.5 Ornstein-Uhlenbeck Process

In this section, we fix a probability space $(\Omega,\mathcal{A},\mathbb{P})$ and $N\in\mathbb{N}.$ We consider an ( $N$ -dimensional) Ornstein Uhlenbeck (OU) process on this probability space, which is an $N$ -dimensional stochastic process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ satisfying the linear stochastic differential equation [13]:

\displaystyle d\mathbf{x}(t)

\displaystyle=-\mathbf{B}\ \mathbf{x}(t)\ dt+\boldsymbol{\Sigma}\ d\mathbf{w}(% t).

(1.29)

Here, $\mathbf{B}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is a fixed linear operator such that $-\mathbf{B}$ is Hurwitz (stable) i.e. all eigenvalues of $-\mathbf{B}$ have negative real parts, $\boldsymbol{\Sigma}:\mathbb{R}^{M}\rightarrow\mathbb{R}^{N}$ is a fixed linear operator for some $M\in\mathbb{N},$ and $\left\{\mathbf{w}(t)\right\}_{t\geq 0}$ is the standard $M$ -dimensional Wiener process independent of the initial random vector $\mathbf{x}(0)$ . We refer to $\mathbf{B}$ as the friction (drift) operator and $\boldsymbol{\Sigma}$ as the volatility operator, and we collectively refer to $\mathbf{B}$ and $\boldsymbol{\Sigma}$ as the OU model parameters. We refer to the matrix

\displaystyle\mathbf{D}

\displaystyle=\frac{\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\text{T}}}{2}

(1.30)

as the diffusion operator.

The OU process is an ergodic Gaussian diffusion process [14] that prominently appears in Financial Mathematics and Biology [15, 16, 17]. In engineering disciplines, the OU Process is heuristically characterized as the solution to the Langevin equation [18]:

\displaystyle\dot{\mathbf{x}}(t)

\displaystyle=-\mathbf{B}\ \mathbf{x}(t)+\boldsymbol{\eta}(t),

(1.31)

in which $\left\{\boldsymbol{\eta}(t)\right\}_{t\geq 0}$ is the $N$ -dimensional Gaussian “white noise” process whose autocovariance function satisfies

\displaystyle\boldsymbol{\Gamma}_{\boldsymbol{\eta}}(s,t)

\displaystyle=2\mathbf{D}\ \delta(s-t),

for all $s,t\geq 0.$ Here, $\delta(\cdot)$ is the Dirac measure. The OU process has a unique stationary distribution [14, 19], namely the Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S},$ in which $\mathbf{S}$ satisfies the Lyapunov equation [20]:

\displaystyle\mathbf{B}\mathbf{S}+\mathbf{S}\mathbf{B}^{\text{T}}

\displaystyle=2\mathbf{D}.

(1.32)

For this reason, we refer to $\mathbf{S}$ as the stationary covariance operator. Any realization of the OU process will fluctuate around the origin $\mathbf{0}$ in the long term [18]. These long term fluctuations vary according to the stationary Gaussian distribution. Finally, if $\mathbf{x}(0)$ is the Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S},$ then the OU process is stationary.

1.6 Problem Statement

In this thesis, we will investigate Cyclicity Analysis for the stationary OU process with model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ in which $\mathbf{x}(0)$ is the Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S}.$ We define the auxiliary lead process, which is the $N\times N$ matrix-valued stochastic process $\left\{\mathbf{A}(t)\right\}_{t\geq 0},$ where

\displaystyle\mathbf{A}(t)

\displaystyle=\frac{1}{2}\int_{0}^{t}\mathbf{x}(s)\ d\mathbf{x}^{\text{T}}(s)-% d\mathbf{x}(s)\ \mathbf{x}^{\text{T}}(s)

(1.33)

for each $t\geq 0.$ The right hand side is a random matrix whose entries are one-dimensional Ito stochastic integrals [21] with respect to the component OU processes $\left\{x_{n}(t)\right\}_{t\geq 0}$ . In particular, the $(m,n)$ -th entry of $\mathbf{A}(t)$ is a one-dimensional Ito integral of the form

\displaystyle A_{m,n}(t)

\displaystyle=\frac{1}{2}\int_{0}^{t}x_{m}(s)\ dx_{n}(s)-x_{n}(s)\ dx_{m}(s).

By definition, $\mathbf{A}(t)$ represents the lead matrix corresponding to a realization of the OU process over the interval $[0,t].$

We will prove the lead process obeys the strong law of large numbers identity:

\displaystyle\lim_{t\rightarrow\infty}\frac{\mathbf{A}(t)}{t}

\displaystyle=\frac{\mathbf{B}\mathbf{S}-\mathbf{S}\mathbf{B}^{\text{T}}}{2},

(1.34)

almost surely. The identity (1.34) states we can infer the long term average pairwise leader follower dynamics amongst the $N$ component OU process through a time-averaged lead matrix corresponding to one single realization of the OU process. For this reason, we refer to the matrix

\displaystyle\mathbf{Q}

\displaystyle=\frac{\mathbf{B}\mathbf{S}-\mathbf{S}\mathbf{B}^{\text{T}}}{2}

(1.35)

as the lead matrix of the OU process.

Next, we consider a cyclic OU process, in which the friction matrix $\mathbf{B}$ is specifically a circulant matrix: if $(b_{1}\ ,\ \dots\ ,\ b_{N})\in\mathbb{R}^{N}$ is the first row of $\mathbf{B},$ then the $(m,n)$ -th entry of $\mathbf{B}$ is $b_{n-m+1},$ where $n-m+1$ is indexed mod $N.$ We assume $b_{2}\ ,\ \dots\ ,\ b_{N}\leq 0$ and $b_{1}>0$ is large enough so that $\mathbf{B}$ is a valid friction matrix. We consider a signal propagation model governed by this cyclic OU process. In our signal propagation model, we assume that there is a network of $N$ sensors located in space and an underlying one-dimensional signal propagating throughout the network. Similar signal propagation models are considered in [22, 4, 23], for example. Letting $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ be the OU process, we interpret the $n$ -th component of $\mathbf{x}(t)$ as a measurement made by the $n$ -th sensor of the propagating effect at time $t\geq 0.$ The friction matrix $\mathbf{B}$ represents the cyclic sensor network structure. We refer to the constant $b_{p}$ as a propagation coefficient for each $p>1,$ and we interpret its magnitude to represent how receptive the $n$ -th sensor is to the activity within the neighboring $(n+p-1)$ -th sensor for each $1\leq n\leq N,$ in which $n+p-1$ is indexed mod $N.$ We refer to the constant $b_{1}$ as a suppression coefficient, and we interpret its magnitude to represent how quickly the signal dissipates as it is propagating throughout the network. Meanwhile, the diffusion matrix $\mathbf{D}=\frac{\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\text{T}}}{2}$ is the covariance matrix whose $(m,n)$ -th entry is the covariance of the $m$ -th and $n$ -th component noises injected into the $m$ -th and $n$ -th sensors, respectively.

For large dimension $N$ and for different choices of cyclic OU model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ we derive the lead matrix $\mathbf{Q}$ of the cyclic OU process. Then, we investigate whether Cyclicity Analysis enables us to recover the network structure given by $\mathbf{B}.$ More specifically, we investigate whether the structure of the leading eigenvector $\mathbf{v}_{1}(\mathbf{Q})$ recovers the structure of the network.

For example, let $\mathbf{B}$ be the $N\times N$ circulant friction matrix such that the suppression coefficient $b_{1}>1,$ the propagation coefficient $b_{2}=-1,$ and all other propagation coefficients are $0.$ The cyclic network structure of $\mathbf{B}$ is such that the $n$ -th sensor is only linked to the $(n+1)$ -th sensor, the neighboring sensor to its right. Does the structure of the leading eigenvector $\mathbf{v}_{1}(\mathbf{Q})$ reflect this ground truth network structure in any way ? For example, does the (cyclic) order of the leading eigenvector component phases match the expected (cyclic) order of the sensors receiving the signal as it propagates throughout the network ? In addition, do the leading eigenvector component moduli and the lead matrix eigenvalues indicate anything about the network structure ?

We will investigate our problem statement under two specific regimes: the first regime is when the suppression coefficient $b_{1}$ is large i.e. the propagating signal quickly dissipates, and the second regime is when $b_{1}$ is small i.e. the propagating signal does not quickly dissipate. In the former regime, the circulant friction matrix $\mathbf{B}$ is perturbed significantly from an unstable matrix, while in the latter regime, $\mathbf{B}$ is perturbed slightly away from an unstable matrix. We investigate these regimes because they yield different types of asymptotic structures pertaining to the lead matrix $\mathbf{Q}$ and its leading eigenvector $\mathbf{v}_{1}(\mathbf{Q}).$

1.7 Literature Review

We conclude the introduction by conducting a literature review of the main topics of the thesis.

1.7.1 Lead Lag Dynamics

There are various tools to identify lead lag dynamics between two signals or between two time series. Such tools have been utilized to study financial markets [24] and fMRI signals [25]. We already discussed one tool, namely the cross-correlation function. Another algorithmic tool is Dynamic Time Warping (DTW) [26]. We briefly describe it. Consider two one-dimensional time series of the form $\left\{x_{m}\right\}_{m=1}^{M}$ and $\left\{y_{m}\right\}_{n=1}^{N}.$ The goal of DTW is to find the optimal alignment of such time-series. To this end, one constructs an $M\times N$ matrix $\mathbf{C}$ whose $(j,k)$ -th entry corresponds to the optimal cost between the restricted time series $\left\{x_{m}\right\}_{m=1}^{j}$ and $\left\{y_{m}\right\}_{n=1}^{k}.$ The optimal cost is defined via the formula

\displaystyle C_{j,k}

\displaystyle=\begin{cases}d(x_{j},y_{k})&j=k=1\\ d(x_{j},y_{k})+\min\left(C_{j-1,k}\ ,\ C_{j-1,k-1}\ ,\ C_{j,k-1}\right)&\text{% else},\end{cases}

(1.36)

where $d(\cdot,\cdot)$ is a specified distance metric. Then, a warping path is any collection of indices $\left\{(j_{1},k_{1})\ ,\ \dots\ ,\ (j_{P},k_{P})\right\}$ satisfying $(j_{1},k_{1})=(1,1)$ and $(j_{P},k_{P})=(M,N)$ and $0\leq j_{p+1}-j_{p}\ ,\ k_{p+1}-k_{p}\leq 1$ for all $1\leq p<P,$ where $P\geq 2$ is a fixed length. Finally, the optimal alignment between the original time-series $\left\{x_{m}\right\}_{m=1}^{M}$ and $\left\{y_{m}\right\}_{n=1}^{N}$ is the warping path $\left\{(j_{1},k_{1})\ ,\ \dots\ ,\ (j_{P},k_{P})\right\}$ minimizing the quantity

\displaystyle\sum_{p=1}^{P}C_{j_{p},k_{p}}^{2}.

DTW, however, has drawbacks. Seeking the warping path to guarantee optimal alignment is a computationally expensive problem. Moreover, DTW is prone to noise.

1.7.2 Cyclicity Analysis

We review the literature on Cyclicity Analysis. As mentioned earlier, Cyclicity Analysis was officially introduced in [3] as a means to study cyclic signals. Later on, a dissertation [23] was published demonstrating how Cyclicity Analysis of fMRI signals yields new insights into understanding the brain processes of patients with a neurological condition known as tinnitus. Moreover, the authors in [4] published an experimental article comparing the results given by Cyclicity Analysis on fMRI signals to results given by other standard tools.

We emphasize that in this thesis, we are applying Cyclicity Analysis in a completely different setting. We are using it to investigate signals generated via a stochastic process.

1.7.3 OU Process

As mentioned earlier, the OU process is one of the most common stochastic processes with applications in Financial Mathematics and Biology.

Firstly, many authors have studied the estimation of OU process model parameters. In [27], for example, the authors showed using MAP estimation how to recover OU model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ in which $\boldsymbol{\Sigma}$ is an $N\times N$ invertible matrix, given a time series that approximates a realization of an OU process. In [14], the authors showed how to approximate the stationary distribution of a generalized OU process with a modified drift coefficient. We emphasize that our problem statement is completely different from model parameter estimation. We are not trying to recover the friction matrix $\mathbf{B}$ itself from realizations of the OU process; rather, we interpret $\mathbf{B}$ as a weighted adjacency matrix associated with a network of sensors, and we are investigating whether we can recover the network structure induced by $\mathbf{B}$ from any realization of the OU process.

Moreover, the lead matrix $\mathbf{Q}$ of the OU process defined in (3.40) has been studied in the literature as a means to characterize the time-irreversibility of the OU process [18, 28, 29]. A stationary stochastic process $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ is said to be time-reversible [29] if for all $K\in\mathbb{N}$ for all times $0\leq t_{1}\ ,\ \dots\ ,\ t_{K}\leq\tau,$ the probability distributions of $\left(\mathbf{y}(t_{1})\ ,\ \dots\ ,\ \mathbf{y}(t_{K})\right)$ and of $\left(\mathbf{y}(\tau-t_{1})\ ,\ \dots\ ,\ \mathbf{y}(\tau-t_{K})\right)$ are equal. The main result is the stationary OU process is time reversible if and only if $\mathbf{B}\mathbf{S}$ is symmetric, where $\mathbf{S}$ is the stationary covariance matrix of the OU process. Equivalently, the OU process is time reversible if and only if the matrix $\mathbf{Q}$ vanishes.

To study the degree to which the OU process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is time irreversible, we briefly outline the approach described in [28]. Let $P:\mathbb{R}^{N}\times[0,\infty)\rightarrow\mathbb{R}$ be the probability density function of $\mathbf{x}(t)$ at time $t\geq 0.$ Since we are considering the OU stationary process, for each $t\geq 0,$ the random vector $\mathbf{x}(t)$ has a Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S}.$ Explicitly, we have

\displaystyle P(\mathbf{x},t)

\displaystyle=\frac{1}{\sqrt{(2\pi)^{N}|\det(\mathbf{S})|}}\ e^{-\frac{\mathbf% {x}^{\text{T}}\mathbf{S}^{-1}\mathbf{x}}{2}}

(1.37)

for each $t\geq 0.$ The density function $P$ satisfies the Fokker-Planck (Kolmogorov forward) equation [28], which is the parabolic partial differential equation

\displaystyle\frac{\partial P}{\partial t}

\displaystyle=\nabla\cdot\left(\mathbf{B}\ \mathbf{x}\ P+\mathbf{D}\ \nabla P% \right),

(1.38)

where $\mathbf{D}$ is the diffusion matrix, $\nabla P$ is the gradient of $P$ with respect to variable $\mathbf{x}=(x_{1}\ ,\ \dots\ ,\ x_{N}),$ and $\nabla\cdot\mathbf{F}$ is the divergence of the vector field $\mathbf{F}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}.$ Define the map $\mathbf{J}:\mathbb{R}^{N}\times[0,\infty)\rightarrow\mathbb{R}^{N}$ by the expression within the divergence operator on the right hand side of (1.38), namely

\displaystyle\mathbf{J}(\mathbf{x},t)

\displaystyle=-\mathbf{B}\ \mathbf{x}\ P(\mathbf{x},t)-\mathbf{D}\ \nabla P(% \mathbf{x},t).

(1.39)

The map $\mathbf{J}$ is known as the probability flux (current) [18, 28]. In general, the probability flux describes how probabilities associated with different states in the state space of a stochastic process change over time. Intuitively, when the probability flux is nonzero, the realizations follow paths that lead towards regions with higher probability density i.e. those states with higher probability of occurrence.

For each $t\geq 0,$ upon recognizing that $\nabla P(\mathbf{x},t)=-\mathbf{S}^{-1}\ \mathbf{x}\ P(\mathbf{x},t),$ we can rewrite the probability flux in the following way:

	$\displaystyle\mathbf{J}(\mathbf{x},t)$	$\displaystyle=-(\mathbf{B}\mathbf{x}\ P(\mathbf{x},t)+\mathbf{D}\mathbf{S}^{-1% }\mathbf{x}\ P(\mathbf{x},t))$
		$\displaystyle=\left(\mathbf{D}\mathbf{S}^{-1}-\mathbf{B}\right)\mathbf{x}\ P(% \mathbf{x},t)$
		$\displaystyle=\left(\frac{\mathbf{B}\mathbf{S}+\mathbf{S}\mathbf{B}^{\text{T}}% }{2}\ \mathbf{S}^{-1}-\mathbf{B}\right)\mathbf{x}\ P(\mathbf{x},t)$
		$\displaystyle=\left(\frac{\mathbf{S}\mathbf{B}^{\text{T}}\mathbf{S}^{-1}-% \mathbf{B}}{2}\right)\mathbf{x}\ P(\mathbf{x},t)$
		$\displaystyle=\left(\frac{\left(\mathbf{B}\mathbf{S}-\mathbf{2}\mathbf{Q}% \right)\mathbf{S}^{-1}-\mathbf{B}}{2}\right)\mathbf{x}\ P(\mathbf{x},t)$
		$\displaystyle=-\mathbf{Q}\mathbf{S}^{-1}\mathbf{x}\ P(\mathbf{x},t).$

Thus, the stationary OU process is time-reversible if and only if the probability flux $\mathbf{J}$ vanishes. If the flux does not vanish, then observe $\mathbf{Q}$ is involved in the measurement of the flux.

Chapter 2 Preliminaries

In this chapter, we will briefly review the needed notations, definitions, and standard results from Linear Algebra, Probability Theory, and Stochastic Processes that will be used throughout this thesis.

2.1 Linear Algebra

Throughout this section, we let $\mathbb{F}\in\left\{\mathbb{R},\mathbb{C}\right\}$ and we fix $M,N\in\mathbb{N}.$ We denote $\text{Mat}_{M,N}(\mathbb{F})$ as the set of all $M\times N$ matrices with entries in $\mathbb{F},$ and we tacitly identify the $N$ -dimensional vector space $\mathbb{F}^{N}$ with $\text{Mat}_{N,1}(\mathbb{F}).$ Recall if $\mathbf{A}\in\text{Mat}_{M,N}(\mathbb{F}),$ then there is a unique linear operator between $\mathbb{F}^{N}$ and $\mathbb{F}^{M}$ whose transformation matrix is $\mathbf{A}.$ Therefore, we will interchangeably refer to $\mathbf{A}$ as either a matrix or as a linear operator between the appropriate vector spaces.

Let $\mathbf{A}\in\text{Mat}_{M,N}(\mathbb{F})$ and $\mathbf{x}\in\mathbb{F}^{N}.$ We recall some standard linear algebra notation. For each $1\leq m\leq M$ and each $1\leq n\leq N,$ we denote $A_{m,n}$ as the $(m,n)$ -th entry of $\mathbf{A},$ and we denote $x_{n}$ as the $n$ -th component of $\mathbf{x}.$ We denote $\mathbf{A}^{\text{T}}$ as the transpose of $\mathbf{A}$ and denote $\mathbf{A}^{*}=\overline{\mathbf{A}^{\text{T}}}$ as the conjugate-transpose of $\mathbf{A}.$ We denote $\mathbf{I}$ as the $N\times N$ identity matrix, in which the value of $N$ will be clear in context. We denote $\mathbf{0}$ as either the $M\times N$ zero matrix or the zero vector in $\mathbb{F}^{N},$ in which the situation, along with the values of $M$ and $N$ , will be clear in context. Finally, we denote $\mathbf{e}_{n}\in\mathbb{R}^{N}$ as the $n$ -th canonical basis vector having $n$ -th component equal to $1$ and all other components equal to $0.$

Now, let $\mathbf{A}\in\text{Mat}_{N,N}(\mathbb{F}).$ We denote $\det(\mathbf{A}),\text{tr}(\mathbf{A}),$ and $\mathbf{A}^{-1}$ as the determinant, trace, and inverse of $\mathbf{A},$ respectively. Recall the matrix exponential of $\mathbf{A}$ is the power series

\displaystyle\exp(\mathbf{A})

\displaystyle=\sum_{k\geq 0}\frac{\mathbf{A}^{k}}{k!}.

(2.1)

We state the basic properties of the matrix exponential that we will be using later in the thesis. Proofs of these properties can be found in [30], for example.

Theorem 2.1 (Matrix Exponential Properties).

If $\mathbf{A}\in\text{Mat}_{N,N}(\mathbb{R}),$ then the following properties hold.

i)

$\exp(\mathbf{A})$ commutes with $\mathbf{A}$ under matrix multiplication.
ii)

$\exp(\mathbf{A})$ is invertible with $(\exp(\mathbf{A}))^{-1}=\exp(\mathbf{-}\mathbf{A}).$

iii)

We have

\displaystyle\left(\exp{\mathbf{A}}\right)^{\text{T}}=\exp\left(\mathbf{A}^{% \text{T}}\right).

(2.2)

iv)

The matrix-valued map $t\mapsto\exp(t\mathbf{A})$ is differentiable on $\mathbb{R}$ . In particular,

\displaystyle(\exp{t\mathbf{A}})^{\prime}

\displaystyle=\mathbf{A}\ \exp(t\mathbf{A})=\exp(t\mathbf{A})\ \mathbf{A}.

(2.3)

Let $\mathbf{A}\in\text{Mat}_{N,N}(\mathbb{F}).$ The matrix $\mathbf{A}$ is Hurwitz if all of its eigenvalues have negative real parts. Note if $\mathbf{A}$ is Hurwitz, then $\lim_{t\rightarrow\infty}\exp(t\mathbf{A})=\mathbf{0}.$ The matrix $\mathbf{A}$ is Hermitian if $\mathbf{A}^{*}=\mathbf{A}$ and skew-Hermitian (anti-Hermitian) if $\mathbf{A}^{*}=-\mathbf{A}.$ Note that Hermitian real matrices are the same as symmetric real matrices, while skew-Hermitian real matrices are the same as skew-symmetric real matrices. We state the important properties of Hermitian and skew-Hermitian matrices. One may consult [30] for proofs.

Theorem 2.2 (Hermitian and skew-Hermitian Matrix Properties).

If $\mathbf{A}\in\text{Mat}_{N,N}(\mathbb{F})$ is either Hermitian or skew-Hermitian, then the following properties hold

i)

$\mathbf{A}$ is skew-Hermitian if and only if $i\mathbf{A}$ is Hermitian.
ii)

If $\mathbf{A}$ is Hermitian or skew-Hermitian, then $\mathbf{A}$ is diagonalizable; there is an orthonormal basis of $\mathbb{C}^{N}$ consisting of the eigenvectors of $\mathbf{A}.$
iii)

If $\mathbf{A}$ is Hermitian, then all its eigenvalues are real. If $\mathbf{A}$ is skew-Hermitian, then all eigenvalues of $\mathbf{A}$ are either $0$ or purely imaginary.

Finally, let $\mathbf{A}\in\text{Mat}_{N,N}(\mathbb{F}).$ The matrix $\mathbf{A}$ is circulant if there exists a sequence of numbers $\left\{a_{n}\right\}_{1\leq n<N}$ such that $A_{m,n}=a_{m-n+1},$ where $m-n+1$ is indexed mod $N.$ Circulant matrices obey several properties [31]. Circulant matrices form a subspace of $\text{Mat}_{N,N}(\mathbb{F}).$ They commute under matrix multiplication. The product of circulant matrices is circulant, and the transpose of a circulant matrix is circulant.

We state and prove the properties pertaining to the eigenvalues and eigenvectors of circulant matrices that are essential for this thesis.

Theorem 2.3 (Circulant Matrix Properties).

If $\mathbf{A},\mathbf{B}\in\text{Mat}_{N,N}(\mathbb{F})$ are circulant matrices, then the following properties hold.

$\mathbf{A}$ is diagonalizable. In particular, if the first row of $\mathbf{A}$ is of the form $(a_{1}\ ,\ \dots\ ,\ a_{N}),$ then for each $1\leq n\leq N,$ the vector

\displaystyle\mathbf{w}_{n}

\displaystyle=\sqrt{\frac{1}{N}}\left(1\ ,\ \omega_{n}\ ,\ \dots\ ,\ \omega_{n% }^{N-1}\right),

(2.4)

is an eigenvector of $\mathbf{A},$ associated with the eigenvalue

\displaystyle\mu_{n}

\displaystyle=\sum_{p=1}^{N}a_{p}\ \omega_{n}^{p-1},

(2.5)

where $\omega_{n}=e^{\frac{2\pi in}{N}}$ is an $N$ -th root of unity. Furthermore, the collection of eigenvectors $\left\{\mathbf{w}_{n}\right\}_{n=1}^{N}$ forms an orthonormal basis of $\mathbb{C}^{N},$ and we have $\mathbf{w}_{n}^{*}=\mathbf{w}_{N-n}^{T}$ for each $1\leq n\leq N.$

ii)

For each $k\geq 0,$ we have

\displaystyle\mathbf{A}^{k}

\displaystyle=\sum_{n=1}^{N}\mu_{n}^{k}\ \mathbf{w}_{n}\mathbf{w}_{N-n}^{\text% {T}}.

(2.6)

iii)

For each $t\in\mathbb{R},$ the matrix exponential of $t\mathbf{A}$ is explicitly

\displaystyle\exp(t\mathbf{A})

\displaystyle=\sum_{n=1}^{N}e^{\mu_{n}t}\ \mathbf{w}_{n}\mathbf{w}_{N-n}^{% \text{T}}.

(2.7)

Moreover, $\exp(t\mathbf{A})$ is a circulant matrix.

Proof.

We prove the first statement. For each $1\leq j,n\leq N,$ the $j$ -th component of $\mathbf{A}\mathbf{w}_{n}$ is of the form

	$\displaystyle\sqrt{\frac{1}{N}}\left(\sum_{p=1}^{N}a_{p-j+1}\ \omega_{n}^{p-1}\right)$	$\displaystyle=\sqrt{\frac{1}{N}}\left(\sum_{p=1}^{N-j+1}a_{p}\ \omega_{n}^{p+j% -2}+\sum_{p=N-j+2}^{N}a_{p}\ \omega_{n}^{p-(N-j+2)}\right)$
		$\displaystyle=\sqrt{\frac{1}{N}}\left(\sum_{p=1}^{N-j+1}a_{p}\ \omega_{n}^{p+j% -2}+\sum_{p=N-j+2}^{N}a_{p}\ \omega_{n}^{p+j-2}\right)$
		$\displaystyle=\sqrt{\frac{1}{N}}\sum_{p=1}^{N}a_{p}\ \omega_{n}^{p+j-2}$
		$\displaystyle=\left(\sum_{p=1}^{N}a_{p}\ \omega_{n}^{p-1}\right)\ \frac{\omega% _{n}^{j-1}}{\sqrt{N}}$
		$\displaystyle=\mu_{n}\ \frac{\omega_{n}^{j-1}}{\sqrt{N}},$

which is the $j$ -th component of $\mu_{n}\mathbf{w}_{n}.$ This shows $\mu_{n}$ is an eigenvalue of $\mathbf{A}$ with associated eigenvector $\mathbf{w}_{n}.$ Next, note that

	$\displaystyle\omega_{N-n}$	$\displaystyle=e^{\frac{2\pi i(N-n)}{N}}$
		$\displaystyle=e^{-\frac{2\pi in}{N}}=\overline{\omega_{n}},$		(2.8)

for each $1\leq n\leq N,$ which implies $\mathbf{w}_{n}^{*}=\mathbf{w}^{\text{T}}_{N-n}.$ Finally, for all $1\leq m,n\leq N$ with $m\neq n,$ we use the finite geometric series identity

\displaystyle\sum_{j=1}^{N}r^{j}

\displaystyle=\frac{r-r^{N+1}}{1-r}

for all $r\neq 1,$ to see that

	$\displaystyle\mathbf{w}_{m}^{*}\mathbf{w}_{n}$	$\displaystyle=\frac{1}{N}\sum_{p=1}^{N}\omega_{N-m}^{p-1}\omega_{n}^{p-1}$
		$\displaystyle=\frac{1}{N}\sum_{p=1}^{N}\omega_{N-m+n}^{p-1}$
		$\displaystyle=\frac{1}{N}\ \left(\frac{1-\omega_{N-m+n}^{N}}{1-\omega_{N-m+n}}\right)$
		$\displaystyle=0.$

If $n=m,$ then

	$\displaystyle\mathbf{w}_{n}^{*}\mathbf{w}_{n}$	$\displaystyle=\frac{1}{N}\sum_{p=1}^{N}\omega_{N-n}^{p-1}\omega_{n}^{p-1}$
		$\displaystyle=\frac{1}{N}\sum_{p=1}^{N}\omega_{N}^{p-1}$
		$\displaystyle=\frac{N}{N}=1.$

Since $\left\{\mathbf{w}_{n}\right\}_{n=1}^{N}$ is a set of $N$ orthonormal vectors in $\mathbb{C}^{N},$ we deduce it is a basis for $\mathbb{C}^{N}.$ This proves the first statement.

We prove the second statement. Let $k\geq 0$ be fixed. Because $\mu_{n}$ is an eigenvalue of $\mathbf{A},$ note that $\mu_{n}^{k}$ is an eigenvalue of $\mathbf{A}^{k}$ with associated eigenvector $\mathbf{w}_{n}.$ For each $1\leq m\leq N,$ the orthonormality of the collection of eigenvectors $\left\{\mathbf{w}_{n}\right\}_{n=1}^{N}$ yields

	$\displaystyle\sum_{n=1}^{N}\mu_{n}^{k}\ \mathbf{w}_{n}\mathbf{w}_{N-n}^{\text{% T}}\mathbf{w}_{m}$	$\displaystyle=\sum_{n=1}^{N}\mu_{n}^{k}\ \mathbf{w}_{n}\ \delta_{m,n}$
		$\displaystyle=\mu_{m}^{k}\mathbf{w}_{m}=\mathbf{A}^{k}\mathbf{w}_{m}.$

This shows that the linear operator $\mathbf{A}^{k}$ and the linear operator $\sum_{n=1}^{N}\lambda_{n}^{k}\ \mathbf{w}_{n}\mathbf{w}_{N-n}^{\text{T}}$ agree on the basis of eigenvectors for each $k\geq 0$ . Because every vector in $\mathbb{C}^{N}$ is a unique linear combination of vectors from this basis, the operators agree on all of $\mathbb{C}^{N}$ and thus are identical. This proves the second statement.

Now, we prove the last statement. By the definition of the matrix exponential in (2.1), we have

$\displaystyle\sum_{k\geq 0}\frac{(t\mathbf{A})^{k}}{k!}$	$\displaystyle=\sum_{k\geq 0}\frac{t^{k}\mathbf{A}^{k}}{k!}$
	$\displaystyle=\sum_{k\geq 0}\frac{t^{k}}{k!}\sum_{n=1}^{N}\mu_{n}^{k}\ \mathbf% {w}_{n}\mathbf{w}_{N-n}^{\text{T}}$
	$\displaystyle=\sum_{n=1}^{N}\sum_{k\geq 0}\frac{(\mu_{n}t)^{k}}{k!}\mathbf{w}_% {n}\mathbf{w}_{N-n}^{\text{T}}$	(2.9)
	$\displaystyle=\sum_{n=1}^{N}e^{\mu_{n}t}\ \mathbf{w}_{n}\mathbf{w}_{N-n}^{% \text{T}},$	(2.10)

in which we obtained (2.10) by simplifying the inner sum in (2.9) using the Taylor series of the ordinary exponential function:

\displaystyle e^{z}

\displaystyle=\sum_{k\geq 0}\frac{z^{k}}{k!}.

The product of circulant matrices is circulant and the scalar multiple of a circulant matrix is also circulant. Therefore, for each $t\in\mathbb{R},$ the matrix $(t\mathbf{A})^{k}$ is a circulant matrix. Since the sum of circulant matrices is circulant, the partial sum

\displaystyle\sum_{k=0}^{K}\frac{(t\mathbf{A})^{k}}{k!}

is circulant for each $K\in\mathbb{N}.$ Thus, upon taking the limit as $K\rightarrow\infty,$ we deduce $\exp(t\mathbf{A})$ is a circulant matrix. ∎

We will adopt the following notation in this thesis. If $c_{1}\ ,\ \dots\ ,\ c_{N}\in\mathbb{F},$ then we denote $\text{Circ}(c_{1}\ ,\ \dots\ ,\ c_{N})$ as the $N\times N$ circulant matrix whose first row is $(c_{1}\ ,\ \dots\ ,\ c_{N}).$

2.2 Probability Theory

Throughout this section, we let $\left(\Omega,\mathcal{A},\mathbb{P}\right)$ be a probability space, and we also fix $M,N\in\mathbb{N}.$ Recall a random vector on $\Omega$ is a measurable function between $\Omega$ and $\mathbb{R}^{N},$ in which $\mathbb{R}^{N}$ is equipped with its Borel sigma algebra.

Let $\mathbf{x}$ and $\mathbf{y}$ be two $M$ -dimensional and $N$ -dimensional random vectors on $\Omega,$ respectively. We denote $\mathbb{E}[\mathbf{x}]$ as the expectation (mean) of $\mathbf{x}$ and recall

\displaystyle\text{Cov}(\mathbf{x},\mathbf{y})

\displaystyle=\mathbb{E}\left[\mathbf{x}\mathbf{y}^{\text{T}}\right]-\mathbb{E% }\left[\mathbf{x}\right]\ \left(\mathbb{E}\left[\mathbf{y}\right]\right)^{% \text{T}}

is the $M\times N$ cross-covariance matrix between $\mathbf{x}$ and $\mathbf{y}.$ The covariance matrix of $\mathbf{x}$ is the positive semidefinite matrix $\text{Cov}(\mathbf{x},\mathbf{x}).$

Let $\mathbf{x}$ be an $M$ -dimensional random vector on $\Omega.$ Recall $\mathbf{x}$ is ( $M$ -dimensional) Gaussian if its probability density function satisfies

\displaystyle f_{\mathbf{x}}(\mathbf{z})

\displaystyle=\frac{1}{\left(2\pi\right)^{\frac{M}{2}}\ \sqrt{|\det(\mathbf{S}% )|}}\ e^{\frac{-(\mathbf{z}-\boldsymbol{\mu})^{\text{T}}\mathbf{S}^{-1}(% \mathbf{z}-\boldsymbol{\mu})}{2}}.

(2.11)

for some $\boldsymbol{\mu}\in\mathbb{R}^{M}$ and positive definite matrix $\mathbf{S}\in\text{Mat}_{M,M}(\mathbb{R}).$ Note that $\mathbf{x}$ has mean $\boldsymbol{\mu}$ and covariance matrix $\mathbf{S},$ and we write $\mathbf{x}\sim\mathcal{N}(\boldsymbol{\mu},\mathbf{S}).$ Recall that any linear combination of independent Gaussian random vectors is also a Gaussian random vector. Moreover, any linear transformation of a Gaussian random vector is also a Gaussian random vector.

2.3 Stochastic Processes

Throughout this section, we let $\left(\Omega,\mathcal{A},\mathbb{P}\right)$ be probability space, $\left\{\mathcal{F}_{t}\right\}_{t\geq 0}$ be a filtration on $\Omega,$ and we fix $N\in\mathbb{N}.$ Recall an event $A\in\mathcal{A}$ occurs ( $\mathbb{P}$ )-almost surely (with probability one) if $\mathbb{P}(A^{c})=0.$

Recall an ( $N$ -dimensional) stochastic process on $\Omega$ is a map $\mathbf{x}:[0,\infty)\times\Omega\rightarrow\mathbb{R}^{N}$ such that for each fixed $t\geq 0,$ the map $\mathbf{x}(t,\cdot):\Omega\rightarrow\mathbb{R}^{N}$ is an $N$ -dimensional random vector on $\Omega.$ A realization (sample path/trajectory) of $\mathbf{x}$ is a map of the form $\mathbf{x}(\cdot,\omega):[0,\infty)\rightarrow\mathbb{R}^{N},$ where $\omega\in\Omega$ is fixed. From henceforth, we adopt the notation $\left\{\mathbf{x}(t)\right\}_{t\geq 0},$ in which we write $\mathbf{x}(t)=\mathbf{x}(t,\cdot)$ as the random vector at the time $t.$ The process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is adapted to the filtration $\left\{\mathcal{F}_{t}\right\}_{t\geq 0}$ if $\mathbf{x}(t)$ is a random vector on $\Omega$ specifically equipped with the sigma algebra $\mathcal{F}_{t}$ for each $t\geq 0.$ In particular, the standard (natural) filtration of the process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is the filtration $\left\{\mathcal{S}_{t}\right\}_{t\geq 0},$ such that $\mathcal{S}_{t}$ is the sigma algebra generated by the collection of random vectors $\left\{\mathbf{x}(s)\right\}_{s=0}^{t}$ for each $t\geq 0.$ Note that $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is always adapted to its own standard filtration.

Let $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ be a stochastic process. We recall the following time dependent statistical functions associated with the process. The mean function, autocovariance function, and the covariance function of $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ are defined respectively as the maps $\boldsymbol{\mu}_{\mathbf{x}}:[0,\infty)\rightarrow\mathbb{R}^{N}\ ,\ % \boldsymbol{\Gamma}_{\mathbf{x}}:[0,\infty)^{2}\rightarrow\text{Mat}_{N,N}(% \mathbb{R}),$ and $\mathbf{V}_{\mathbf{x}}:[0,\infty)\rightarrow\text{Mat}_{N,N}(\mathbb{R}),$ where

$\displaystyle\boldsymbol{\mu}_{\mathbf{x}}(t)$	$\displaystyle=\mathbb{E}[\mathbf{x}(t)]$	(2.12)
$\displaystyle\boldsymbol{\Gamma}_{\mathbf{x}}(s,t)$	$\displaystyle=\text{Cov}(\mathbf{x}(s),\mathbf{x}(t))$	(2.13)
$\displaystyle\mathbf{V}_{\mathbf{x}}(t)$	$\displaystyle=\boldsymbol{\Gamma}_{\mathbf{x}}(t,t).$	(2.14)

The process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is (strictly) stationary if the cumulative distribution functions of $\mathbf{x}(t)$ and $\mathbf{x}(t+h)$ are identical for all $t,h\geq 0.$ The process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is weakly (wide sense) stationary if its mean function $\boldsymbol{\mu}_{\mathbf{x}}$ is identically constant and the map $(t,h)\mapsto\boldsymbol{\Gamma}_{\mathbf{x}}(t,t+h)$ between $[0,\infty)^{2}$ and $\text{Mat}_{N,N}(\mathbb{R})$ depends only on the variable $h.$ Note that all stationary processes are weakly stationary.

Let $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ be a stochastic process adapted to the filtration $\left\{\mathcal{F}_{t}\right\}_{t\geq 0}.$ The process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is Markov with respect to the filtration $\left\{\mathcal{F}_{t}\right\}_{t\geq 0}$ if for all $s,t\geq 0$ and any measurable, bounded function $f:\mathbb{R}^{N}\rightarrow\mathbb{R},$ we have

\displaystyle\mathbb{E}\left[f(\mathbf{x}(s+t))\ |\ \mathcal{F}_{s}\right]

\displaystyle=\mathbb{E}\left[f(\mathbf{x}(s+t))\ |\ \mathbf{x}(s)\right]

almost surely. We simply say $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is Markov if it is Markov with respect to its standard filtration. Moreover, the process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is a diffusion process if it is a Markov process having continuous sample paths with probability one. The process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is a martingale with respect to the filtration $\left\{\mathcal{F}_{t}\right\}_{t\geq 0}$ if for all $0\leq s\leq t,$ the expectation $\mathbb{E}\left[\left|\mathbf{x}(t)\right|\right]$ exists and $\mathbb{E}\left[\mathbf{x}(t)\ |\ \mathcal{F}_{s}\right]=\mathbf{x}(s).$ We simply say $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is a martingale if it is a martingale with respect to its standard filtration.

Let $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ be a stochastic process. The process is said to be ( $N$ -dimensional) Gaussian if for all $K\in\mathbb{N}$ and all times $t_{1}\ ,\ \dots\ ,\ t_{K}\geq 0$ and all $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{K}\in\mathbb{R},$ the linear combination $\sum_{k=1}^{K}\alpha_{k}\ \mathbf{x}(t_{k})$ is a Gaussian random vector. For Gaussian processes, the notions of stationarity and weak stationarity are equivalent. The process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is the ( $N$ -dimensional) standard Wiener process (Brownian Motion) if it is a diffusion process such that $\mathbf{x}(0)=\mathbf{0}$ almost surely with independent, Gaussian increments: for all $0\leq s_{1}<t_{1}\leq s_{2}<t_{2},$ both the increments $\mathbf{x}(t_{2})-\mathbf{x}(s_{2})$ and $\mathbf{x}(t_{1})-\mathbf{x}(s_{1})$ are independent and $\mathbf{x}(t_{1})-\mathbf{x}(s_{1})\sim\mathcal{N}(\mathbf{0},(t_{1}-s_{1})\ % \mathbf{I}),$ where $\mathbf{I}$ is the $N\times N$ identity matrix. Note that the Wiener process is a Gaussian process and a martingale. Furthermore, its component processes $\left\{w_{m}(t)\right\}_{t\geq 0}$ are independent for each $1\leq m\leq M.$

Let $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ be a Markov process adapted to the filtration $\left\{\mathcal{F}_{t}\right\}_{t\geq 0},$ and let the map $(t,\omega)\mapsto\mathbf{x}(t)(\omega)$ defined between $[0,\infty)\times\Omega$ and $\mathbb{R}^{N}$ be $\mathcal{B}([0,\infty))\times\mathcal{A}$ measurable, where $\mathcal{B}([0,\infty))$ is the Borel sigma algebra on $[0,\infty).$ Let $\left\{\theta_{t}\right\}_{t\geq 0}$ be a collection of shift operators on $\Omega$ such that $\mathbf{x}(t)(\theta_{s})=\mathbf{x}(s+t)$ for all $s,t\geq 0.$ For each $t\geq 0$ and $B\subset\mathbb{R}^{N}$ and $\mathbf{y}\in\mathbb{R}^{N},$ let

p_{t}(\mathbf{y},B)=\mathbb{P}(\mathbf{x}(t)\in B\ |\ \mathbf{x}(0)=\mathbf{y}).

A probability measure $\pi$ on $\mathbb{R}^{N}$ equipped with the Borel sigma algebra is invariant with respect to $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ if

\displaystyle\int_{\mathbb{R}^{N}}p_{t}(\mathbf{y},B)\ d\pi(\mathbf{y})

\displaystyle=\pi(B)

for each $t\geq 0$ and each $B\subset\mathbb{R}^{N}.$ A set $A\in\mathcal{A}$ is shift invariant if $\theta_{t}^{-1}(A)=A$ for all $t\geq 0.$ The Markov process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is ergodic if it admits an invariant probability measure $\pi$ satisfying $\mathbb{P}_{\pi}(A)\in\{0,1\}$ for every shift invariant set $A,$ in which $\mathbb{P}_{\pi}$ is the probability measure defined by

\displaystyle\mathbb{P}_{\pi}(B)

\displaystyle=\int_{\mathbb{R}^{N}}p_{t}(\mathbf{y},B)\ d\pi(\mathbf{y}).

(2.15)

The most important theorem pertaining to Markov ergodic processes is Birkhoff’s Ergodic Theorem [32].

Theorem 2.4 (Birkhoff’s Ergodic Theorem).

If $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ is a Markov ergodic process, then for any $p\in\mathbb{N}$ and for any measurable $f:\mathbb{R}^{N}\rightarrow\mathbb{R}$ satisfying the condition that $\int_{\mathbb{R}^{N}}(f(\mathbf{y}))^{p}\ d\pi(\mathbf{y})$ exists, we have

\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}f(\mathbf{x}(s))% \ ds

\displaystyle=\mathbb{E}_{\pi}[f(\mathbf{x}(0))].

(2.16)

$\mathbb{P}_{\pi}$ -almost surely, where

\displaystyle\mathbb{E}_{\pi}[f(\mathbf{y})]

\displaystyle=\int_{\mathbb{R}^{N}}f(\mathbf{y})\ d\pi(\mathbf{y}).

(2.17)

2.4 Stochastic Differential Equations

As before, let $(\Omega,\mathcal{A},\mathbb{P})$ be a probability space and fix $M,N\in\mathbb{N}.$ Let $\left\{\mathbf{w}(t)\right\}_{t\geq 0}$ be the standard $M$ -dimensional Wiener process. In this section, we briefly review the theory of stochastic differential equations.

Let $\mathbf{F}:[0,\infty)\times\Omega\rightarrow\text{Mat}_{N,M}(\mathbb{R})$ be a map such that the $MN$ -dimensional stochastic process $\left\{\text{vec}(\mathbf{F}(t,\cdot))\right\}_{t\geq 0}$ is adapted to the standard filtration of the Wiener process, where vec is the vectorization operator that maps $N\times M$ matrix to the $NM$ -dimensional vector formed by stacking the columns of the matrix one after another. The time integral of $\mathbf{F}$ and stochastic integral of $\mathbf{F}$ with respect to the Wiener process are respectively defined as

	$\displaystyle\int_{0}^{t}\mathbf{F}(s,\omega)\ ds$	$\displaystyle=\lim_{K\rightarrow\infty}\sum_{k=1}^{K-1}\mathbf{F}(s_{k},\omega% )\ \left(s_{k+1}-s_{k}\right)$		(2.18)
	$\displaystyle\int_{0}^{t}\mathbf{F}(s,\omega)\ d\mathbf{w}(s)$	$\displaystyle=\lim_{K\rightarrow\infty}\sum_{k=1}^{K-1}\mathbf{F}(s_{k},\omega% )\ \left(\mathbf{w}(s_{k+1})-\mathbf{w}(s_{k})\right),$		(2.19)

where $\left\{s_{k}\right\}_{k=1}^{K}$ is an increasing sequence in $[0,t]$ satisfying $s_{1}=0$ and $s_{K}=t.$ Note the limits in (2.18) and (2.19) are taken in the sense of $L^{2}$ i.e. the mean-squared limit.

An Ito process is a process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ adapted to the standard filtration of the Wiener process satisfying

\displaystyle\mathbf{x}(t)

\displaystyle=\mathbf{x}(0)+\int_{0}^{t}\mathbf{b}(s,\mathbf{x}(s))\ ds+\int_{% 0}^{t}\boldsymbol{\sigma}(s,\mathbf{x}(s))\ d\mathbf{w}(s),

(2.20)

where $\mathbf{b}:[0,\infty)\times\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ and $\boldsymbol{\sigma}:[0,\infty)\times\mathbb{R}^{N}\rightarrow\text{Mat}_{N,M}(% \mathbb{R})$ are measurable functions. Informally, we write (2.20) as an equation of the form

\displaystyle d\mathbf{x}(t)

\displaystyle=\mathbf{b}(t,\mathbf{x}(t))\ dt+\boldsymbol{\sigma}(t,\mathbf{x}% (t))\ d\mathbf{w}(t),

(2.21)

which we refer to as a stochastic differential equation (SDE). We interpret any Ito process in (2.20) as a solution to the SDE (2.21).

We state the existence and uniqueness criteria [21, Theorem 5.2.1] pertaining to the solution of the SDE (2.21).

Theorem 2.5 (SDE Solution Existence and Uniqueness Conditions).

Consider the SDE (2.21). Suppose there exist constants $C,D$ satisfying

	$\displaystyle\\|\mathbf{b}(t,\mathbf{y})\\|+\\|\boldsymbol{\sigma}(t,\mathbf{y})% \\|_{F}$	$\displaystyle\leq C(1+\\|\mathbf{y}\\|)$		(2.22)
	$\displaystyle\\|\mathbf{b}(t,\mathbf{z})-\mathbf{b}(t,\mathbf{y})\\|$	$\displaystyle\leq D(\\|\mathbf{z}-\mathbf{y}\\|),$		(2.23)

for all $\mathbf{y},\mathbf{z}\in\mathbb{R}^{N}$ and all $t\geq 0,$ where $\|\cdot\|_{F}$ is the Frobenius norm and $\|\cdot\|$ is the Euclidean $2$ -norm. Suppose further that the initial random vector $\mathbf{x}(0)$ satisfies $\mathbb{E}[\|\mathbf{x}(0)\|^{2}]$ being finite and $\mathbf{x}(0)$ is independent of the Wiener process. Then, the SDE (2.21) has a unique solution for $\mathbf{x}(t),$ which is a diffusion process adapted to the standard filtration of the Wiener process. Moreover, the solution is square integrable i.e. $\mathbb{E}\left[\int_{0}^{t}\|\mathbf{x}(s)\|^{2}\ ds\right]$ is finite for each $t\geq 0.$

Chapter 3 The Ornstein-Uhlenbeck Process

Let $(\Omega,\mathcal{A},\mathbb{P})$ be a probability space. Recall in the introduction, we defined an $N$ -dimensional OU Process on this probability space, namely the $N$ -dimensional stochastic process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ satisfying the SDE

\displaystyle d\mathbf{x}(t)=-\mathbf{B}\ \mathbf{x}(t)\ dt+\boldsymbol{\Sigma% }\ d\mathbf{w}(t).

(3.1)

The linear operator $\mathbf{B}\in\text{Mat}_{N,N}(\mathbb{R})$ is called the friction operator which satisfies the condition $-\mathbf{B}$ is Hurwitz. The linear operator $\boldsymbol{\Sigma}\in\text{Mat}_{M,N}(\mathbb{R})$ is called the volatility operator. The process $\left\{\mathbf{w}(t)\right\}_{t\geq 0}$ is the $M$ -dimensional standard Wiener Process. Moreover, we assume that $\mathbf{x}(0)\sim\mathcal{N}(\mathbf{0},\mathbf{S})$ is independent of the Wiener process, in which we recall $\mathbf{S}\in\text{Mat}_{N,N}(\mathbb{R})$ is the stationary covariance operator satisfying the Lyapunov equation

\displaystyle\mathbf{B}\mathbf{S}+\mathbf{S}\mathbf{B}^{\text{T}}

\displaystyle=2\mathbf{D},

(3.2)

and

\displaystyle\mathbf{D}

\displaystyle=\frac{\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\text{T}}}{2}

(3.3)

is the diffusion operator. We view the friction operator $\mathbf{B}$ and volatility operator $\boldsymbol{\Sigma}$ as OU model parameters.

In this thesis, we will utilize the matrix-valued function $\mathbf{G}:\mathbb{R}\rightarrow\text{Mat}_{N,N}(\mathbb{R})$ defined by

\displaystyle\mathbf{G}(t)

\displaystyle=\exp(-t\mathbf{B}),

(3.4)

which is commonly known as the Green’s function (propagator) in the literature [18]. We note that $\mathbf{G}$ is the fundamental solution operator to the linear ordinary differential equation $\dot{\mathbf{y}}=-\mathbf{B}\mathbf{y}.$

The stationary covariance matrix $\mathbf{S}$ has an explicit formula as a matrix-valued integral involving $\mathbf{G}(t).$ We will be using this explicit formula throughout this thesis.

Theorem 3.1 (OU Process Stationary Covariance Matrix Explicit Formula).

The stationary covariance matrix $\mathbf{S}$ is explicitly

\displaystyle\mathbf{S}

\displaystyle=2\int_{0}^{\infty}\mathbf{G}(t)\ \mathbf{D}\ \mathbf{G}^{\text{T% }}(t)\ dt.

(3.5)

Moreover, $\mathbf{S}$ is positive semidefinite.

Proof.

First, we verify that (3.5) is indeed a solution to the Lyapunov equation (3.2). Substituting the right hand side of (3.5) for $\mathbf{S}$ into the left hand side of the equation (3.2) and utilizing the property $\mathbf{G}^{\prime}=-\mathbf{B}\mathbf{G}$ from Theorem 2.1, we have

	$\displaystyle\mathbf{B}\mathbf{S}+\mathbf{S}\mathbf{B}^{\text{T}}$	$\displaystyle=2\mathbf{B}\left(\int_{0}^{\infty}\mathbf{G}(t)\ \mathbf{D}\ % \mathbf{G}^{\text{T}}(t)\ dt\right)+2\left(\int_{0}^{\infty}\mathbf{G}(t)\ % \mathbf{D}\ \mathbf{G}^{\text{T}}(t)\ dt\right)\mathbf{B}^{\text{T}}$
		$\displaystyle=2\int_{0}^{\infty}\mathbf{B}\ \mathbf{G}(t)\ \mathbf{D}\ \mathbf% {G}^{\text{T}}(t)\ dt+2\int_{0}^{\infty}\mathbf{G}(t)\ \mathbf{D}\ \left(% \mathbf{B}\ \mathbf{G}(t)\right)^{\text{T}}dt$
		$\displaystyle=-2\int_{0}^{\infty}\mathbf{G}^{\prime}(t)\ \mathbf{D}\ \mathbf{G% }^{\text{T}}(t)\ dt-2\int_{0}^{\infty}\mathbf{G}(t)\ \mathbf{D}\ \left(\mathbf% {G}^{\prime}(t)\right)^{\text{T}}dt$
		$\displaystyle=-2\int_{0}^{\infty}\left(\mathbf{G}(t)\ \mathbf{D}\ \left(% \mathbf{G}^{\text{T}}(t)\right)\right)^{\prime}\ dt$
		$\displaystyle=\lim_{t\rightarrow\infty}-2\ \mathbf{G}(t)\ \mathbf{D}\ \mathbf{% G}^{\text{T}}(t)+2\ \mathbf{G}(0)\ \mathbf{D}\ \mathbf{G}^{\text{T}}(0)$
		$\displaystyle=\mathbf{0}+2\mathbf{D}=2\mathbf{D}.$

Thus, the right hand side of (3.5) is a solution of (3.2).

Next, we verify our solution to the Lyapunov equation (3.2) is unique. To this end, we let $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ be two solutions to (3.2). Via the product rule, we obtain

	$\displaystyle\left(\mathbf{G}(t)\ (\mathbf{S}_{1}-\mathbf{S}_{2})\ \mathbf{G}^% {\text{T}}(t)\right)^{\prime}$	$\displaystyle=-\mathbf{G}(t)\ \mathbf{B}\ (\mathbf{S}_{1}-\mathbf{S}_{2})\ % \mathbf{G}^{\text{T}}(t)$
		$\displaystyle\qquad\qquad-\mathbf{G}(t)\ (\mathbf{S}_{1}-\mathbf{S}_{2})\ % \mathbf{B}^{\text{T}}\ \mathbf{G}^{\text{T}}(t)$
		$\displaystyle=-\mathbf{G}(t)\ \left(\mathbf{B}\mathbf{S}_{1}+\mathbf{S}_{1}% \mathbf{B}^{\text{T}}-\mathbf{B}\mathbf{S}_{2}-\mathbf{S}_{2}\mathbf{B}^{\text% {T}}\right)\ \mathbf{G}^{\text{T}}(t)$
		$\displaystyle=-\mathbf{G}(t)\ \left(2\mathbf{D}-2\mathbf{D}\right)\ \mathbf{G}% ^{\text{T}}(t)=\mathbf{0},$

which implies $\mathbf{G}(t)\ (\mathbf{S}_{1}-\mathbf{S}_{2})\ \mathbf{G}^{\text{T}}(t)$ must be a constant function of $t.$ Since $-\mathbf{B}$ is Hurwitz, we have

\displaystyle\lim_{t\rightarrow\infty}\mathbf{G}(t)\ (\mathbf{S}_{1}-\mathbf{S% }_{2})\ \mathbf{G}^{\text{T}}(t)

\displaystyle=\mathbf{0},

which forces

\displaystyle\mathbf{G}(t)\ (\mathbf{S}_{1}-\mathbf{S}_{2})\ \mathbf{G}^{\text% {T}}(t)

\displaystyle=\mathbf{0}

(3.6)

for all $t\geq 0.$ Since $\mathbf{G}(t)$ and $\mathbf{G}^{\text{T}}(t)$ are invertible, we left multiply both sides of the equation by $\mathbf{G}(-t)$ and right multiply both sides of the equation by $\mathbf{G}^{\text{T}}(-t)$ to obtain the equality $\mathbf{S}_{1}-\mathbf{S}_{2}=\mathbf{0}.$ This shows $\mathbf{S}_{1}=\mathbf{S}_{2},$ which means the Lyapunov equation (3.2) has only one solution.

Finally, we show $\mathbf{S}$ is positive semidefinite. Note that

	$\displaystyle\mathbf{B}\mathbf{S}^{\text{T}}+\mathbf{S}^{\text{T}}\mathbf{B}^{% \text{T}}$	$\displaystyle=\left(\mathbf{B}\mathbf{S}+\mathbf{S}\mathbf{B}^{\text{T}}\right% )^{\text{T}}$
		$\displaystyle=2\mathbf{D}^{\text{T}}=2\mathbf{D},$		(3.7)

which implies $\mathbf{S}^{\text{T}}$ solves the Lyapunov equation (3.2). But by uniqueness, we have $\mathbf{S}=\mathbf{S}^{\text{T}},$ which means $\mathbf{S}$ is symmetric. Now, for all nonzero $\mathbf{y}\in\mathbb{R}^{N},$ we have

	$\displaystyle\mathbf{y}^{\text{T}}\mathbf{S}\mathbf{y}$	$\displaystyle=2\int_{0}^{\infty}\mathbf{y}^{\text{T}}\ \mathbf{G}(t)\ \mathbf{% D}\ \mathbf{G}^{\text{T}}(t)\ \mathbf{y}\ dt$
		$\displaystyle=\int_{0}^{\infty}\mathbf{y}^{\text{T}}\ \mathbf{G}(t)\ % \boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\text{T}}\ \mathbf{G}^{\text{T}}(t)\ % \mathbf{y}\ dt$
		$\displaystyle=\int_{0}^{\infty}\left\\|\boldsymbol{\Sigma}^{\text{T}}\ \mathbf{% G}^{\text{T}}(t)\ \mathbf{y}\right\\|_{2}^{2}\ dt,$

which is non-negative because the integrand $\left\|\boldsymbol{\Sigma}^{\text{T}}\ \mathbf{G}^{\text{T}}(t)\ \mathbf{y}% \right\|_{2}^{2}$ is non-negative for all $t\geq 0.$ Hence, $\mathbf{S}$ is positive semidefinite. ∎

3.1 Statistical Properties

In this section, we will recall the important statistical properties of the OU Process. We first derive the explicit representation of the OU process i.e. the solution to the governing stochastic differential equation (3.1). The solution also reveals the essential statistics of the OU process.

Theorem 3.2 (OU Process Explicit Representation).

The unique solution of the governing OU process stochastic differential equation (3.1) is explicitly

\displaystyle\mathbf{x}(t)

\displaystyle=\mathbf{G}(t)\ \mathbf{x}(0)+\int_{0}^{t}\mathbf{G}(t-s)\ % \boldsymbol{\Sigma}\ d\mathbf{w}(s).

(3.8)

Moreover, the OU process is an Ito diffusion process that is square integrable.

Proof.

First, we show that the right hand side of (3.8) solves the stochastic differential equation (3.1). Let

\displaystyle\mathbf{y}(t)

\displaystyle=\mathbf{G}(-t)\ \mathbf{x}(t).

Using Ito’s product rule and the Ito formalisms [21, Theorem 4.2.1], namely

\displaystyle(dt)^{2}=0\ ,\ d\mathbf{w}(t)\ dt=\mathbf{0},

we obtain

$\displaystyle d\mathbf{y}(t)$	$\displaystyle=d\left(\mathbf{G}(-t)\ \mathbf{x}(t)\right)$
	$\displaystyle=d\mathbf{G}(-t)\ \mathbf{x}(t)+\mathbf{G}(-t)\ d\mathbf{x}(t)+d% \mathbf{G}(-t)\ d\mathbf{x}(t)$
	$\displaystyle=\mathbf{B}\mathbf{G}(-t)\ \mathbf{x}(t)\ dt+\mathbf{G}(-t)\left(% -\mathbf{B}\ \mathbf{x}(t)\ dt+\boldsymbol{\Sigma}\ d\mathbf{w}(t)\right)$
	$\displaystyle\qquad+\left(\mathbf{B}\mathbf{G}(-t)\ dt\right)(-\mathbf{B}\ % \mathbf{x}(t)\ dt+\boldsymbol{\Sigma}\ d\mathbf{w}(t))$
	$\displaystyle=\mathbf{G}(-t)\ \boldsymbol{\Sigma}\ d\mathbf{w}(t)+\left(% \mathbf{B}\mathbf{G}(-t)\ dt\right)(-\mathbf{B}\ \mathbf{x}(t)\ dt+\boldsymbol% {\Sigma}\ d\mathbf{w}(t))$
	$\displaystyle=\mathbf{G}(-t)\ \boldsymbol{\Sigma}\ d\mathbf{w}(t),$	(3.9)

Formally, (3.9) is the integral equation

\displaystyle\mathbf{y}(t)-\mathbf{y}(0)

\displaystyle=\int_{0}^{t}\mathbf{G}(-s)\ \boldsymbol{\Sigma}\ d\mathbf{w}(s).

(3.10)

Upon recognizing that $\mathbf{y}(0)=\mathbf{x}(0)$ and recalling that $\mathbf{G}(-t)$ is invertible with inverse $\mathbf{G}(t)$ by Theorem 2.1, we have

	$\displaystyle\mathbf{x}(t)$	$\displaystyle=\mathbf{G}(t)\ \mathbf{y}(t)$
		$\displaystyle=\mathbf{G}(t)\left(\mathbf{y}(0)+\int_{0}^{t}\mathbf{G}(-s)\ % \boldsymbol{\Sigma}\ d\mathbf{w}(s)\right)$
		$\displaystyle=\mathbf{G}(t)\left(\mathbf{x}(0)+\int_{0}^{t}\mathbf{G}(-s)\ % \boldsymbol{\Sigma}\ d\mathbf{w}(s)\right)$
		$\displaystyle=\mathbf{G}(t)\ \mathbf{x}(0)+\int_{0}^{t}\mathbf{G}(t-s)\ % \boldsymbol{\Sigma}\ d\mathbf{w}(s).$

Thus, the right side of (3.8) is a solution to the governing stochastic differential equation (3.1).

Next, we show that (3.1) has a unique solution that is an Ito diffusion process that is square integrable. We need to verify the existence and uniqueness criteria listed in Theorem 2.5 for the OU process. Define $\mathbf{b}:[0,\infty)\times\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ and $\boldsymbol{\sigma}:[0,\infty)\times\mathbb{R}^{N}\rightarrow\text{Mat}_{N,M}(% \mathbb{R})$ by

	$\displaystyle\mathbf{b}(t,\mathbf{y})$	$\displaystyle=-\mathbf{B}\mathbf{y}$		(3.11)
	$\displaystyle\boldsymbol{\sigma}(t,\mathbf{y})$	$\displaystyle=\boldsymbol{\Sigma}.$		(3.12)

Since $\mathbf{b}(t,\cdot)$ is a linear operator that is independent of the variable $t,$ it is continuous and is hence measurable. Since $\boldsymbol{\sigma}$ is constant function in both $\mathbf{y}$ and $t,$ it is automatically measurable. Furthermore, $\mathbf{b}$ is bounded; by the definition of boundedness, there exists a constant $D$ directly satisfying (2.23). Now, let $C=\max(D,\|\boldsymbol{\Sigma}\|_{F}).$ Then

	$\displaystyle\\|\mathbf{b}(\mathbf{y},t)\\|+\\|\boldsymbol{\sigma}(\mathbf{y},t)% \\|_{F}$	$\displaystyle=\\|\mathbf{b}(\mathbf{y},t)\\|+\\|\boldsymbol{\Sigma}\\|_{F}$
		$\displaystyle\leq D\\|\mathbf{y}\\|+\\|\boldsymbol{\Sigma}\\|_{F}$
		$\displaystyle\leq C\\|\mathbf{y}\\|+C=C(1+\\|\mathbf{y}\\|).$

Hence, our choice of $C$ satisfies (2.22). Moreover, recall by assumption that $\mathbf{x}(0)\sim\mathcal{N}(\mathbf{0},\mathbf{S})$ is independent of the Wiener process, and $\mathbb{E}[\|\mathbf{x}(0)\|^{2}]$ exists by virtue of $\mathbf{x}(0)$ having a finite covariance matrix. Hence, the SDE (3.1) satisfies the assumptions listed in Theorem 2.5. This completes the proof. ∎

Using the solution of the OU process, we derive its important statistical properties.

Theorem 3.3 (OU Process is Gaussian).

The OU process is a Gaussian process.

Proof.

Consider the stochastic processes $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ and $\left\{\mathbf{z}(t)\right\}_{t\geq 0},$ where

	$\displaystyle\mathbf{y}(t)$	$\displaystyle=\mathbf{G}(t)\ \mathbf{x}(0)$
	$\displaystyle\mathbf{z}(t)$	$\displaystyle=\mathbf{G}(t)\ \int_{0}^{t}\mathbf{G}(-s)\ \boldsymbol{\Sigma}\ % d\mathbf{w}(s)$

for each $t\geq 0.$ Note that every finite linear combination of random vectors from the OU process is the sum of a finite linear combination of random vectors from $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ and a finite linear combination of random vectors from $\left\{\mathbf{z}(t)\right\}_{t\geq 0}.$ Furthermore, any linear combination of random vectors from $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ is independent of any linear combination of random vectors from $\left\{\mathbf{z}(t)\right\}_{t\geq 0}$ since $\mathbf{x}(0)$ is independent of the Wiener process.

We show that the process $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ is Gaussian. Every finite linear combination of random vectors from $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ is of the form

\displaystyle\sum_{k=1}^{K}\alpha_{k}\ \mathbf{y}(t_{k})

\displaystyle=\left(\sum_{k=1}^{K}\alpha_{k}\ \mathbf{G}(t_{k})\right)\ % \mathbf{x}(0)

(3.13)

for some $K\in\mathbb{N}$ and times $t_{1}\ ,\ \dots\ ,\ t_{K}\geq 0$ and constants $\alpha_{1}\ ,\ \dots\ ,\ \alpha_{K}\in\mathbb{R}.$ Observe that the right hand side of (3.13) is a deterministic linear transformation of the Gaussian random vector $\mathbf{x}(0),$ which means it is also Gaussian. Hence, $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ is a Gaussian process.

Next, we show $\left\{\mathbf{z}(t)\right\}_{t\geq 0}$ is a Gaussian process. Consider the map $(s,\omega)\mapsto\mathbf{G}(t-s)\ \boldsymbol{\Sigma}$ defined between $[0,\infty)\times\Omega$ and $\text{Mat}_{N,M}(\mathbb{R}).$ Since this map does not depend on the variable $\omega,$ it is a deterministic matrix-valued function whose $NM$ entry functions are deterministic real-valued functions. These entry functions are automatically adapted to the standard filtration of the Wiener process and are continuous in $s,$ which means they are square integrable over any closed, finite interval. Hence, by [21, Theorem 3.2.1] , the process $\left\{\mathbf{z}(t)\right\}_{t\geq 0}$ is a Gaussian process.

Finally, every finite linear combination of random vectors from the OU process is the sum of two independent Gaussian random vectors, one of which is a linear combination of random vectors from the process $\left\{\mathbf{y}(t)\right\}_{t\geq 0},$ while the other is a linear combination of random vectors from the process $\left\{\mathbf{z}(t)\right\}_{t\geq 0}.$ Therefore, every linear combination of random vectors from the OU process is a Gaussian random vector. This proves the claim.

∎

Theorem 3.4 (OU Process Statistics).

The OU process is a Gaussian process whose mean, auto-covariance functions are

	$\displaystyle\boldsymbol{\mu}_{\mathbf{x}}(t)$	$\displaystyle=\mathbf{0}$		(3.14)
	$\displaystyle\boldsymbol{\Gamma}_{\mathbf{x}}(s,t)$	$\displaystyle=\mathbf{G}(s)\ \mathbf{S}\ \mathbf{G}^{\text{T}}(t)+2\int_{0}^{% \min(s,t)}\mathbf{G}(s-u)\ \mathbf{D}\ \mathbf{G}^{\text{T}}(t-u)\ du,$		(3.15)

respectively.

Proof.

First, we determine the mean function of the OU process. Let $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ and $\left\{\mathbf{z}(t)\right\}_{t\geq 0}$ be the processes as defined in the proof of the previous theorem. Since $\mathbf{x}(t)=\mathbf{y}(t)+\mathbf{z}(t)$ for each $t\geq 0,$ we use the linearity of the expectation operator to obtain

$\displaystyle\boldsymbol{\mu}_{\mathbf{x}}(t)$	$\displaystyle=\boldsymbol{\mu}_{\mathbf{y}}(t)+\boldsymbol{\mu}_{\mathbf{z}}(t)$
	$\displaystyle=\mathbf{G}(t)\ \mathbb{E}[\mathbf{x}(0)]+\mathbb{E}\left[\int_{0% }^{t}\mathbf{G}(t-s)\ \boldsymbol{\Sigma}\ d\mathbf{w}(s)\right]$	(3.16)
	$\displaystyle=\mathbf{0}+\mathbf{0}=\mathbf{0},$

in which we recognized the expectation of the stochastic integral term in (3.16) vanishes due to the process $\left\{\mathbf{z}(t)\right\}_{t\geq 0}$ being a martingale.

Next, we calculate the autocovariance function of the OU process. Using the definition (2.13) and the bilinearity of the covariance operator, we have

$\displaystyle\boldsymbol{\Gamma}_{\mathbf{x}}(s,t)$	$\displaystyle=\text{Cov}\left(\mathbf{x}(s),\mathbf{x}(t)\right)$
	$\displaystyle=\text{Cov}\left(\mathbf{y}(s)+\mathbf{z}(s),\mathbf{y}(t)+% \mathbf{z}(t)\right)$	(3.17)
	$\displaystyle=\text{Cov}\left(\mathbf{y}(s),\mathbf{y}(t)\right)+\text{Cov}% \left(\mathbf{z}(s),\mathbf{z}(t)\right)$	(3.18)
	$\displaystyle=\boldsymbol{\Gamma}_{\mathbf{y}}(s,t)+\boldsymbol{\Gamma}_{% \mathbf{z}}(s,t)$
	$\displaystyle=\mathbf{G}(s)\ \text{Cov}(\mathbf{x}(0),\mathbf{x}(0))\ \mathbf{% G}^{\text{T}}(t)+\boldsymbol{\Gamma}_{\mathbf{z}}(s,t)$
	$\displaystyle=\mathbf{G}(s)\ \mathbf{S}\ \mathbf{G}^{\text{T}}(t)+\boldsymbol{% \Gamma}_{\mathbf{z}}(s,t),$

in which (3.18) is the result of expanding out the covariance expression in (3.17) using bilinearity and the fact $\mathbf{y}(s)$ and $\mathbf{z}(t)$ are independent and hence uncorrelated. We now determine $\boldsymbol{\Gamma}_{\mathbf{z}}(s,t).$ For all $s\leq t,$ we use the additivity of the stochastic integral and the bilinearity of the covariance operator to obtain

$\displaystyle\boldsymbol{\Gamma}_{\mathbf{z}}(s,t)$	$\displaystyle=\text{Cov}\left(\int_{0}^{s}\mathbf{G}(s-u)\ \boldsymbol{\Sigma}% \ d\mathbf{w}(u)\ ,\ \int_{0}^{t}\mathbf{G}(t-v)\ \boldsymbol{\Sigma}\ d% \mathbf{w}(v)\right)$
	$\displaystyle=\text{Cov}\left(\int_{0}^{s}\mathbf{G}(s-u)\ \boldsymbol{\Sigma}% \ d\mathbf{w}(u)\ ,\ \int_{0}^{s}\mathbf{G}(t-v)\ \boldsymbol{\Sigma}\ d% \mathbf{w}(v)\right)$	(3.19)
	$\displaystyle\qquad+\text{Cov}\left(\int_{0}^{s}\mathbf{G}(s-u)\ \boldsymbol{% \Sigma}\ d\mathbf{w}(u),\int_{s}^{t}\mathbf{G}(t-v)\ \boldsymbol{\Sigma}\ d% \mathbf{w}(v)\right)$
	$\displaystyle=\text{Cov}\left(\int_{0}^{s}\mathbf{G}(s-u)\ \boldsymbol{\Sigma}% \ d\mathbf{w}(u)\ ,\ \int_{0}^{s}\mathbf{G}(t-v)\ \boldsymbol{\Sigma}\ d% \mathbf{w}(v)\right)$	(3.20)
	$\displaystyle=\int_{0}^{s}\mathbf{G}(s-u)\ \boldsymbol{\Sigma}\ \left(\mathbf{% G}(t-u)\ \boldsymbol{\Sigma}\right)^{\text{T}}\ du$	(3.21)
	$\displaystyle=2\int_{0}^{s}\mathbf{G}(s-u)\ \mathbf{D}\ \mathbf{G}^{\text{T}}(% t-u)\ du,$

in which the second covariance term in (3.19) vanishes due to the Wiener process having independent increments, and (3.21) follows from evaluating the covariance term in (3.20) with Ito’s Isometry [21, Corollary 3.1.7]. Similar reasoning gives

\displaystyle\boldsymbol{\Gamma}_{\mathbf{z}}(s,t)

\displaystyle=2\int_{0}^{t}\mathbf{G}(s-u)\ \mathbf{D}\ \mathbf{G}^{\text{T}}(% t-u)\ du

for $s\geq t.$ Upon combining the results corresponding to the cases $s\leq t$ and $s\geq t,$ we have

\displaystyle\boldsymbol{\Gamma}_{\mathbf{z}}(s,t)

\displaystyle=2\int_{0}^{\min(s,t)}\mathbf{G}(s-u)\ \mathbf{D}\ \mathbf{G}^{% \text{T}}(t-u)\ du

for all $s,t\geq 0.$ Hence, we obtain the OU process autocovariance function. ∎

Based on the time-dependent statistics, we can deduce the OU process is stationary.

Theorem 3.5 (OU Process Stationary).

The OU process is stationary, and its covariance function is

\displaystyle\mathbf{V}_{\mathbf{x}}(t)

\displaystyle=\mathbf{S}

(3.22)

for all $t\geq 0.$

Proof.

Since the OU Process is a Gaussian process, it suffices to show the process is weakly stationary. In particular, we show the map $(t,h)\mapsto\boldsymbol{\Gamma}_{\mathbf{x}}(t+h,t)$ defined between $[0,\infty)^{2}$ and $\text{Mat}_{N,N}(\mathbb{R})$ depends only on $h.$ Differentiating the map $\boldsymbol{\Gamma}_{\mathbf{x}}(t+h,t)$ with respect to $t,$ we obtain

$\displaystyle\frac{\partial\ \boldsymbol{\Gamma}_{\mathbf{x}}(t+h,t)}{\partial t}$	$\displaystyle=\frac{\partial}{\partial t}\left(\mathbf{G}(t+h)\ \mathbf{S}\ % \mathbf{G}^{\text{T}}(t)+2\int_{0}^{t}\mathbf{G}(t+h-u)\ \mathbf{D}\ \mathbf{G% }^{\text{T}}(t-u)\ du\right)$	(3.23)
	$\displaystyle=-\mathbf{B}\ \mathbf{G}(t+h)\ \mathbf{S}\ \mathbf{G}^{\text{T}}(% t)-\mathbf{G}(t+h)\ \mathbf{S}\ \left(\mathbf{B}\ \mathbf{G}(t)\right)^{\text{% T}}$
	$\displaystyle\qquad\qquad+2\ \frac{\partial}{\partial t}\left(\int_{0}^{t}% \mathbf{G}(u+h)\ \mathbf{D}\ \mathbf{G}^{\text{T}}(u)\ du\right)$	(3.24)
	$\displaystyle=\mathbf{G}(t+h)\ (-\mathbf{B}\mathbf{S}-\mathbf{S}\mathbf{B}^{% \text{T}})\ \mathbf{G}^{\text{T}}(t)+\mathbf{G}(t+h)\ (2\mathbf{D})\ \mathbf{G% }^{\text{T}}(t)$	(3.25)
	$\displaystyle=\mathbf{G}(t+h)\ (-\mathbf{B}\mathbf{S}+\mathbf{S}\mathbf{B}^{% \text{T}}+2\mathbf{D})\ \mathbf{G}^{\text{T}}(t)=\mathbf{0},$	(3.26)

where the second term of (3.24) follows from applying the change of variables $u\mapsto t-u$ on the integral term within the partial derivative expression of (3.23), and the first term of (3.25) follows from rewriting the first term of (3.24) using the the fact that $\mathbf{B}$ and $\mathbf{G}(t)$ commute according to Theorem 2.1. Finally, (3.26) follows from the Lyapunov equation (3.2). Since the partial derivative of $\boldsymbol{\Gamma}_{\mathbf{x}}(t+h,t)$ vanishes on the connected set $[0,\infty)^{2},$ we conclude that $\boldsymbol{\Gamma}_{\mathbf{x}}(t+h,t)$ depends only on $h.$

Next, for each $t\geq 0,$ we have

	$\displaystyle\mathbf{V}_{\mathbf{x}}(t)$	$\displaystyle=\boldsymbol{\Gamma}_{\mathbf{x}}(t,t)$
		$\displaystyle=\boldsymbol{\Gamma}_{\mathbf{x}}(0,0)$
		$\displaystyle=\mathbf{V}_{\mathbf{x}}(0)=\mathbf{S}.$

∎

The previous theorem implies $\mathbf{x}(t)\sim\mathcal{N}(\mathbf{0},\mathbf{S})$ for each $t\geq 0,$ which trivially implies the stationary distribution of the OU process is the Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S}.$

We provide a sufficient condition for the OU process to be ergodic.

Theorem 3.6 (OU Process Ergodicity).

If the diffusion matrix $\mathbf{D}$ is invertible, then the OU process is ergodic.

Proof.

In order to prove ergodicity, we need to show that the stationary distribution of the OU process, namely the Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S},$ is unique [32].

Consider a generic (time homogeneous) Ito diffusion process $\left\{\mathbf{y}(t)\right\}_{t\geq 0}$ satisfying the stochastic differential equation

\displaystyle d\mathbf{y}(t)

\displaystyle=\mathbf{b}(\mathbf{y}(t))\ dt+\boldsymbol{\sigma}(\mathbf{y}(t))% \ d\mathbf{w}(t)

(3.27)

with $\mathbf{b}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ and $\boldsymbol{\sigma}:\mathbb{R}^{N}\rightarrow\text{Mat}_{N,M}(\mathbb{R})$ being Borel measurable functions. A sufficient condition [33, Theorem 5.3.2] for this diffusion process to have a unique stationary distribution is that there is an open, bounded domain $U\subset\mathbb{R}^{N}$ with smooth boundary $\partial U$ such that

i)

For each $\mathbf{x}\in U,$ the smallest eigenvalue of the diffusion coefficient $\frac{\boldsymbol{\sigma}(\mathbf{x})\ \boldsymbol{\sigma}^{\text{T}}(\mathbf{% x})}{2}$ is bounded away from $0.$
ii)

The process is recurrent relative to $U$ : for each $\mathbf{z}\in\mathbb{R}^{N}\setminus U,$ the mean time $\tau$ at which a path issuing from $\mathbf{z}$ reaches the set $U$ is finite and $\sup_{\mathbf{z}\in K}\mathbb{E}[\tau|\mathbf{z}]$ exists for each compact set $K$ in $\mathbb{R}^{N}.$

Moreover, in [34, Theorem 3.9], it was shown that the condition ii) is equivalent to the Foster Lyapunov drift condition, which states there is a non-negative, twice continuously differentiable function $V:\mathbb{R}^{N}\rightarrow\mathbb{R}$ for which $\mathcal{L}V$ is negative on $\mathbb{R}^{N}\setminus U.$ Here, $\mathcal{L}$ is the infinitesimal generator defined by

\displaystyle(\mathcal{L}f)(\mathbf{z})

\displaystyle=\left(\nabla f\right)^{\text{T}}\mathbf{b}(\mathbf{z})+\text{tr}% \left(\nabla^{2}f\ \frac{\boldsymbol{\sigma}(\mathbf{z})\boldsymbol{\sigma}^{% \text{T}}(\mathbf{z})}{2}\right).

(3.28)

for any twice continuously-differentiable $f:\mathbb{R}^{N}\rightarrow\mathbb{R},$ where $\nabla f$ is the gradient of $f$ and $\nabla^{2}f$ is the Hessian matrix of $f.$

First, we verify statement i) holds for the OU process. Let

\displaystyle U=\left\{\mathbf{z}\in\mathbb{R}^{N}\ :\ \|\mathbf{z}\|<\sqrt{3% \text{tr}(\mathbf{P})\ \text{tr}(\mathbf{D})}\right\},

where $\mathbf{P}\in\text{Mat}_{N,N}(\mathbb{R})$ is the unique positive definite solution to the Lyapunov equation

\displaystyle\mathbf{B}^{\text{T}}\mathbf{P}+\mathbf{P}\mathbf{B}

\displaystyle=\mathbf{I},

(3.29)

where $\mathbf{I}$ is the $N\times N$ identity matrix. Since $\mathbf{P}$ and $\mathbf{D}$ are positive definite, their traces are positive. Hence, $U$ is nonempty. It is an open ball in $\mathbb{R}^{N}$ centered around $\mathbf{0}.$ Note that the diffusion coefficient of the OU process is the constant matrix $\mathbf{D}.$ Since $\mathbf{D}$ is positive definite, its eigenvalues are all positive. So automatically, its smallest eigenvalue is bounded away from $0$ on $U.$ This means the OU process satisfies statement $\text{i}).$

Next, we verify the second statement ii) by verifying the Foster Lyapunov drift condition holds. For the OU process, the infinitesimal generator is explicitly

\displaystyle(\mathcal{L}f)(\mathbf{y})

\displaystyle=-(\nabla f)^{\text{T}}\ \mathbf{B}\mathbf{y}+\text{tr}(\nabla^{2% }f\ \mathbf{D}),

(3.30)

for each twice-continuously differentiable $f:\mathbb{R}^{N}\rightarrow\mathbb{R}.$ Consider the Lyapunov function $V:\mathbb{R}^{N}\rightarrow\mathbb{R}$ defined by

\displaystyle V(\mathbf{z})=\mathbf{z}^{\text{T}}\mathbf{P}\mathbf{z}.

Then, we compute

	$\displaystyle\nabla V(\mathbf{z})$	$\displaystyle=(\mathbf{P}+\mathbf{P}^{\text{T}})\ \mathbf{z}=2\mathbf{P}% \mathbf{z}$
	$\displaystyle\nabla^{2}V(\mathbf{z})$	$\displaystyle=2\mathbf{P}.$

For all $\mathbf{z}\in\mathbb{R}^{N}\setminus U,$ by the Lyapunov equation (3.29), we have

$\displaystyle(\mathcal{L}V)(\mathbf{z})$	$\displaystyle=-\left(2\mathbf{P}\mathbf{z}\right)^{\text{T}}\mathbf{B}\mathbf{% z}+\text{tr}(2\mathbf{P}\mathbf{D})$
	$\displaystyle=-2\mathbf{z}^{\text{T}}(\mathbf{P}\mathbf{B})\ \mathbf{z}+2\text% {tr}(\mathbf{P}\mathbf{D})$
	$\displaystyle=-\mathbf{z}^{\text{T}}(\mathbf{P}\mathbf{B}+\left(\mathbf{P}% \mathbf{B}\right)^{\text{T}})\ \mathbf{z}+2\text{tr}(\mathbf{P}\mathbf{D})$
	$\displaystyle=-\mathbf{z}^{\text{T}}(\mathbf{P}\mathbf{B}+\mathbf{B}^{\text{T}% }\mathbf{P})\ \mathbf{z}+2\text{tr}(\mathbf{P}\mathbf{D})$
	$\displaystyle=-\\|\mathbf{z}\\|^{2}+2\text{tr}(\mathbf{P}\mathbf{D})$
	$\displaystyle\leq-3\text{tr}(\mathbf{P})\ \text{tr}(\mathbf{D})+2\text{tr}(% \mathbf{P})\ \text{tr}(\mathbf{D})$	(3.31)
	$\displaystyle=-\text{tr}(\mathbf{P})\ \text{tr}(\mathbf{D})<0,$

in which the second term of (3.31) follows from the Cauchy Schwarz inequality applied to the trace of the product of two positive semidefinite matrices. Hence, the OU process satisfies statement $\text{ii)},$ and we conclude the process is ergodic. ∎

3.2 The Lead Process

In this section, we consider the OU process generated by model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ and we assume the process is ergodic. We recall the auxiliary lead process that we defined in the introduction, namely the $N\times N$ matrix-valued stochastic process $\left\{\mathbf{A}(t)\right\}_{t\geq 0},$ in which

\displaystyle\mathbf{A}(t)

\displaystyle=\frac{1}{2}\int_{0}^{t}\mathbf{x}(s)\ d\mathbf{x}^{\text{T}}(s)-% d\mathbf{x}(s)\ \mathbf{x}^{\text{T}}(s).

(3.32)

for each $t\geq 0.$ More explicitly, for each $1\leq m,n\leq N,$ the $(m,n)$ -th entry of $\mathbf{A}(t)$ is an Ito stochastic integral of the form

\displaystyle A_{m,n}(t)

\displaystyle=\frac{1}{2}\int_{0}^{t}x_{m}(s)\ dx_{n}(s)-x_{n}(s)\ dx_{m}(s),

(3.33)

where $\left\{x_{n}(t)\right\}_{t\geq 0}$ is the $n$ -th component OU process. By definition, each $\mathbf{A}(t)$ is the lead matrix corresponding to a realization of the OU process over the finite time interval $[0,t].$

We prove a strong law of large numbers identity pertaining to the lead process.

Theorem 3.7 (Lead Process Strong Law of Large Numbers Identity).

The lead process satisfies the following strong law of large numbers identity:

\displaystyle\lim_{t\rightarrow\infty}\frac{\mathbf{A}(t)}{t}

\displaystyle=\frac{\mathbf{B}\mathbf{S}-\mathbf{S}\mathbf{B}^{\text{T}}}{2}

(3.34)

almost surely.

Proof.

Using the governing stochastic differential equation (3.1), we write

	$\displaystyle\frac{1}{2t}\int_{0}^{t}\mathbf{x}(s)\ d\mathbf{x}^{\text{T}}(s)$	$\displaystyle=\frac{1}{2t}\int_{0}^{t}\mathbf{x}(s)\ \left(-\mathbf{B}\ % \mathbf{x}(s)\ ds+\boldsymbol{\Sigma}\ d\mathbf{w}(s)\right)^{\text{T}}$
		$\displaystyle=-\frac{1}{2t}\int_{0}^{t}\mathbf{x}(s)\ \mathbf{x}^{\text{T}}(s)% \ ds\ \mathbf{B}^{\text{T}}+\frac{1}{2t}\int_{0}^{t}\mathbf{x}(s)\ d\mathbf{w}% ^{\text{T}}(s)\ \boldsymbol{\Sigma}^{\text{T}}.$		(3.35)

We need to determine the almost sure limits of both integral terms in (3.35) as $t\rightarrow\infty.$

First, we show that the almost sure limit of the first term of (3.35) is $-\frac{\mathbf{S}\mathbf{B}^{\text{T}}}{2}$ as $t\rightarrow\infty.$ To this end, consider the map $f:\mathbb{R}^{N}\rightarrow\text{Mat}_{N,N}(\mathbb{R})$ defined by

\displaystyle f(\mathbf{y})

\displaystyle=\mathbf{y}\mathbf{y}^{\text{T}}.

(3.36)

This map has $N^{2}$ real-valued entry functions of the form $f_{m,n}:\mathbb{R}^{N}\rightarrow\mathbb{R}$ defined by $f_{m,n}(\mathbf{y})=y_{m}y_{n}$ for $1\leq m,n\leq N.$ Each $f_{m,n}$ is continuous and hence measurable between the spaces $\mathbb{R}^{N}$ and $\mathbb{R}$ equipped with their respective Borel sigma-algebras. Hence, $f$ itself is measurable and by Theorem 2.4, we obtain

	$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{2t}\int_{0}^{t}\mathbf{x}(s)\ % \mathbf{x}^{\text{T}}(s)\ ds\ \mathbf{B}^{\text{T}}$	$\displaystyle=\lim_{t\rightarrow\infty}\frac{1}{2t}\int_{0}^{t}f(\mathbf{x}(s)% )\ ds\ \mathbf{B}^{\text{T}}$
		$\displaystyle=\frac{\mathbb{E}[f(\mathbf{x}(0))]\ \mathbf{B}^{\text{T}}}{2}$
		$\displaystyle=\frac{\mathbf{S}\mathbf{B}^{\text{T}}}{2}$

almost surely.

Next, we show that the almost sure limit of the second term of (3.37) is $\mathbf{0}.$ To this end, consider the $N\times M$ matrix-valued stochastic process $\left\{\mathbf{R}(t)\right\}_{t\geq 0}$ defined by

\displaystyle\mathbf{R}(t)

\displaystyle=\int_{0}^{t}\mathbf{x}(s)\ d\mathbf{w}^{\text{T}}(s).

for each $t\geq 0.$ We need to show that

\displaystyle\lim_{t\rightarrow\infty}\frac{\mathbf{R}(t)}{t}

\displaystyle=\mathbf{0}

almost surely. Equivalently, we need to show each entry process of $\left\{\mathbf{R}(t)\right\}_{t\geq 0},$ namely the process $\left\{R_{n,m}(t)\right\}_{t\geq 0},$ in which

\displaystyle R_{n,m}(t)

\displaystyle=\int_{0}^{t}x_{n}(s)\ dw_{m}(s)

for each $t\geq 0,$ satisfies

\displaystyle\lim_{t\rightarrow\infty}\frac{R_{n,m}(t)}{t}

\displaystyle=0

(3.37)

almost surely. Note that $\left\{R_{n,m}(t)\right\}_{t\geq 0}$ is a continuous time square integrable martingale with respect to the standard filtration of the Wiener process because the component OU process $\left\{x_{n}(t)\right\}_{t\geq 0}$ is square integrable and is adapted to the component Wiener process $\left\{w_{m}(t)\right\}_{t\geq 0}.$ By Ito’s Isometry, for each $t\geq 0,$ we have

	$\displaystyle\mathbb{E}\left[R^{2}_{n,m}(t)\right]$	$\displaystyle=\mathbb{E}\left[\int_{0}^{t}x^{2}_{n}(s)\ ds\right]$
		$\displaystyle=\int_{0}^{t}\mathbb{E}\left[x^{2}_{n}(s)\right]\ ds$
		$\displaystyle=tS_{n,n}.$

Moreover, the discrete process $\left\{R_{n,m}(K)\right\}_{K\in\mathbb{N}\cup\left\{0\right\}}$ is a discrete time square integrable martingale.

We show that the discrete time martingale $\left\{R_{n,m}(K)\right\}_{K\in\mathbb{N}\cup\left\{0\right\}}$ satisfies the strong law of large numbers identity

\displaystyle\lim_{K\rightarrow\infty}\frac{R_{n,m}(K)}{K}

\displaystyle=0

(3.38)

almost surely, in which $K$ runs over the integers. We consider the martingale difference sequence $\left\{Y_{n,m}(K)\right\}_{K\in\mathbb{N}},$ where

\displaystyle Y_{n,m}(K)

\displaystyle=R_{n,m}(K)-R_{n,m}(K-1)

for each $K\in\mathbb{N}.$ Because $\left\{R_{n,m}(K)\right\}_{K\in\mathbb{N}\cup\{0\}}$ is a martingale, we have

	$\displaystyle\mathbb{E}\left[R_{n,m}(K)\ R_{n,m}(K-1)\right]$	$\displaystyle=\mathbb{E}\left[R_{n,m}(K-1)\ \mathbb{E}\left[R_{n,m}(K)\ \|\ % \mathcal{S}_{K-1}\right]\ \right]$
		$\displaystyle=\mathbb{E}\left[R^{2}_{n,m}(K)\right],$

where $\left\{\mathcal{S}_{t}\right\}_{t\geq 0}$ is the standard filtration of the Wiener process. Therefore,

	$\displaystyle\sum_{k=1}^{K}\frac{\mathbb{E}\left[(Y_{n,m}(k))^{2}\right]}{k^{2}}$	$\displaystyle=\sum_{k=1}^{K}\frac{\mathbb{E}\left[\left(R_{n,m}(k)-R_{n,m}(k-1% )\right)^{2}\right]}{k^{2}}$
		$\displaystyle=\sum_{k=1}^{K}\frac{\mathbb{E}\left[(R_{n,m}(k))^{2}\right]-2% \mathbb{E}\left[R_{n,m}(k)R_{n,m}(k-1)\right]+\mathbb{E}\left[(R_{n,m}(k-1))^{% 2}\right]}{k^{2}}$
		$\displaystyle=\sum_{k=1}^{K}\frac{\mathbb{E}\left[(R_{n,m}(k))^{2}\right]-2% \mathbb{E}\left[(R_{n,m}(k-1))^{2}\right]+\mathbb{E}\left[(R_{n,m}(k-1))^{2}% \right]}{k^{2}}$
		$\displaystyle=\sum_{k=1}^{K}\frac{\mathbb{E}\left[(R_{n,m}(k))^{2}\right]-% \mathbb{E}\left[(R_{n,m}(k-1))^{2}\right]}{k^{2}}$
		$\displaystyle=\sum_{k=1}^{K}\frac{kS_{n,n}-(k-1)S_{n,n}}{k^{2}}$
		$\displaystyle=\sum_{k=1}^{K}\frac{S_{n,n}}{k^{2}},$

which converges as $K\rightarrow\infty$ . By [35, Exercise 4.4.10], we deduce $\frac{R_{n,m}(K)}{K}$ converges almost surely to some limit as $K\rightarrow\infty.$ Because the sequence $\left\{k\right\}_{k\in\mathbb{N}}$ is an increasing sequence diverging to infinity, we deduce by Kronecker’s Lemma [35, Theorem 2.5.9] that

\displaystyle\lim_{K\rightarrow\infty}\frac{R_{n,m}(K)}{K}

\displaystyle=0

almost surely, which proves the identity (3.38). Now, by Doob’s submartingale inequality [35, Theorem 4.4.2], we have

	$\displaystyle\sum_{K\geq 0}\mathbb{P}\left(\sup_{t\in[K,K+1]}\left\|R_{n,m}(t)-% R_{n,m}(K)\right\|>K^{\frac{2}{3}}\right)$	$\displaystyle\leq\sum_{K\geq 0}\frac{\mathbb{E}[\left\|R_{n,m}(K+1)-R_{n,m}(K)% \right\|^{2}]}{K^{\frac{4}{3}}}$
		$\displaystyle\leq\sum_{K\geq 0}\frac{S_{n,n}}{K^{\frac{4}{3}}},$

which we note converges as $K\rightarrow\infty$ . Therefore, by the Borel Cantelli Lemma, we have

\displaystyle\mathbb{P}\left(\left\{\sup_{t\in[K,K+1]}\left|R_{n,m}(t)-R_{n,m}% (K)\right|>K^{\frac{2}{3}}\right\}\text{ infinitely often}\right)

\displaystyle=0.

This means for large $t\geq 0,$ there is a unique integer $K(t)$ with $K(t)-1\leq t\leq K(t)$ such that

$\displaystyle\left\|\frac{R_{n,m}(t)}{t}\right\|$	$\displaystyle\leq\left\|\frac{R_{n,m}(t)-R_{n,m}(K(t))}{t}\right\|+\left\|\frac{R% _{n,m}(K(t))}{t}\right\|$
	$\displaystyle=\left\|\frac{R_{n,m}(t)-R_{n,m}(K(t))}{t}\right\|+\left\|\frac{R_{n% ,m}(K(t))}{K(t)}\right\|$
	$\displaystyle\leq\frac{K(t)^{\frac{2}{3}}}{K(t)}+\left\|\frac{R_{n,m}(K(t))}{K(% t)}\right\|$
	$\displaystyle=\frac{1}{K(t)^{\frac{1}{3}}}+\left\|\frac{R_{n,m}(K(t))}{K(t)}% \right\|,$	(3.39)

almost surely. Therefore, as $t\rightarrow\infty,$ we have $K(t)\rightarrow\infty,$ which means (3.39) converges almost surely to $0.$ Hence, we deduce the identity (3.37) and conclude the second integral in (3.35) almost surely converges to $\mathbf{0}\boldsymbol{\Sigma}^{\text{T}}=\mathbf{0}.$

Thus, we conclude

	$\displaystyle\lim_{t\rightarrow\infty}\frac{\mathbf{A}(t)}{t}$	$\displaystyle=-\frac{\mathbf{S}\mathbf{B}^{\text{T}}}{2}-\left(-\frac{\mathbf{% S}\mathbf{B}^{\text{T}}}{2}\right)^{\text{T}}$
		$\displaystyle=\frac{\mathbf{B}\mathbf{S}-\mathbf{S}\mathbf{B}^{\text{T}}}{2}.$

almost surely. ∎

Because of this strong law of large numbers identity, we use the matrix

\displaystyle\mathbf{Q}

\displaystyle=\frac{\mathbf{B}\mathbf{S}-\mathbf{S}\mathbf{B}^{\text{T}}}{2}

(3.40)

as our main tool for the Cyclicity Analysis of the OU process. We refer to $\mathbf{Q}$ as the lead matrix of the OU process. In the next chapter, we will investigate the structure of this matrix in depth.

Chapter 4 Cyclicity Analysis of the Cyclic Ornstein-Uhlenbeck Process

In this chapter, we consider the $N$ -dimensional cyclic OU process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ with model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ in which $\mathbf{B}$ is a circulant friction matrix. We write $\mathbf{B}=\text{Circ}(b_{1}\ ,\ \dots\ ,\ b_{N}),$ in which $b_{p}\leq 0$ for each $p>1$ and $b_{1}>0$ is large enough so that $\mathbf{B}$ is a friction matrix. We let $\mathbf{D}$ be the $N\times N$ diffusion matrix corresponding to this OU process, as defined in (3.3). Recall by Theorem 2.3, for each $1\leq n\leq N,$

\displaystyle\mu_{n}

\displaystyle=\sum_{p=1}^{N}b_{p}\ \omega_{n}^{p-1}

(4.1)

is an eigenvalue of $\mathbf{B}$ with associated eigenvector

\displaystyle\mathbf{w}_{n}

\displaystyle=\frac{1}{\sqrt{N}}\left(1\ ,\ \omega_{n}\ ,\ \dots\ ,\ \omega_{n% }^{N-1}\right),

(4.2)

in which $\omega_{n}=e^{\frac{2\pi in}{N}}.$

Using the eigenvalues of $\mathbf{B}$ , we can determine the explicit lower bound of $b_{1}$ to ensure that $\mathbf{B}$ is a valid friction matrix.

Theorem 4.1 (Cyclic OU Process Friction Matrix Stability Condition).

$\mathbf{B}$ is a valid friction matrix if and only if $b_{1}>-\sum_{p=2}^{N}b_{p}.$

Proof.

In order for $\mathbf{B}$ to be a valid friction matrix, $-\mathbf{B}$ must be Hurwitz i.e. $\Re(\mu_{n})>0$ for all $1\leq n\leq N,$ where $\Re(z)$ is the real part of the complex number $z.$ This is equivalent to the condition that $\min_{1\leq n\leq N}\Re(\mu_{n})>0.$

We determine $\min_{1\leq n\leq N}\Re(\mu_{n})$ explicitly. For each $1\leq n\leq N,$ we extract the real part of $\mu_{n}$ to obtain

$\displaystyle\Re(\mu_{n})$	$\displaystyle=\sum_{p=1}^{N}b_{p}\ \cos\left(\frac{2\pi n(p-1)}{N}\right)$
	$\displaystyle=b_{1}+\sum_{p=2}^{N}b_{p}\ \cos\left(\frac{2\pi n(p-1)}{N}\right)$
	$\displaystyle\geq b_{1}+\sum_{p=2}^{N}b_{p}$	(4.3)
	$\displaystyle=\Re(\mu_{N}),$

in which the inequality (4.3) follows from the fact that $b_{p}$ is non-positive and $b_{p}\cos\left(\frac{2\pi n(p-1)}{N}\right)\geq b_{p}$ for each $1<p\leq N.$ Hence, we see

\displaystyle\min_{1\leq n\leq N}\Re(\mu_{n})

\displaystyle=\Re(\mu_{N}),

which means $\Re(\lambda_{N})>0$ if and only if

\displaystyle b_{1}+\sum_{p=2}^{N}b_{p}>0.

Rearranging the above inequality yields the claim. ∎

In this chapter, we consider the signal propagation model that we introduced Chapter 1 governed by the aforementioned cyclic OU process. In this model, there is a one-dimensional signal propagating throughout a network consisting of $N$ sensors. Writing out $\mathbf{x}(t)=\left(x_{1}(t)\ ,\ \dots\ ,\ x_{N}(t)\right)$ for each $t\geq 0,$ we interpret $x_{n}(t)$ as the measurement made by the $n$ -th sensor of the propagating effect at time $t$ . The circulant friction matrix $\mathbf{B}$ represents the sensor network structure. More specifically, for each $p>1,$ the constant $b_{p}$ is a propagation coefficient, whose magnitude quantifies how receptive the $n$ -th sensor is to the activity within the $(n+p-1)$ -th sensor for each $1\leq n\leq N,$ in which $n+p-1$ is indexed mod $N.$ The constant $b_{1}$ as a suppression coefficient, whose magnitude quantifies the rate at which the propagating signal in our model dissipates. Finally, the diffusion matrix $\mathbf{D}$ is the covariance matrix of the noise injected into the sensors. In particular, $D_{m,n}$ is the covariance between $m$ -th component noise injected into the $m$ -th sensor and the $n$ -th component noise injected into the $n$ -th sensor for each $1\leq m,n\leq N.$

For large dimension $N\in\mathbb{N}$ and different choices of cyclic OU model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ we compute the lead matrix $\mathbf{Q}$ defined in (3.40). We investigate whether Cyclicity Analysis enables us to recover the structure of the cyclic network structure induced by $\mathbf{B}.$ More specifically, we investigate whether the structure of $\mathbf{v}_{1}(\mathbf{Q}),$ the leading eigenvector of $\mathbf{Q},$ reflects the cyclic network structure. We investigate this problem under two regimes: the first regime is when $b_{1}$ is far from $\sum_{p=2}^{N}b_{p}$ and the second regime is when $b_{1}$ is close to $\sum_{p=2}^{N}b_{p}.$

4.1 Lead Matrix Explicit Representation

In this section, we derive the explicit formula of the lead matrix $\mathbf{Q}$ corresponding to the cyclic OU Process generated by the model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma}$ defined at the beginning of the chapter.

Theorem 4.2 (Cyclic OU Process Lead Matrix Explicit Formula).

The lead matrix $\mathbf{Q}$ of the cyclic OU process is explicitly

\displaystyle\mathbf{Q}

\displaystyle=\sum_{m,n=1}^{N}\frac{\sum_{p=2}^{N}b_{p}(\omega_{m}^{p-1}-% \omega_{n}^{p-1})}{2b_{1}+\sum_{q=2}^{N}b_{q}(\omega_{m}^{q-1}+\omega_{n}^{q-1% })}\ \mathbf{w}_{m}\ \mathbf{w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}% \ \mathbf{w}_{n}^{\text{T}}.

(4.4)

Proof.

First, we explicitly compute the stationary covariance matrix $\mathbf{S}$ of the OU process in (3.5). By Theorem 2.3, the Green’s Function $\mathbf{G}(t),$ defined in (3.4), is explicitly

\displaystyle\mathbf{G}(t)

\displaystyle=\sum_{n=1}^{N}e^{-\mu_{n}t}\ \mathbf{w}_{n}\ \mathbf{w}_{N-n}^{% \text{T}}.

Substituting this expression for $\mathbf{G}(t)$ into the integrand of (3.5), we have

	$\displaystyle\mathbf{S}$	$\displaystyle=2\int_{0}^{\infty}\sum_{m,n=1}^{N}e^{-\mu_{m}t}\ \mathbf{w}_{m}% \ \mathbf{w}_{N-m}^{\text{T}}\ \mathbf{D}\ \left(e^{-\mu_{n}t}\ \mathbf{w}_{n}% \ \mathbf{w}_{N-n}^{\text{T}}\right)^{\text{T}}\ dt$
		$\displaystyle=2\sum_{m,n=1}^{N}\int_{0}^{\infty}e^{-t(\mu_{m}+\mu_{n})}\ % \mathbf{w}_{m}\ \mathbf{w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ % \mathbf{w}_{n}^{\text{T}}\ dt$
		$\displaystyle=2\sum_{m,n=1}^{N}\frac{\mathbf{w}_{m}\ \mathbf{w}_{N-m}^{\text{T% }}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}_{n}^{\text{T}}}{\mu_{m}+\mu_{n}}.$

Now, upon substituting this explicit representation of $\mathbf{S}$ into the lead matrix definition in (3.40), we obtain

	$\displaystyle\mathbf{Q}$	$\displaystyle=\mathbf{B}\left(\sum_{m,n=1}^{N}\frac{\mathbf{w}_{m}\ \mathbf{w}% _{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}_{n}^{\text{T}}}{% \mu_{m}+\mu_{n}}\right)-\left(\sum_{m,n=1}^{N}\frac{\mathbf{w}_{m}\ \mathbf{w}% _{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}_{n}^{\text{T}}}{% \mu_{m}+\mu_{n}}\right)\mathbf{B}^{\text{T}}$
		$\displaystyle=\left(\sum_{m,n=1}^{N}\frac{\mu_{m}\mathbf{w}_{m}\ \mathbf{w}_{N% -m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}_{n}^{\text{T}}}{\mu_{% m}+\mu_{n}}\right)-\left(\sum_{m,n=1}^{N}\frac{\mathbf{w}_{m}\ \mathbf{w}_{N-m% }^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ (\mathbf{B}\ \mathbf{w}_{n})^{% \text{T}}}{\mu_{m}+\mu_{n}}\right)$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{\mu_{m}-\mu_{n}}{\mu_{m}+\mu_{n}}\ \mathbf% {w}_{m}\ \mathbf{w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}% _{n}^{\text{T}}$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{\sum_{p=1}^{N}b_{p}(\omega_{m}^{p-1}-% \omega_{n}^{p-1})}{\sum_{q=1}^{N}b_{q}(\omega_{m}^{q-1}+\omega_{n}^{q-1})}\ % \mathbf{w}_{m}\ \mathbf{w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ % \mathbf{w}_{n}^{\text{T}}.$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{\sum_{p=2}^{N}b_{p}(\omega_{m}^{p-1}-% \omega_{n}^{p-1})}{2b_{1}+\sum_{q=2}^{N}b_{q}(\omega_{m}^{q-1}+\omega_{n}^{q-1% })}\ \mathbf{w}_{m}\ \mathbf{w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}% \ \mathbf{w}_{n}^{\text{T}}.$

This completes the proof. ∎

As mentioned in the beginning of the chapter, we consider two regimes for our signal propagation model: the first regime in which the propagating signal exhibits significant dissipation and the second regime in which the signal does not exhibit significant dissipation. In order to better distinguish these regimes and make it suitable for future calculations, we modify our setup. For each $\epsilon>0,$ we consider the circulant friction matrix

\displaystyle\mathbf{B}(\epsilon)=\text{Circ}\left(\epsilon-\sum_{p=2}^{N}b_{p% }\ ,\ b_{2}\ ,\ \dots\ ,\ b_{N}\right)

(4.5)

in place of the original friction matrix $\mathbf{B},$ in which the constant $\epsilon>0$ is a perturbation. In this setup, the former regime occurs when $\epsilon$ is large i.e. the friction matrix is significantly perturbed from an unstable matrix, while the latter regime occurs when $\epsilon$ is small i.e. the friction matrix is slightly perturbed from an unstable matrix. So in future sections, we will consider the cyclic OU process with model parameters $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma}$ with $\mathbf{B}(\epsilon)$ defined as in (4.5). Consequently, we let $\mathbf{Q}(\epsilon)$ be the lead matrix of this cyclic OU process and $\mathbf{v}_{1}(\mathbf{Q}(\epsilon))$ be its leading eigenvector.

In future sections, we will also extend the big-Oh notation to matrices in the following way: if $g:[0,\infty)\rightarrow\mathbb{R}$ is a function, then we write $O(g(\epsilon))$ to represent an $M\times N$ matrix-valued function of $\epsilon$ such that there exists a $K\geq 0$ for which all $MN$ entry functions of $\epsilon$ are bounded above by $Kg(\epsilon).$ We will specify in context as to whether $\epsilon$ tends to infinity or tends to $0.$

4.2 Only One Nonzero Propagation Coefficient and Injection of Noise into Only One Sensor

In this section, for each perturbation $\epsilon>0,$ we consider the signal propagation model governed by the Cyclic OU process with model parameters $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma},$ in which $\mathbf{B}(\epsilon)$ is defined as in (4.5). Here, we assume exactly one propagation coefficient is nonzero i.e. $b_{p}$ is nonzero for only one index $1<p\leq N.$ The ground truth network structure is such that the $n$ -th sensor is only linked to the $(n+p-1)$ -th sensor for each $1\leq n\leq N,$ where $n+p-1$ is indexed mod $N.$ We further assume $p-1$ and $N$ are coprime. Note that this condition ensures the signal reaches every sensor in the network. To see this, $\gcd(p-1,N)=1$ implies there exist integers $x,y$ for which $Nx+(p-1)y=1$ i.e. $(p-1)y\equiv 1\mod N$ has an integer solution $y>0.$ Hence, for each $1\leq k<N,$ we have

\displaystyle N-k(p-1)y\equiv k\mod N.

In other words, if the signal is broadcast from the $N$ -th sensor, then it reaches the $k$ -th sensor after passing through $ky$ intermediate sensors.

Moreover, we assume that $\boldsymbol{\Sigma}$ is a volatility matrix such that $\mathbf{D}$ is a diagonal diffusion matrix whose $s$ -th diagonal entry $d_{s}$ is nonzero for only one index $1\leq s\leq N.$ This corresponds to the situation where noise is only injected into only one sensor, namely the $s$ -th sensor.

We derive the lead matrix $\mathbf{Q}(\epsilon)$ corresponding to this particular cyclic OU process with the assumptions imposed on $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma}.$

Theorem 4.3 (One Nonzero Propagation Coefficient, Noise in One Sensor Lead Matrix Explicit Formula).

For each $\epsilon>0,$ with the assumptions imposed on the model parameters $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma},$ the lead matrix $\mathbf{Q}(\epsilon)$ is explicitly

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\frac{d_{s}b_{p}}{N}\sum_{m,n=1}^{N}\frac{(\omega_{m}^{p-1}-% \omega_{n}^{p-1})\ \omega_{m+n}^{1-s}}{2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega% _{n}^{p-1}-2)}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}.

(4.6)

Proof.

Substituting $b_{1}=\epsilon-b_{p}$ and $b_{q}=0$ for all $1<q\leq N$ with $q\neq p$ into the general cyclic OU process lead matrix formula (4.6), we simplify the resulting summand using matrix multiplication:

	$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \mathbf% {w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \mathbf% {w}_{N-m}^{\text{T}}\ \left(\frac{d_{s}\ \omega_{N-n}^{s-1}}{\sqrt{N}}\ % \mathbf{e}_{s}\right)\ \mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \left(% \frac{d_{s}\ \omega_{N-m}^{s-1}\ \omega_{N-n}^{s-1}}{N}\right)\ \mathbf{w}_{n}% ^{\text{T}}$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \left(% \frac{d_{s}\ \omega_{-m-n}^{s-1}}{N}\right)\ \mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\frac{d_{s}b_{p}}{N}\sum_{m,n=1}^{N}\frac{(\omega_{m}^{p-1}-% \omega_{n}^{p-1})\ \omega_{m+n}^{1-s}}{2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega% _{n}^{p-1}-2)}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}.$

This proves the claim. ∎

Now, under the two regimes, we determine whether the structure of the leading eigenvector of $\mathbf{Q}(\epsilon)$ enables us to recover the expected network structure induced by $\mathbf{B}(\epsilon)$ .

4.2.1 First Regime

We start with the first regime, in which the perturbation $\epsilon$ is large. We derive the asymptotic expansion of the lead matrix under this regime.

Theorem 4.4 (One Nonzero Propagation Coefficient, Noise in One Sensor First Regime Lead Matrix Asymptotic Expansion).

If $\mathbf{Q}(\epsilon)$ is the lead matrix in (4.6), then as $\epsilon\rightarrow\infty,$ we have

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\mathbf{A}(\epsilon)+O\left(\frac{1}{\epsilon^{2}}\right),

(4.7)

where $\mathbf{A}(\epsilon)$ is the $N\times N$ skew-symmetric matrix whose $(j,k)$ -th entry is

\displaystyle A_{j,k}(\epsilon)

\displaystyle=\begin{cases}\frac{d_{s}b_{p}}{2\epsilon}&j\equiv s+1-p\mod N\ ,% \ k=s\\ -\frac{d_{s}b_{p}}{2\epsilon}&j=s\ ,\ k\equiv s+1-p\mod N\\ 0&\text{else}\end{cases}.

(4.8)

Proof.

Upon dividing the numerator and denominator of the summand in (4.6) by $2\epsilon,$ we obtain

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\frac{d_{s}b_{p}}{N}\sum_{m,n=1}^{N}\frac{\frac{(\omega_{m}^{p-1% }-\omega_{n}^{p-1})\ \omega_{m+n}^{1-s}}{2\epsilon}}{1+\frac{b_{p}(\omega_{m}^% {p-1}+\omega_{n}^{p-1}-2)}{2\epsilon}}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}.

(4.9)

Suppose $\epsilon$ satisfies

\displaystyle\left|\frac{b_{p}\ \left(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2% \right)}{2\epsilon}\right|

\displaystyle<1

for all $1\leq m,n\leq N.$ We can then rewrite the summand on the right hand side of (4.9) as a geometric series:

$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\frac{d_{s}b_{p}}{N}\sum_{m,n=1}^{N}\frac{(\omega_{m}^{p-1}-% \omega_{n}^{p-1})\ \omega_{m+n}^{1-s}}{2\epsilon}\ \sum_{r\geq 0}\left(-\frac{% b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}{2\epsilon}\right)^{r}\ \mathbf{w}_% {m}\mathbf{w}_{n}^{\text{T}}$
	$\displaystyle=\frac{d_{s}b_{p}}{N}\sum_{r\geq 0}\sum_{m,n=1}^{N}\frac{\left(-b% _{p}\left(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2\right)\right)^{r}(\omega_{m}^{p-% 1}-\omega_{n}^{p-1})\ \omega_{m+n}^{1-s}}{(2\epsilon)^{r+1}}\ \mathbf{w}_{m}% \mathbf{w}_{n}^{\text{T}}$	(4.10)
	$\displaystyle=\frac{d_{s}b_{p}}{2\epsilon N}\sum_{m,n=1}^{N}(\omega_{m}^{p-1}-% \omega_{n}^{p-1})\ \omega_{m+n}^{1-s}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}% +O\left(\frac{1}{\epsilon^{2}}\right),$	(4.11)

in which (4.11) is the result of extracting the inner double sum corresponding to the $r=0$ term in (4.10) and replacing the remaining sum over all $r>0$ with $O\left(\frac{1}{\epsilon^{2}}\right).$ Extracting the $(j,k)$ -th entry of (4.11), we get

	$\displaystyle Q_{j,k}(\epsilon)$	$\displaystyle=\frac{d_{s}b_{p}}{2\epsilon N^{2}}\sum_{m,n=1}^{N}(\omega_{m}^{p% -1}-\omega_{n}^{p-1})\ \omega_{m+n}^{1-s}\ \omega_{m}^{j-1}\ \omega_{n}^{k-1}+% O\left(\frac{1}{\epsilon^{2}}\right)$
		$\displaystyle=\frac{d_{s}b_{p}}{2\epsilon N^{2}}\sum_{m,n=1}^{N}\left(\omega_{% m}^{p+j-s-1}\omega_{n}^{k-s}-\omega_{m}^{j-s}\omega_{n}^{p+k-s-1}\right)+O% \left(\frac{1}{\epsilon^{2}}\right).$		(4.12)

To complete the proof, we evaluate the double sum appearing on the right hand side of (4.12). In order to do so, we utilize the following geometric series identity for all $q\in\mathbb{Z}$ :

\displaystyle\frac{1}{N}\sum_{n=1}^{N}\omega_{n}^{q}

\displaystyle=\begin{cases}1&q\equiv 0\mod N\\ 0&\text{else}\end{cases}.

(4.13)

Using (4.13), we obtain

	$\displaystyle Q_{j,k}(\epsilon)$	$\displaystyle=\begin{cases}\frac{d_{s}b_{p}}{2\epsilon}&j\equiv s+1-p\mod N\ ,% \ k\equiv s\mod N\\ 0&\text{else}\end{cases}$
		$\displaystyle\qquad\qquad-\begin{cases}\frac{d_{s}b_{p}}{2\epsilon}&k\equiv s+% 1-p\mod N\ ,\ j\equiv s\mod N\\ 0&\text{else}\end{cases}\ +O\left(\frac{1}{\epsilon^{2}}\right)$		(4.14)

Now, for $1\leq j,k\leq N,$ consider the systems of congruence equalities appearing in (4.14), namely

	$\displaystyle j$	$\displaystyle\equiv s+1-p\mod N\ ,\ k\equiv s\mod N$		(4.15)
	$\displaystyle k$	$\displaystyle\equiv s+1-p\mod N\ ,\ j\equiv s\mod N.$		(4.16)

Observe the second congruence equation in (4.15) holds if and only if $k-s$ divides $N.$ But since $k-s<N,$ the second congruence equation in (4.15) holds if and only if $k=s.$ Similarly, the second congruence equation in (4.16) holds if and only if $j=s.$ Moreover, note that (4.15) and (4.16) cannot simultaneously hold if $p>2$ . Otherwise, we would have both $j\equiv s\mod N$ and $j\equiv s+1-p\mod N,$ which would imply $s\equiv s+1-p\mod N$ i.e. $p\equiv 1\mod N.$ This would mean $p-1$ divides $N,$ and if $p>2,$ this would contradict the assumption $\gcd(p-1,N)=1.$ Hence, we can simplify (4.14) in the following way:

\displaystyle Q_{j,k}(\epsilon)

\displaystyle=\left(\begin{cases}\frac{d_{s}b_{p}}{2\epsilon}&j\equiv s+1-p% \mod N\ ,\ k=s\\ -\frac{d_{s}b_{p}}{2\epsilon}&k\equiv s+1-p\mod N\ ,\ j=s\\ 0&\text{else}\end{cases}\right)+\ O\left(\frac{1}{\epsilon^{2}}\right),

which proves the claim. ∎

Knowing the asymptotic structure of the lead matrix under the first regime, we now investigate the asymptotic structure of its leading eigenvector.

Theorem 4.5 (One Nonzero Propagation Coefficient, Noise in One Sensor First Regime Leading Eigenvector Asymptotic Expansion).

The leading eigenvector of $\mathbf{Q}(\epsilon)$ satisfies

\displaystyle\lim_{\epsilon\rightarrow\infty}\mathbf{v}_{1}(\mathbf{Q}(% \epsilon))

\displaystyle=\frac{i\ \mathbf{e}_{s}+\mathbf{e}_{s+1-p}}{\sqrt{2}},

(4.17)

where $\mathbf{e}_{s}$ is the $s$ -th canonical basis vector in $\mathbb{C}^{N}$ and $s+1-p$ is indexed mod $N$ as before.

Proof.

Consider $N\times N$ matrix $\mathbf{A}(\epsilon)$ whose $(j,k)$ -th entry is defined in (4.8) for all $1\leq j,k\leq N.$ Note that $\mathbf{A}(\epsilon)$ is a skew-symmetric matrix that has only two nonzero rows, namely $\frac{d_{s}b_{p}}{2\epsilon}\ \mathbf{e}_{s}$ and $-\frac{d_{s}b_{p}}{2\epsilon}\ \mathbf{e}_{s+1-p}.$ Since these rows are orthogonal, we deduce $\mathbf{A}(\epsilon)$ has rank $2.$ Moreover, the previous theorem states as $\epsilon\rightarrow\infty,$ we have

\displaystyle\|\mathbf{Q}(\epsilon)-\mathbf{A}(\epsilon)\|_{F}

\displaystyle=O\left(\frac{1}{\epsilon^{2}}\right),

which implies $\mathbf{A}(\epsilon)$ a low rank approximation of $\mathbf{Q}(\epsilon)$ as $\epsilon\rightarrow\infty.$

We now determine the leading eigenvector of $\mathbf{A}(\epsilon).$ Define

	$\displaystyle\mathbf{a}(\epsilon)$	$\displaystyle=\sqrt{\frac{d_{s}b_{p}}{2\epsilon}}\ \mathbf{e}_{s+1-p}$
	$\displaystyle\mathbf{b}(\epsilon)$	$\displaystyle=\sqrt{\frac{d_{s}b_{p}}{2\epsilon}}\ \mathbf{e}_{s}.$

Note that

\displaystyle\mathbf{A}(\epsilon)

\displaystyle=\mathbf{a}(\epsilon)\ \mathbf{b}^{\text{T}}(\epsilon)-\mathbf{b}% (\epsilon)\ \mathbf{a}^{\text{T}}(\epsilon).

(4.18)

Because $\mathbf{a}(\epsilon)$ and $\mathbf{b}(\epsilon)$ are orthogonal, they are linearly independent, and the angle between them is $\frac{\pi}{2}.$ Hence, by the formula for the eigenvalue and eigenvector of a rank $2$ skew-symmetric matrix stated in [3], the largest eigenvalue of $\mathbf{A}(\epsilon)$ is

	$\displaystyle\lambda_{1}(\mathbf{A}(\epsilon))$	$\displaystyle=i\sin\left(\frac{\pi}{2}\right)\\|\mathbf{a}(\epsilon)\\|\\|\mathbf% {b}(\epsilon)\\|$		(4.19)
		$\displaystyle=\frac{d_{s}b_{p}}{2\epsilon}\ i.$

Its associated eigenvector is the leading eigenvector of $\mathbf{A}(\epsilon),$ which is explicitly

	$\displaystyle\mathbf{v}_{1}(\mathbf{A}(\epsilon))$	$\displaystyle=\frac{-e^{-\frac{\pi}{2}i}\frac{\mathbf{a}(\epsilon)}{\\|\mathbf{% a}(\epsilon)\\|}+\frac{\mathbf{b}(\epsilon)}{\\|\mathbf{b}(\epsilon)\\|}}{\left\\|% -e^{-\frac{\pi}{2}i}\frac{\mathbf{a}(\epsilon)}{\\|\mathbf{a}(\epsilon)\\|}+% \frac{\mathbf{b}(\epsilon)}{\\|\mathbf{b}(\epsilon)\\|}\right\\|}$
		$\displaystyle=\frac{1}{\left\\|\frac{i\ \mathbf{e}_{s}-\mathbf{e}_{s+1-p}}{% \sqrt{\frac{d_{s}b_{2}}{2\epsilon}}}\right\\|}\left(\frac{i\ \mathbf{e}_{s}+% \mathbf{e}_{s+1-p}}{\sqrt{\frac{d_{s}b_{2}}{2\epsilon}}}\right)$
		$\displaystyle=\frac{1}{\sqrt{\frac{2}{\frac{d_{s}b_{p}}{2\epsilon}}}}\left(% \frac{i\ \mathbf{e}_{s}+\mathbf{e}_{s+1-p}}{\sqrt{\frac{d_{s}b_{p}}{2\epsilon}% }}\right)$
		$\displaystyle=\frac{i\ \mathbf{e}_{s}+\mathbf{e}_{s+1-p}}{\sqrt{2}}.$

∎

We interpret the result of this previous theorem. In the first regime, the leading eigenvector of the lead matrix $\mathbf{Q}(\epsilon)$ converges to a vector with exactly two nonzero components with equal moduli, namely the $s$ -th and $(s+1-p)$ -th components. The $s$ -th component lies on the real axis, while the $(s+1-p)$ -th component lies on the imaginary axis. All other components are equal to $0.$ Therefore, the structure of the asymptotic leading eigenvector under the first regime suggests that only the $s$ -th and $(s+1-p)$ -th sensors receive the dissipating signal that is propagating throughout the network. Thus, Cyclicity Analysis does not enable us to fully recover the expected network structure under the first regime.

4.2.2 Second Regime

We now investigate the second regime, in which the perturbation $\epsilon$ is small. We derive the asymptotic expansion of the lead matrix under this regime.

Theorem 4.6 (One Nonzero Propagation Coefficient and Noise in One Sensor Second Regime Lead Matrix Asymptotic Formula).

Consider the lead matrix $\mathbf{Q}(\epsilon)$ in (4.6). Then, as $\epsilon\rightarrow 0,$ the matrix $\mathbf{Q}(\epsilon)$ has an asymptotic expansion of the form

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\frac{d_{s}}{N}\sum_{m,n=1}^{N}\left(\frac{\left(\omega_{m}^{p-1% }-\omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-s}}{\omega_{m}^{p-1}+\omega_{n}^{p-% 1}-2}\right)\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}+O(\epsilon),

(4.20)

in which the summand inside the double sum is defined to be $0$ if $m=n=N.$

Proof.

We use the same technique that was used to prove Theorem 4.4. We divide the numerator and denominator of the summand in (4.6) by the quantity $b_{p}\left(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2\right)$ to obtain

	$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\frac{d_{s}b_{p}}{N}\sum_{m,n=1}^{N}\frac{\frac{\left(\omega_{m}% ^{p-1}-\omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-s}}{b_{p}\left(\omega_{m}^{p-1% }+\omega_{n}^{p-1}-2\right)}}{1+\frac{2\epsilon}{b_{p}\left(\omega_{m}^{p-1}+% \omega_{n}^{p-1}-2\right)}}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\frac{d_{s}}{N}\sum_{m,n=1}^{N}\frac{\frac{\left(\omega_{m}^{p-1% }-\omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-s}}{\omega_{m}^{p-1}+\omega_{n}^{p-% 1}-2}}{1+\frac{2\epsilon}{b_{p}\left(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2\right% )}}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}.$		(4.21)

Suppose

\displaystyle\left|\frac{2\epsilon}{b_{p}\left(\omega_{m}^{p-1}+\omega_{n}^{p-% 1}-2\right)}\right|<1,

for all $1\leq m,n\leq N.$ Then, we can rewrite (4.21) as a geometric series. We have

$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\frac{d_{s}}{N}\sum_{m,n=1}^{N}\frac{\left(\omega_{m}^{p-1}-% \omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-s}}{\omega_{m}^{p-1}+\omega_{n}^{p-1}% -2}\ \sum_{r\geq 0}\left(-\frac{2\epsilon}{b_{p}\ (\omega_{m}^{p-1}+\omega_{n}% ^{p-1}-2)}\right)^{r}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}$
	$\displaystyle=\frac{d_{s}}{N}\sum_{r\geq 0}\sum_{m,n=1}^{N}\left(\frac{\left(% \omega_{m}^{p-1}-\omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-s}}{\omega_{m}^{p-1}% +\omega_{n}^{p-1}-2}\right)\ \left(-\frac{2\epsilon}{b_{p}(\omega_{m}^{p-1}+% \omega_{n}^{p-1}-2)}\right)^{r}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}$	(4.22)
	$\displaystyle=\frac{d_{s}}{N}\sum_{m,n=1}^{N}\left(\frac{\left(\omega_{m}^{p-1% }-\omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-s}}{\omega_{m}^{p-1}+\omega_{n}^{p-% 1}-2}\right)\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}+O(\epsilon),$	(4.23)

in which we obtained (4.23) by extracting the term $r=0$ in the outer sum of (4.22) and replacing the sum over all $r>0$ with $O(\epsilon).$

Now, we need to justify defining the summand in the double sum of (4.20) to be $0$ when $m=n=N.$ Suppose $1\leq m,n\leq N$ are two indices satisfying

\displaystyle\omega_{m}^{p-1}+\omega_{n}^{p-1}=2.

(4.24)

Upon taking the moduli of both sides and using the triangle inequality, we obtain

	$\displaystyle 2$	$\displaystyle=\|\omega_{m}^{p-1}+\omega_{n}^{p-1}\|$
		$\displaystyle\leq\|\omega_{m}^{p-1}\|+\|\omega_{n}^{p-1}\|$
		$\displaystyle=1+1=2.$

Therefore, we have

\displaystyle|\omega_{m}^{p-1}+\omega_{n}^{p-1}|=|\omega_{m}^{p-1}|+|\omega_{n% }^{p-1}|.

(4.25)

Recall if $z,w\in\mathbb{C}$ satisfy $|z+w|=|z|+|w|$ i.e. $z,w$ are such that equality is obtained in the triangle inequality for complex numbers, then either one of $z,w$ must be zero or $\text{Arg}(z)=\text{Arg}(w)$ (see [36]). Since $\omega_{m}^{p-1},\omega_{n}^{p-1}$ are both nonzero, (4.25) implies $\text{Arg}(\omega_{m}^{p-1})=\text{Arg}(\omega_{n}^{p-1})$ i.e. $\frac{2\pi m(p-1)}{N}=\frac{2\pi n(p-1)}{N},$ and so $m=n.$ Now, substituting $m=n$ into (4.24), we obtain $\omega_{m}^{p-1}=1,$ which implies

m(p-1)\equiv 0\mod N.

Because $\gcd(p-1,N)=1$ divides $0,$ we see the congruence equation has a solution for $m.$ In particular, by [37, Section 2.6], we have

\displaystyle m

\displaystyle=N(k+1)

(4.26)

for $k\in\mathbb{Z}.$ But since $1\leq m\leq N,$ this forces $m=N.$ This shows (4.21) implies $m=n=N.$ Finally, note that

\displaystyle\frac{\omega_{m}^{p-1}-\omega_{m}^{p-1}}{\omega_{m}^{p-1}+\omega_% {m}^{p-1}-2}

\displaystyle=0

if $m<N.$ But taking the limit of the left hand side as $m\rightarrow N,$ we obtain via L’Hospital’s rule:

\displaystyle\lim_{m\rightarrow N}\frac{\omega_{m}^{p-1}-\omega_{m}^{p-1}}{% \omega_{m}^{p-1}+\omega_{m}^{p-1}-2}

\displaystyle=\lim_{m\rightarrow N}\frac{\frac{2\pi i(p-1)}{N}\left(\omega_{m}% ^{p-1}-\omega_{m}^{p-1}\right)}{\frac{2\pi i(p-1)}{N}\left(\omega_{m}^{p-1}+% \omega_{m}^{p-1}\right)}=0.

(4.27)

This completes the proof. ∎

4.3 Second Regime Special Case: First Propagation Coefficient Nonzero and Noise Injected into Last Sensor

At the current time of writing, we are not aware of any closed formulas for the eigenvalues and eigenvectors of the matrix defined on the right hand side of (4.20). Nevertheless, in this section, we investigate the second regime in the situation where $p=2$ and $s=N.$ That is, the ground truth network structure is such that $n$ -th sensor is linked only to the $(n+1)$ -th sensor, and noise is injected only into the $N$ -th sensor. To simplify the analysis, we will also assume that $b_{2}=-1$ and $d_{N}=1.$

We investigate the asymptotic structure of the lead matrix in the second regime.

Theorem 4.7 (First Nonzero Propagation Coefficient, Noise in Last Sensor Second Regime Lead Matrix).

Consider the lead matrix $\mathbf{Q}(\epsilon)$ in (4.6). Then, as $\epsilon\rightarrow 0,$ we have an alternate asymptotic expansion

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\mathbf{A}+O(\epsilon),

(4.28)

where $\mathbf{A}$ is the $N\times N$ skew-symmetric matrix whose $(j,k)$ -th entry is

\displaystyle A_{j,k}

\displaystyle=\begin{cases}\left(\frac{1}{2}\right)^{2N-j-k}\frac{k-j}{2N-j-k}% \binom{2N-j-k}{N-k}&j,k\neq N\\ 0&\text{else}\end{cases},

(4.29)

and

\displaystyle\binom{m}{n}

\displaystyle=\begin{cases}\frac{m!}{(m-n)!\ n!}&0\leq n\leq m\\ 0&\text{else}\end{cases}

is the binomial coefficient.

Proof.

Substituting $b_{2}=-1$ and $d_{N}=1$ into (4.20), we get

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\frac{1}{N}\sum_{m,n=1}^{N}\left(\frac{(\omega_{m}-\omega_{n})\ % \omega_{m+n}}{\omega_{m}+\omega_{n}-2}\right)\ \mathbf{w}_{m}\mathbf{w}_{n}^{% \text{T}}+O(\epsilon).

(4.30)

Extracting the $(j,k)$ -th entry of (4.30), we obtain

$\displaystyle Q_{j,k}(\epsilon)$	$\displaystyle=\frac{1}{N^{2}}\sum_{m,n=1}^{N}\frac{(\omega_{m}-\omega_{n})\ % \omega_{m+n}}{\omega_{m}+\omega_{n}-2}\ \omega_{m}^{j-1}\omega_{n}^{k-1}+O(\epsilon)$
	$\displaystyle=-\frac{1}{2N^{2}}\sum_{m,n=1}^{N}\frac{\omega_{m}^{j+1}\omega_{n% }^{k}-\omega_{m}^{j}\omega_{n}^{k+1}}{1-\frac{\omega_{m}+\omega_{n}}{2}}+O(\epsilon)$	(4.31)
	$\displaystyle=-\frac{1}{2N^{2}}\sum_{r\geq 0}\sum_{m,n=1}^{N}\left(\omega_{m}^% {j+1}\omega_{n}^{k}-\omega_{m}^{j}\omega_{n}^{k+1}\right)\left(\frac{\omega_{m% }+\omega_{n}}{2}\right)^{r}+O(\epsilon).$	(4.32)

Upon expanding the expression $\left(\frac{\omega_{m}+\omega_{n}}{2}\right)^{r}$ within the summand of (4.32) with the binomial theorem, we get

$\displaystyle Q_{j,k}(\epsilon)$	$\displaystyle=-\frac{1}{2N^{2}}\sum_{q,r\geq 0}\sum_{m,n=1}^{N}\left(\omega_{m% }^{j+1}\omega_{n}^{k}-\omega_{m}^{j}\omega_{n}^{k+1}\right)\ \left(\frac{1}{2}% \right)^{r}\ \binom{r}{q}\ \omega_{m}^{r-q}\ \omega_{n}^{q}+O(\epsilon)$
	$\displaystyle=-\frac{1}{2N^{2}}\sum_{q,r\geq 0}\sum_{m,n=1}^{N}\ \left(\frac{1% }{2}\right)^{r}\ \binom{r}{q}\left(\omega_{m}^{r+j-q+1}\ \omega_{n}^{k+q}-% \omega_{m}^{r+j-q}\ \omega_{n}^{k+q+1}\right)$	(4.33)
	$\displaystyle\qquad\qquad\qquad+O(\epsilon)$
	$\displaystyle=-\frac{1}{2}\sum_{q,r\geq 0}\begin{cases}\left(\frac{1}{2}\right% )^{r}\binom{r}{q}&r\equiv q-j-1\mod N\ ,\ q\equiv-k\mod N\\ -\left(\frac{1}{2}\right)^{r}\binom{r}{q}&r\equiv q-j\mod N\ ,\ q\equiv-k-1% \mod N\\ 0&\text{else}\end{cases}$	(4.34)
	$\displaystyle\qquad\qquad\qquad\qquad+O(\epsilon),$

where (4.34) follows from evaluating the inner double sum within (4.33) using the geometric series identity (4.13).

Now, consider the system of congruence equations

	$\displaystyle r$	$\displaystyle\equiv q-j-1\mod N\ ,\ q\equiv-k\mod N$		(4.35)
	$\displaystyle r$	$\displaystyle\equiv q-j\mod N\ ,\ q\equiv-k-1\mod N,$		(4.36)

where $r,q\geq 0$ and $1\leq j,k\leq N.$ The second congruence equation in (4.35) implies $q=vN-k$ for some integer $v\in\mathbb{N},$ while the first implies $r=uN-j-k-1$ for some integer $u\geq 2.$ Similarly, The congruence equations in (4.36) imply $q=vN-k-1$ and $r=uN-j-k-1$ for some integers $v\in\mathbb{N}$ and $u\geq 2.$ Hence, we can rewrite (4.34) as a double sum over all $u\geq 2$ and $v\geq 1$ in the following way:

$\displaystyle Q_{j,k}(\epsilon)$	$\displaystyle=-\sum_{u\geq 2,v\geq 1}\left(\frac{1}{2}\right)^{uN-j-k}\left(% \binom{uN-j-k-1}{vN-k-1}-\binom{uN-j-k-1}{vN-k}\right)$	(4.37)
	$\displaystyle\qquad\qquad\qquad+O(\epsilon)$
	$\displaystyle=\left(\frac{1}{2}\right)^{2N-j-k}\left(\binom{2N-j-k-1}{N-k}-% \binom{2N-j-k-1}{N-k-1}\right)+O(\epsilon)$	(4.38)
	$\displaystyle=\left(\frac{1}{2}\right)^{2N-j-k}\frac{2N-j-k-2(N-k)}{2N-j-k}% \binom{2N-j-k}{N-k}+O(\epsilon)$	(4.39)
	$\displaystyle=\left(\frac{1}{2}\right)^{2N-j-k}\frac{k-j}{2N-j-k}\binom{2N-j-k% }{N-k}+O(\epsilon),$

in which (4.38) follows from extracting the term corresponding $u=2,v=1$ in (4.37) and replacing the remaining sum over all other $u,v$ with $O(\epsilon)$ and (4.39) follows from rewriting (4.38) using the binomial coefficient identity [38]:

\displaystyle\binom{n-1}{k}-\binom{n-1}{k-1}

\displaystyle=\frac{n-2k}{n}\binom{n}{k}.

This completes the proof. ∎

4.3.1 Structure of Eigenvalues

We investigate the spectrum of $\mathbf{A},$ whose $(j,k)$ -th entry is defined in (4.29), for large $N$ . To make things suitable for future calculations, we will consider the Hermitian matrix $i\mathbf{A}$ in place of the original $\mathbf{A}.$ Note that the eigenvalues of $i\mathbf{A}$ are the eigenvalues of $\mathbf{A}$ multiplied by $i.$ However, the eigenvectors of $i\mathbf{A}$ are the same as those of $\mathbf{A}.$ Furthermore, we parameterize $i\mathbf{A}$ by its size $N.$ In the end, we define $\mathbf{A}_{N}$ to be the $N\times N$ Hermitian matrix whose $(j,k)$ -th entry is the quantity (4.29) multiplied by $i$ for all $1\leq j,k\leq N.$ Throughout, we let $\lambda_{n}(\mathbf{A}_{N})$ be the $n$ -th largest eigenvalue of $\mathbf{A}_{N}$ and $\left\{\mathbf{v}_{n}(\mathbf{A}_{N})\right\}_{n=1}^{N}$ be the orthonormal basis of eigenvectors such that $\mathbf{v}_{n}(\mathbf{A}_{N})$ is the associated unit eigenvector of $\mathbf{A}_{N}$ with $\lambda_{n}(\mathbf{A}_{N}).$ Finally, we denote $A_{N,j,k}$ as the $(j,k)$ -th entry of $\mathbf{A}_{N}.$

Recall the $j$ -th Gershgorin disk of $\mathbf{A}_{N}$ is the disk in the complex plane with center $0$ and radius equal to

\displaystyle R_{N,j}

\displaystyle=\sum_{k=1}^{N}|A_{N,j,k}|.

(4.40)

We concentrate on the largest eigenvalue of $\mathbf{A}_{N}$ for large $N.$ By Gershgorin’s Circle Theorem [39], all eigenvalues of $\mathbf{A}_{N}$ lie in the union of the Gershgorin disks. But we conjecture an upper bound on the radii of these disks.

Conjecture 4.1 (Gershgorin Radii Upper Bound).

The Gershgorin radii of $\mathbf{A}_{N}$ satisfy

\displaystyle\lim_{N\rightarrow\infty}R_{N,N-j}

\displaystyle=\begin{cases}\frac{\binom{2j}{j}}{2^{2j-1}}&j>0\\ 1&j=0\end{cases}.

(4.41)

To see this conjecture holds for $j=1,$ we evaluate

$\displaystyle R_{N,N-1}$	$\displaystyle=\sum_{k=1}^{N}\|A_{N,N-1,k}\|$
	$\displaystyle=\sum_{k=1}^{N}\left(\frac{1}{2}\right)^{N-k+1}\|N-k-1\|$	(4.42)
	$\displaystyle=\sum_{k=1}^{N-1}\left(\frac{1}{2}\right)^{N-k+1}(N-k-1)+\frac{1}% {2}$	(4.43)
	$\displaystyle=\frac{1}{2}-\frac{N}{2^{N}}+\frac{1}{2}$	(4.44)
	$\displaystyle=1-\frac{N}{2^{N}}$	(4.45)

in which we obtained (4.43) by rewriting (4.42) after extracting the term $k=N$ and used the identity

\displaystyle\sum_{k=1}^{N-1}kr^{k}

\displaystyle=\frac{Nr^{N+1}-r^{N+1}-Nr^{N}+r}{(r-1)^{2}}

to obtain (4.44) from (4.43). Taking the limit of (4.45) as $N\rightarrow\infty,$ we see the conjecture holds true for the case $j=1.$

Assuming Conjecture 4.1 holds, we can show the sequence of largest eigenvalues $\left\{\lambda_{1}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}$ converges.

Theorem 4.8 (Asymptotic Binomial Matrix Largest Eigenvalue Convergence).

The sequence of largest eigenvalues $\left\{\lambda_{1}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}$ converges.

Proof.

To show convergence, we show that the sequence $\left\{\lambda_{1}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}$ is increasing and bounded above in $\mathbb{R}.$

First, we show the sequence is increasing. By induction on $N,$ we can write $\mathbf{A}_{N}$ in the block form

\displaystyle\mathbf{A}_{N}

\displaystyle=\begin{bmatrix}0&\mathbf{a}_{N}^{*}\\ \mathbf{a}_{N}&\mathbf{A}_{N-1}\end{bmatrix},

(4.46)

where $\mathbf{a}_{N}\in\mathbb{C}^{N-1}$ is the first column of $\mathbf{A}_{N}$ after excluding its first component $0$ . Notice that $\mathbf{A}_{N-1}$ is the first principal minor of $\mathbf{A}_{N},$ that is, we obtain $\mathbf{A}_{N-1}$ by deleting the first row and first column of $\mathbf{A}_{N}.$ By the Cauchy Interlacing Theorem [40], the eigenvalues of $\mathbf{A}_{N-1}$ interlace those of $\mathbf{A}_{N},$ which specifically means

\displaystyle\lambda_{N}(\mathbf{A}_{N})\leq\lambda_{N-1}(\mathbf{A}_{N-1})% \leq\lambda_{N-1}(\mathbf{A}_{N})\leq\dots\ \leq\lambda_{1}(\mathbf{A}_{N-1})% \leq\lambda_{1}(\mathbf{A}_{N})

(4.47)

The above inequality directly implies $\left\{\lambda_{1}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}$ is an increasing sequence.

Next, we show that $\left\{\lambda_{1}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}$ is a sequence that is bounded above. Consider $R_{N,N-j},$ the $(N-j)$ -th Gershgorin disk radius of $\mathbf{A}_{N}.$ Assuming Conjecture 4.1 holds, for $j>0,$ we have

	$\displaystyle\lim_{N\rightarrow\infty}R_{N,N-j}$	$\displaystyle=\frac{\binom{2j}{j}}{2^{2j-1}}$
		$\displaystyle=\frac{\frac{2j}{j}\binom{2(j-1)}{j-1}}{2^{2j-1}}$
		$\displaystyle=\frac{\binom{2(j-1)}{j-1}}{2^{2j-2}}$
		$\displaystyle<\frac{\binom{2(j-1)}{j-1}}{2^{2(j-1)-1}}=\lim_{N\rightarrow% \infty}R_{N,N-(j-1)}.$

Hence, as $N\rightarrow\infty,$ we see the Gershgorin radii of $\mathbf{A}_{N}$ are bounded above by $1.$ This implies the sequence $\left\{\lambda_{1}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}$ is bounded above by $1.$ ∎

Now, we make a conjecture about the limit of the sequence $\left\{\lambda_{n}(\mathbf{A}_{N})\right\}_{N\in\mathbb{N}}.$

Conjecture 4.2 (Asymptotic Binomial Matrix Largest Eigenvalue Convergence).

For small $n\in\mathbb{N},$ we have

\displaystyle\lim_{N\rightarrow\infty}\lambda_{2n-1}(\mathbf{A}_{N})

\displaystyle=\frac{2}{(2n-1)\pi}.

(4.48)

We provide numerical evidence supporting Conjecture 4.2. In Table 4.1, we tabulate the distance between $\lambda_{1}(\mathbf{A}_{N})$ and $\frac{2}{\pi}$ for $20\leq n\leq 24.$ We observe such distances are close to zero.

$N$	$\left\|\lambda_{1}(\mathbf{A}_{N})-\frac{2}{\pi}\right\|$
$20$	$\approx 0.0000232513$
$21$	$\approx 0.0000158971$
$22$	$\approx 0.0000109305$
$23$	$\approx 7.55596\times 10^{-6}$
$24$	$\approx 5.24969\times 10^{-6}$

Table 4.1: We tabulate the distance between the largest eigenvalue of

\mathbf{A}_{N}

and the conjectured limit as

N\rightarrow\infty

, namely

\frac{2}{\pi},

for

20\leq N\leq 24.

Now, we pose a conjecture about the largest eigenvalues corresponding to the principal minors of $\mathbf{A}_{N}.$ Recall the $j$ -th principal minor of $\mathbf{A}_{N}$ is the $(N-1)\times(N-1)$ matrix $\widetilde{\mathbf{A}}_{N,j}$ formed by deleting the $j$ -th row and $j$ -th column of $\mathbf{A}_{N}.$

Conjecture 4.3 (Asymptotic Binomial Matrix Principal Minor Largest Eigenvalues).

For each $1\leq j<N,$ we have

\displaystyle\lambda_{1}(\widetilde{\mathbf{A}}_{N,j+1})<\lambda_{1}(% \widetilde{\mathbf{A}}_{N,j})<\lambda_{1}(\mathbf{A}_{N}),

(4.49)

where $\lambda_{n}(\widetilde{\mathbf{A}}_{N,j})$ is the $n$ -th largest eigenvalue of $\widetilde{\mathbf{A}}_{N,j}.$

We provide some numerical evidence supporting the conjecture. In Table 4.2, we tabulate the largest eigenvalues of the principal minors of $\mathbf{A}_{N}$ in the situation $N=10$ . We observe as $j$ gets larger, the largest eigenvalue of the $j$ -th principal minor precipitously decreases.

$j$	$\lambda_{1}(\widetilde{\mathbf{A}}_{10,j})$
$1$	$\approx 0.634179$
$2$	$\approx 0.633523$
$3$	$\approx 0.632329$
$4$	$\approx 0.630125$
$5$	$\approx 0.626035$
$6$	$\approx 0.618269$
$7$	$\approx 0.6018$
$8$	$\approx 0.555167$
$9$	$\approx 0.347293$
$10$	$\approx 0.298889$

Table 4.2: We tabulate the largest eigenvalue for each principal minor of

\mathbf{A}_{10}.

4.3.2 Structure of Leading Eigenvector

For large $N,$ we now investigate the structure of $\mathbf{v}_{1}(\mathbf{A}_{N}),$ the leading eigenvector of $\mathbf{A}_{N},$ assuming all the previous conjectures we stated pertaining to the spectrum of $\mathbf{A}_{N}$ hold.

First, we describe the behavior of $|v_{1,n}(\mathbf{A}_{N})|,$ the modulus of the $n$ -th component $\mathbf{v}_{1}(\mathbf{A}_{N}).$

Theorem 4.9 (Asymptotic Binomial Matrix Leading Eigenvector Component Moduli Structure).

For each $1\leq n<N,$ we have

\displaystyle|v_{1,n}(\mathbf{A}_{N})|<|v_{1,n+1}(\mathbf{A}_{N})|.

(4.50)

Moreover, we have

\displaystyle\lim_{N\rightarrow\infty}v_{1,1}(\mathbf{A}_{N})

\displaystyle=0.

(4.51)

Proof.

We first prove the first statement. We utilize the eigenvalue-eigenvector identity [40] which states

\displaystyle\left|v_{1,j}(\mathbf{A}_{N})\right|^{2}

\displaystyle=\frac{\prod_{n=1}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda% _{n}(\widetilde{\mathbf{A}}_{N,j})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(% \mathbf{A}_{N})-\lambda_{n}(\mathbf{A}_{N})\right)}.

(4.52)

Assuming Conjecture 4.3 holds, we have

\displaystyle\lambda_{1}(\mathbf{A}_{N})-\lambda_{n}(\widetilde{\mathbf{A}}_{N% ,j})

\displaystyle<\lambda_{1}(\mathbf{A}_{N})-\lambda_{n}(\widetilde{\mathbf{A}}_{% N,j+1}).

This means we can bound the numerator on the right hand side of (4.52) in the following way:

	$\displaystyle\left\|v_{1,j}(\mathbf{A}_{N})\right\|^{2}$	$\displaystyle=\frac{\prod_{n=1}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda% _{n}(\widetilde{\mathbf{A}}_{N,j})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(% \mathbf{A}_{N})-\lambda_{n}(\mathbf{A}_{N})\right)}$
		$\displaystyle<\frac{\prod_{n=1}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda% _{n}(\widetilde{\mathbf{A}}_{N,j+1})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(% \mathbf{A}_{N})-\lambda_{n}(\mathbf{A}_{N})\right)}=\left\|v_{1,j+1}(\mathbf{A}% _{N})\right\|^{2},$

which proves the first statement.

Now, we prove the second statement, recall that $\lambda_{n}(\mathbf{A}_{N})\leq\lambda_{n}(\mathbf{A}_{N-1})$ by the Cauchy Interlacing inequality. We have

	$\displaystyle\|v_{1,1}(\mathbf{A}_{N})\|^{2}$	$\displaystyle=\frac{\prod_{n=1}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda% _{n}(\widetilde{\mathbf{A}}_{N,1})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(% \mathbf{A}_{N})-\lambda_{n}(\mathbf{A}_{N})\right)}$
		$\displaystyle=\frac{\prod_{n=1}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda% _{n}(\mathbf{A}_{N-1})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(\mathbf{A}_{N}% )-\lambda_{n}(\mathbf{A}_{N})\right)}$
		$\displaystyle=\frac{(\lambda_{1}(\mathbf{A}_{N})-\lambda_{1}(\mathbf{A}_{N-1})% )\ \prod_{n=2}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda_{n}(\mathbf{A}_{% N-1})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda_{n}(% \mathbf{A}_{N})\right)}$
		$\displaystyle\leq\frac{(\lambda_{1}(\mathbf{A}_{N})-\lambda_{1}(\mathbf{A}_{N-% 1}))\ \prod_{n=2}^{N-1}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda_{n}(\mathbf{A% }_{N})\right)}{\prod_{n=2}^{N}\left(\lambda_{1}(\mathbf{A}_{N})-\lambda_{n}(% \mathbf{A}_{N})\right)}$
		$\displaystyle=\frac{\lambda_{1}(\mathbf{A}_{N})-\lambda_{1}(\mathbf{A}_{N-1})}% {\lambda_{1}(\mathbf{A}_{N})-\lambda_{N}(\mathbf{A}_{N})}$
		$\displaystyle=\frac{\lambda_{1}(\mathbf{A}_{N})-\lambda_{1}(\mathbf{A}_{N-1})}% {2\lambda_{1}(\mathbf{A}_{N})},$

which approaches $0$ as $N\rightarrow\infty$ by Theorem 4.8. This proves the second statement. ∎

We interpret the previous theorem. The previous theorem states for large $N$ , the component moduli $\mathbf{v}_{1}(\mathbf{A}_{N})$ decrease from right to left. We interpret this in the context of our signal propagation model. As the signal propagates throughout the network, the left most sensors are less receptive to the signal than the right most sensors.

Now, we pose a conjecture about the component phases of $\mathbf{v}_{1}(\mathbf{A}_{N}).$

Conjecture 4.4 (Asymptotic Binomial Matrix Leading Eigenvector Component Phase Structure).

For large $N\in\mathbb{N},$ there exists large $K\leq N$ for which

0=\text{Arg}(v_{1,N}(\mathbf{A}_{N}))<\ \dots\ <\text{Arg}(v_{1,K}(\mathbf{A}_% {N}))<\pi.

In Table 4.3, we show numerical evidence of Conjecture 4.4 when $N=50.$ We observe the last $6$ component phases are in increasing order, supporting the conjecture. But at the current time of writing, we do not have a general conjecture pertaining to the behavior of the other phases.

We interpret the result of the conjecture. Suppose the signal in our propagation model were to be broadcast from the last sensor, namely the $N$ -th sensor. Then, Conjecture 4.4 states the order of the last few eigenvector component phases reflects the expected order of the next few sensors that receive the signal. Therefore, Cyclicity Analysis enables us to partially recover the expected network structure.

$n$	$\text{Arg}(v_{1,n}(\mathbf{A}_{50}))$
$1$	$\approx 2.11988$
$2$	$\approx 2.04186$
$3$	$\approx 1.96213$
$44$	$\approx-3.09417$
$45$	$\approx 2.93144$
$46$	$\approx 2.6331$
$47$	$\approx 2.27014$
$48$	$\approx 1.83497$
$49$	$\approx 1.38311$
$50$	$0$

Table 4.3: We tabulate some component phases of

\mathbf{v}_{1}(\mathbf{A}_{50}).

4.4 Only One Nonzero Propagation Coefficient and Independent and Identical Noise Injected into Few Sensors

Next, let $\mathbf{B}(\epsilon)$ be the friction matrix as defined in the beginning of Section 4.2. But this time, let $\boldsymbol{\Sigma}$ be a volatility matrix such that the diffusion matrix $\mathbf{D}$ is a diagonal diffusion matrix containing multiple nonzero diagonal entries that are equal. In particular, we let the last $L$ diagonal entries of $\mathbf{D}$ be positive and equal i.e. we assume $D_{N-L+1,N-L+1}=\ \dots\ =D_{N,N}=d$ for some $d>0,$ while all other diagonal entries are $0.$ This corresponds to the situation we inject independent and identical noise into the last $L$ sensors.

We explicitly derive the lead matrix with the assumptions imposed on $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma}.$

Theorem 4.10 (Cyclic OU Process One Nonzero Propagation Coefficient and Independent and Identical Noise in Multiple Sensors Lead Matrix Explicit Formula).

For each $\epsilon>0,$ the lead matrix $\mathbf{Q}(\epsilon)$ is explicitly

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\frac{b_{p}d}{N}\sum_{m,n=1}^{N}\sum_{\ell=1}^{L}\frac{(\omega_{% m}^{p-1}-\omega_{n}^{p-1})\ \omega_{m+n}^{\ell}}{2\epsilon+b_{p}(\omega_{m}^{p% -1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}

(4.53)

Proof.

Upon substituting $b_{1}=\epsilon-b_{p}$ and $b_{q}=0$ for all $q\neq 1,p$ into the right hand side of the Cyclic OU lead matrix formula in (4.6) and letting $\mathbf{e}_{n}$ be the $n$ -th canonical basis vector in $\mathbb{C}^{N},$ we have

	$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \mathbf% {w}_{N-m}^{\text{T}}\ \mathbf{D}\ \mathbf{w}_{N-n}\ \mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \mathbf% {w}_{N-m}^{\text{T}}\ \left(\frac{d\sum_{\ell=1}^{L}\omega_{N-n}^{(N-\ell+1)-1% }\ \mathbf{e}_{s_{\ell}}}{\sqrt{N}}\right)\ \mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\frac{b_{p}d}{N}\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-% \omega_{n}^{p-1})}{2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ % \mathbf{w}_{m}\ \ \left(\sum_{\ell=1}^{L}\omega_{N-m-n}^{-\ell}\right)\ % \mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\frac{b_{p}d}{N}\sum_{m,n=1}^{N}\sum_{\ell=1}^{L}\frac{(\omega_{% m}^{p-1}-\omega_{n}^{p-1})\ \omega_{m+n}^{\ell}}{2\epsilon+b_{p}(\omega_{m}^{p% -1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}.$

This proves the claim. ∎

We now investigate the structure of the leading eigenvector of $\mathbf{Q}(\epsilon)$ under the two regimes.

4.4.1 First Regime

We begin with the first regime when $\epsilon>0$ is large.

Theorem 4.11 (One Nonzero Propagation Coefficient and Independent and Identical Noise in Multiple Sensors First Regime Lead Matrix Asymptotic Expansion).

Consider the lead matrix $\mathbf{Q}(\epsilon)$ defined in (4.53). Then, as $\epsilon\rightarrow\infty,$ we have the following asymptotic expansion:

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\mathbf{A}(\epsilon)+O\left(\frac{1}{\epsilon^{2}}\right),

(4.54)

where $\mathbf{A}(\epsilon)$ is the $N\times N$ skew symmetric matrix whose $(j,k)$ -th entry is of the form

\displaystyle A_{j,k}(\epsilon)

\displaystyle=\sum_{\ell=1}^{L}\begin{cases}\frac{b_{p}d}{2\epsilon}&j\equiv% \ell+1-p\mod N\ ,\ k=\ell\\ -\frac{b_{p}d}{2\epsilon}&j=\ell\ ,\ k\equiv\ell+1-p\mod N\\ 0&\text{else}\end{cases}.

(4.55)

Proof.

Note that $\mathbf{Q}(\epsilon)$ is the sum of lead matrices of the form (4.6) with $s=N-\ell+1$ for each $1\leq\ell\leq L.$ Therefore, we substitute $s=N-\ell+1$ into the asymptotic formula (4.7), and the result immediately follows upon summing the $L$ total asymptotic formulas corresponding to each $1\leq\ell\leq L.$ ∎

Under the first regime, we investigate the leading eigenvector $\mathbf{v}_{1}(\mathbf{Q}(\epsilon))$ as $\epsilon\rightarrow\infty.$

At the current time of writing, we are not aware of any explicit formulas for the eigenvalues and eigenvectors of (4.54) for general $p.$ So we concentrate on the case $p=2$ i.e. only the propagation coefficient $b_{2}$ is nonzero. For simplicity, we let $b_{2}=-1$ and $d=1.$

Theorem 4.12 (First Nonzero Propagation Coefficient and Independent and Identical Noise in Multiple Sensors First Regime Leading Eigenvector Asymptotic Structure).

If $b_{2}=-1$ and $d=1,$ then the lead matrix $\mathbf{Q}(\epsilon)$ in (4.53) satisfies

\displaystyle\lim_{\epsilon\rightarrow\infty}\mathbf{v}_{1}(\mathbf{Q}(% \epsilon))

\displaystyle=\frac{\left(\underbrace{0\ ,\ \dots\ ,\ 0}_{\text{$N-L-1$ times}% }\ ,\ \left(-1\right)^{\frac{1}{2}}\sin\left(\frac{(1)\pi}{L+2}\right)\ ,\ % \dots\ ,\ \left(-1\right)^{\frac{L+1}{2}}\sin\left(\frac{(L+1)\pi}{L+2}\right)% \right)}{\sqrt{\sum_{\ell=1}^{L+1}\sin^{2}\left(\frac{\ell\pi}{L+2}\right)}}

(4.56)

Proof.

Let $\mathbf{A}(\epsilon)$ be the $N\times N$ skew-symmetric matrix whose $(j,k)$ -th entry is defined in (4.55) for each $1\leq j,k\leq N.$ Since

\displaystyle\|\mathbf{A}(\epsilon)-\mathbf{Q}(\epsilon)\|_{F}=O\left(\frac{1}% {\epsilon^{2}}\right),

(4.57)

we see that $\mathbf{A}(\epsilon)$ is a low rank approximation of $\mathbf{Q}(\epsilon)$ as $\epsilon\rightarrow\infty.$

We seek the leading eigenvector of $\mathbf{A}(\epsilon).$ Note that $\mathbf{A}(\epsilon)$ can be written in the block form

\displaystyle\mathbf{A}(\epsilon)

\displaystyle=\begin{bmatrix}\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{T}(\epsilon)\end{bmatrix},

(4.58)

in which $\mathbf{T}(\epsilon)$ is the $(L+1)\times(L+1)$ matrix such that $T_{\ell,\ell+1}(\epsilon)=-T_{\ell+1,\ell}=-\frac{b_{2}d}{2\epsilon}$ for each $1\leq\ell\leq L$ and all other entries are $0.$

We note that $\mathbf{T}(\epsilon)$ is a tridiagonal Toeplitz matrix, whose eigenvalues and eigenvectors are known [41]. In particular, if $\tau=T_{\ell,\ell+1},\sigma=T_{\ell+1,\ell}$ for each $1\leq\ell\leq L$ and $\delta=T_{\ell,\ell}$ for all $1\leq\ell\leq L+1.$ Then,

\displaystyle\widetilde{\lambda}_{\ell}

\displaystyle=\delta+2\sqrt{\sigma\tau}\cos\left(\frac{\ell\pi}{L+2}\right)

(4.59)

is an eigenvalue of $\mathbf{T}(\epsilon)$ with associated eigenvector

\displaystyle\widetilde{\mathbf{v}}_{\ell}

\displaystyle=\left(\left(\frac{\sigma}{\tau}\right)^{\frac{1}{2}}\sin\left(% \frac{(1)\ell\pi}{L+2}\right)\ ,\ \dots\ ,\ \left(\frac{\sigma}{\tau}\right)^{% \frac{L}{2}}\sin\left(\frac{(L)\ell\pi}{L+2}\right)\right).

(4.60)

Therefore, upon substituting $\tau=\frac{1}{2\epsilon}$ and $\sigma=-\tau$ and $\delta=0,$ the eigenvalues of the bottom right $(L+1)\times(L+1)$ submatrix $\mathbf{T}(\epsilon)$ are of the form

\displaystyle\widetilde{\lambda}_{\ell}

\displaystyle=-\frac{1}{\epsilon}\cos\left(\frac{\ell\pi}{L+2}\right)\ i

(4.61)

for each $1\leq\ell\leq L+1$ with corresponding eigenvector

\displaystyle\widetilde{\mathbf{v}}_{\ell}

\displaystyle=\left(\left(-1\right)^{\frac{1}{2}}\sin\left(\frac{(1)\ell\pi}{L% +2}\right)\ ,\ \dots\ ,\ \left(-1\right)^{\frac{L+1}{2}}\sin\left(\frac{(L+1)% \ell\pi}{L+2}\right)\right).

(4.62)

Therefore, the eigenvalues of $\mathbf{A}(\epsilon)$ are the eigenvalues of $\mathbf{T}(\epsilon)$ and $0.$ In particular, the largest eigenvalue of $\mathbf{A}(\epsilon)$ is $\widetilde{\lambda}_{1}$ with associated leading eigenvector being the vector formed by prepending $\widetilde{\mathbf{v}}_{1}$ with $N-L-1$ zeros. Hence, upon normalization, we have

\displaystyle\mathbf{v}_{1}(\mathbf{A}(\epsilon))

\displaystyle=\frac{\left(\underbrace{0\ ,\ \dots\ ,\ 0}_{\text{$N-L-1$ times}% }\ ,\ \left(-1\right)^{\frac{1}{2}}\sin\left(\frac{(1)\pi}{L+2}\right)\ ,\ % \dots\ ,\ \left(-1\right)^{\frac{L+1}{2}}\sin\left(\frac{(L+1)\pi}{L+2}\right)% \right)}{\sqrt{\sum_{\ell=1}^{L+1}\sin^{2}\left(\frac{\ell\pi}{L+2}\right)}},

(4.63)

which proves the claim. ∎

We interpret the result given by Theorem 4.12. The last $L+1$ leading eigenvector components are nonzero, suggesting the propagating signal in our model is received by the $L+1$ sensors. However, because the last $L+1$ components either lie on the real or imaginary axes the leading eigenvector component phases are either $0,\pi,$ or $\pm\frac{\pi}{2}.$ So we observe the cyclic order of the component phases does not match the expected cyclic order of sensors receiving the signal. Hence, Cyclicity Analysis does not enable us to fully recover the network structure.

However, we notice the last $L+1$ leading eigenvector component moduli form a sinusoidal pattern, which is interesting. This suggests if the signal were to be broadcast from the $N$ -th sensor, then one of the later sensors is most receptive to the signal.

4.4.2 Second Regime

We investigate the second regime, in which the perturbation $\epsilon>0$ is small.

Theorem 4.13 (One Nonzero Propagation Coefficient and Independent and Identical Noise in Multiple Sensors Second Regime Lead Matrix Asymptotic Expansion).

Consider the lead matrix $\mathbf{Q}(\epsilon)$ in (4.53). Then, as $\epsilon\rightarrow 0,$ the matrix $\mathbf{Q}(\epsilon)$ has an asymptotic expansion of the form

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=\frac{d}{N}\sum_{m,n=1}^{N}\sum_{\ell=1}^{L}\left(\frac{\left(% \omega_{m}^{p-1}-\omega_{n}^{p-1}\right)\ \omega_{m+n}^{1-\ell}}{\omega_{m}^{p% -1}+\omega_{n}^{p-1}-2}\right)\ \mathbf{w}_{m}\mathbf{w}_{n}^{\text{T}}+O(% \epsilon),

(4.64)

in which the summand inside the double sum is defined to be $0$ if $m=n=N.$

We omit the proof, as it involves the exact same reasoning that was used to prove Theorem 4.6. At the current time of writing, we are not aware of any closed formulas for the eigenvalues and eigenvectors of the matrix on the right hand side of (4.64). Nevertheless, we pose a conjecture in the situation where $p=2$ and $d=1.$

Conjecture 4.5 (First Nonzero Propagation Coefficient and Independent and Identical Noise in Multiple Sensors Second Regime Asymptotic Leading Eigenvector Structure).

Let $\mathbf{A}$ be the $N\times N$ matrix on the right hand side of (4.64). If $\mathbf{v}_{1}(\mathbf{A})$ is the leading eigenvector of $\mathbf{A},$ then we have

\displaystyle\text{Arg}(v_{1,N}(\mathbf{A}))<

\displaystyle\ \dots\ <\text{Arg}(v_{1,N-L+1}(\mathbf{A}))=0.

(4.65)

and $|v_{1,\ell}(\mathbf{A})|<|v_{1,\ell+1}(\mathbf{A})|$ for each $1\leq\ell\leq N-L$ and $|v_{1,\ell}(\mathbf{A})|>|v_{1,\ell+1}(\mathbf{A})|$ for each $N-L\leq\ell<N-1.$

We interpret this conjecture. In the situation where we inject independent and identical noise into $L$ nodes, the order of the few $L$ leading eigenvector component phases matches the expected order of those last few $L$ sensors receiving the signal in our propagation model. In addition, even though we inject noise into the last $L$ sensors, the $(N-L)$ -th sensor is most receptive to the propagating signal, which is interesting.

In Table 4.4, we provide evidence supporting Conjecture 4.5 in the situation where $N=10$ and $L=4.$ In particular, we observe the $7$ -th component has the largest modulus. Moreover, we observe the last $3$ component phases are increasing.

$n$	$\|v_{1,n}(\mathbf{A})\|$	$\text{Arg}(v_{1,n}(\mathbf{A}))$
$1$	$\approx 0.084589$	$\approx 2.97575$
$2$	$\approx 0.115083$	$\approx 2.61363$
$3$	$\approx 0.153119$	$\approx 2.22315$
$4$	$\approx 0.205092$	$\approx 1.80388$
$5$	$\approx 0.283689$	$\approx 1.3701$
$6$	$\approx 0.409898$	$\approx 0.94586$
$7$	$\approx 0.52438$	$0$
$8$	$\approx 0.466637$	$\approx-0.859917$
$9$	$\approx 0.354551$	$\approx-1.75523$
$10$	$\approx 0.217176$	$\approx-2.73113$

Table 4.4: Leading eigenvector component moduli and phases for the matrix

\mathbf{A}

defined in Conjecture 4.5 in the situation where

N=10

and

L=4.

4.5 Independent and Identical Noise Injected into All Sensors

We now investigate the signal propagation model under the assumptions the propagation coefficient $b_{p}$ is nonzero for exactly one index $1<p\leq N$ with $\gcd(p-1,N)=1$ and $\boldsymbol{\Sigma}$ is a volatility matrix such that the diffusion matrix $\mathbf{D}$ is proportional to the $N\times N$ identity matrix. This choice of $\boldsymbol{\Sigma}$ corresponds to the situation where we inject independent and identical noise into all $N$ sensors. Since the diagonal entries of $\mathbf{D}$ are all the same, we may assume without loss of generality that $\mathbf{D}$ is the $N\times N$ identity matrix itself.

We first derive the lead matrix with the imposed assumptions on these cyclic OU model parameters.

Theorem 4.14 (Cyclic OU Process One Nonzero Propagation Coefficient and Independent and Identical Noise in All Sensors Lead Matrix Explicit Formula).

For each $\epsilon>0,$ with the assumptions imposed on $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma},$ the lead matrix $\mathbf{Q}(\epsilon)$ is explicitly

\displaystyle\mathbf{Q}(\epsilon)

\displaystyle=-i\ \sum_{n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{\epsilon-2b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\ \mathbf{w}% _{N-n}\mathbf{w}_{n}^{\text{T}}.

(4.66)

Furthermore, $\mathbf{Q}(\epsilon)$ is a circulant matrix.

Proof.

We substitute the $N\times N$ identity matrix for $\mathbf{D}$ and $b_{1}=\epsilon-b_{p}$ and $b_{q}=0$ for all $q\neq 1,p$ into the general Cyclic OU lead matrix formula (4.4). Recalling $\left\{\mathbf{w}_{n}\right\}_{n=1}^{N}$ is an orthonormal basis of $\mathbb{C}^{N},$ we can simplify the summand in (4.4) in the following way:

$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \mathbf% {w}_{N-m}^{\text{T}}\mathbf{w}_{N-n}\mathbf{w}_{n}^{\text{T}}$
	$\displaystyle=\sum_{m,n=1}^{N}\frac{b_{p}(\omega_{m}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{m}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{m}\ \delta_% {m,N-n}\mathbf{w}_{n}^{\text{T}}$
	$\displaystyle=\sum_{n=1}^{N}\frac{b_{p}(\omega_{N-n}^{p-1}-\omega_{n}^{p-1})}{% 2\epsilon+b_{p}(\omega_{N-n}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{N-n}% \mathbf{w}_{n}^{\text{T}}$
	$\displaystyle=\sum_{n=1}^{N}\frac{b_{p}(\omega_{-n}^{p-1}-\omega_{n}^{p-1})}{2% \epsilon+b_{p}(\omega_{-n}^{p-1}+\omega_{n}^{p-1}-2)}\ \mathbf{w}_{N-n}\mathbf% {w}_{n}^{\text{T}}.$	(4.67)

Now, recall the Euler identities

	$\displaystyle\sin(\theta)$	$\displaystyle=\frac{e^{i\theta}-e^{-i\theta}}{2i}$
	$\displaystyle\cos(\theta)$	$\displaystyle=\frac{e^{i\theta}+e^{-i\theta}}{2}.$

Using the Euler identities, we rewrite the numerator and denominator of the summand in (4.67) to obtain

	$\displaystyle\mathbf{Q}(\epsilon)$	$\displaystyle=\sum_{n=1}^{N}\frac{-2ib_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{2\epsilon+b_{p}\left(2\cos\left(\frac{2\pi n(p-1)}{N}\right)-2\right)% }\ \mathbf{w}_{N-n}\mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=\sum_{n=1}^{N}\frac{-2ib_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{2\epsilon-4b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\ \mathbf{w% }_{N-n}\mathbf{w}_{n}^{\text{T}}$
		$\displaystyle=-i\ \sum_{n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{\epsilon-2b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\ \mathbf{w}% _{N-n}\mathbf{w}_{n}^{\text{T}}.$

Next, since $\mathbf{B}(\epsilon)$ is circulant, recall $\exp(-t\mathbf{B}(\epsilon))$ is circulant for each $t\geq 0$ by Theorem 2.3. Examining the explicit formula for the stationary covariance matrix $\mathbf{S}$ in (3.5), we see that the integrand is a product of the circulant matrices $\exp(-t\mathbf{B}(\epsilon))$ and its transpose, which implies $\mathbf{S}$ is circulant. Therefore, $\mathbf{B}(\epsilon)\ \mathbf{S}$ and $\mathbf{S}\ \mathbf{B}^{\text{T}}(\epsilon)$ are both circulant, and since $\mathbf{Q}(\epsilon)$ is proportional to the difference between $\mathbf{B}(\epsilon)\ \mathbf{S}$ and $\mathbf{S}\ \mathbf{B}^{\text{T}}(\epsilon),$ we conclude $\mathbf{Q}(\epsilon)$ is circulant. ∎

4.5.1 Second Regime

We now consider the second regime in which the perturbation constant $\epsilon>0$ is small.

Theorem 4.15 (Cyclic OU Process One Nonzero Propagation Coefficient and Independent and Identical Noise in All Sensors Second Regime Leading Eigenvector Explicit Formula).

Consider the lead matrix $\mathbf{Q}(\epsilon)$ defined in (4.66). Then, as $\epsilon\rightarrow 0,$ the leading eigenvector of $\mathbf{Q}(\epsilon)$ is $\mathbf{w}_{q},$ where $1\leq q\leq N$ satisfies $q(p-1)\equiv 1\mod N.$

Proof.

First, we need to determine the eigenvalues of $\mathbf{Q}(\epsilon).$ In order to do this, we need to compute the first row of $\mathbf{Q}(\epsilon).$ Since $\mathbf{Q}(\epsilon)$ is skew-symmetric, its first row is the opposite of its first column. To obtain its first column, we compute the product of $\mathbf{Q}(\epsilon)$ and $\mathbf{e}_{1},$ where $\mathbf{e}_{1}$ is the first standard basis vector in $\mathbb{C}^{N}.$ This means the first row of $\mathbf{Q}(\epsilon)$ is

	$\displaystyle-\mathbf{Q}(\epsilon)\ \mathbf{e}_{1}$	$\displaystyle=i\ \sum_{n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{\epsilon-2b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\ \mathbf{w}% _{N-n}\mathbf{w}_{n}^{\text{T}}\ \mathbf{e}_{1}$
		$\displaystyle=\frac{i}{\sqrt{N}}\sum_{n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n% (p-1)}{N}\right)}{\epsilon-2b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\ % \mathbf{w}_{N-n}.$		(4.68)

As a result, for each $1\leq q\leq N,$ each eigenvalue of $\mathbf{Q}(\epsilon)$ is of the form

	$\displaystyle\widetilde{\lambda_{q}}(\epsilon)$	$\displaystyle=\sqrt{N}\ \mathbf{w}^{\text{T}}_{q}\left(\frac{i}{\sqrt{N}}\sum_% {n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n(p-1)}{N}\right)}{\epsilon-2b_{p}% \sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\ \mathbf{w}_{N-n}\right)$
		$\displaystyle=i\ \sum_{n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{\epsilon-2b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\mathbf{w}_{% q}^{\text{T}}\ \mathbf{w}_{N-n}$
		$\displaystyle=i\ \sum_{n=1}^{N}\frac{b_{p}\sin\left(\frac{2\pi n(p-1)}{N}% \right)}{\epsilon-2b_{p}\sin^{2}\left(\frac{\pi n(p-1)}{N}\right)}\delta_{N-q,% N-n}$
		$\displaystyle=\frac{b_{p}\sin\left(\frac{2\pi q(p-1)}{N}\right)}{\epsilon-2b_{% p}\sin^{2}\left(\frac{\pi q(p-1)}{N}\right)}\ i.$

Note that

$\displaystyle\lim_{\epsilon\rightarrow 0}\left\|\widetilde{\lambda_{q}}(% \epsilon)\right\|$	$\displaystyle=\left\|\frac{b_{p}\sin\left(\frac{2\pi q(p-1)}{N}\right)}{-2b_{p}% \sin^{2}\left(\frac{\pi q(p-1)}{N}\right)}\right\|$
	$\displaystyle=\left\|\frac{\sin\left(\frac{\pi q(p-1)}{N}\right)\cos\left(\frac% {\pi q(p-1)}{N}\right)}{\sin^{2}\left(\frac{\pi q(p-1)}{N}\right)}\right\|$
	$\displaystyle=\left\|\cot\left(\frac{\pi q(p-1)}{N}\right)\right\|.$	(4.69)

We seek the index $q$ maximizing $\left|\cot\left(\frac{\pi q(p-1)}{N}\right)\right|.$ Note that the map $\theta\mapsto|\cot(\theta)|$ is strictly decreasing on $(0,\frac{\pi}{2})$ and strictly increasing on $\left(\frac{\pi}{2},\pi\right).$ Since $\cot(\theta)=\cot\left(\frac{\pi}{2}-\theta\right)$ for all $\theta\in(0,\frac{\pi}{2}),$ the map $\theta\mapsto|\cot(\theta)|$ restricted to the interval $(0,\pi)$ is symmetric about $\frac{\pi}{2}.$ Finally, note that $\theta\mapsto|\cot(\theta)|$ is $\pi$ -periodic. The integers $u\in\mathbb{Z}$ maximizing $\left|\cot(\frac{u\pi}{N})\right|$ specifically in the domain $(0,\pi)$ are $u=1$ and $u=N-1.$ By periodicity, the integer $u\in\mathbb{Z}$ maximizing $\left|\cot(\frac{u\pi}{N})\right|$ on its entire domain satisfies $u\equiv\pm 1\mod N.$

This means an index $1\leq q\leq N$ that maximizes the eigenvalue limit (4.66) must satisfy either the congruence equation $q(p-1)\equiv 1\mod N$ or the congruence equation $q(p-1)\equiv-1\mod N.$ Both congruence equations have integer solutions because $\gcd(p-1,N),$ which is equal to $1$ by assumption, divides $\pm 1.$

Let $1\leq q\leq N$ be such that $q(p-1)\equiv 1\mod N.$ Then, the associated eigenvector of the eigenvalue $\widetilde{\lambda}_{q}(\epsilon)$ is $\mathbf{w}_{q},$ which is the leading eigenvector of $\mathbf{Q}(\epsilon).$ ∎

Knowing the leading eigenvector of $\mathbf{Q}(\epsilon)$ under the second regime, we investigate whether the structure of the leading eigenvector recovers the structure of the cyclic sensor network.

Theorem 4.16 (Cyclic OU Process One Nonzero Propagation Coefficient and Independent and Identical Noise in All Sensors Second Regime Leading Eigenvector Structure).

Let $\mathbf{w}_{q}$ be the leading eigenvector of $\mathbf{Q}(\epsilon).$ Then, the permutation $\sigma\in S_{N}$ representing the cyclic order of the component phases $\mathbf{w}_{q}$ is defined by $\sigma(n)\equiv(N-(n-1)(p-1))\mod N.$

Proof.

First, we note that $\gcd(q,N)=1.$ To see this, recall the previous theorem states that $q(p-1)\equiv 1\mod N.$ By definition, this means

q(p-1)-Nk=1

for some integer $k.$ Hence, we can write $1$ as an integer linear combination of the numbers $q$ and $N,$ which implies $\gcd(q,N)=1.$

Let $\sigma\in S_{N}$ be the permutation such that

\displaystyle q(\sigma(n)-1)

\displaystyle=n-1\mod N.

(4.70)

for each $1\leq n\leq N.$ Firstly, we note $\sigma$ is well-defined. For each $1\leq n\leq N,$ because $\gcd(q,N)=1$ divides $n-1,$ there is a unique solution to the congruence equation (4.70) for $\sigma(n).$

For each $1\leq n<N,$ we have

	$\displaystyle q(\sigma(n+1)-\sigma(n))$	$\displaystyle\equiv(n+1-n)\mod N$
		$\displaystyle\equiv 1\mod N.$

Therefore, we deduce

\displaystyle q(\sigma(n+1)-\sigma(n))\equiv q(p-1)\mod N.

Because $\gcd(q,N)=1,$ we may divide both sides of the congruence equation by $q$ to deduce that

\displaystyle\sigma(n+1)-\sigma(n)

\displaystyle\equiv p-1\mod N.

(4.71)

Finally, we claim $\sigma$ is a cyclic order permutation corresponding to the component principal arguments of the eigenvector $\mathbf{w}_{q}.$ Recall the $n$ -th component of $\mathbf{w}_{q}$ is explicitly

	$\displaystyle w_{q,n}$	$\displaystyle=\frac{\omega_{q}^{n-1}}{\sqrt{N}}$
		$\displaystyle=\frac{\omega_{1}^{q(n-1)}}{\sqrt{N}}.$

By (4.70), we have

	$\displaystyle w_{q,\sigma(n)}$	$\displaystyle=\frac{\omega_{1}^{q(\sigma(n)-1)}}{\sqrt{N}}$
		$\displaystyle=\frac{\omega_{1}^{n-1}}{\sqrt{N}},$

which has principal argument $\frac{2\pi(n-1)}{N}\in[0,2\pi).$ Therefore, our permutation $\sigma$ ensures the components of $\mathbf{w}_{q}$ are sorted in increasing order by their corresponding principal arguments. Thus, $\sigma$ is a valid cyclic order permutation. ∎

We interpret the result of the previous theorem. If we inject independent and identical noise into all sensors, then under the second regime, the structure of the leading eigenvector of $\mathbf{Q}(\epsilon)$ enables us to recover the expected network structure induced by the friction matrix $\mathbf{B}(\epsilon).$ In particular, the cyclic order of the phases of the leading eigenvector matches the expected cyclic order of the sensors receiving the signal in our propagation model.

Chapter 5 Experimental Results

In this chapter, we discuss the experimental results pertaining to the main topic of our thesis. Throughout, we consider an $N$ -dimensional OU Process $\left\{\mathbf{x}(t)\right\}_{t\geq 0}$ with fixed model parameters $\mathbf{B}\in\text{Mat}_{N,N}(\mathbb{R})$ and $\boldsymbol{\Sigma}\in\text{Mat}_{M,N}(\mathbb{R}).$ We let $\mathbf{D}$ be the diffusion matrix and $\mathbf{Q}$ be the lead matrix of the OU process. For the purposes of this chapter, we specifically refer to $\mathbf{Q}$ as the theoretical lead matrix of the OU process.

5.1 OU Process Realization Time Series Simulation

In this section, we describe how to generate a time series approximating a realization of the OU process over a time interval of the form $[0,T],$ where $T>0$ is a fixed, large constant. We employ the Euler-Maruyama method [42], which we describe in detail.

First, we fix a large integer $K,$ representing the number of iterations, and a small constant $\Delta>0,$ representing the time step size between consecutive times. In the original governing SDE (3.1) of the OU process, we replace the time differential $dt$ with $\Delta,$ the stochastic differential $d\mathbf{x}(t)$ with the difference $\mathbf{x}(t+\Delta)-\mathbf{x}(t),$ and the Wiener differential $d\mathbf{w}(t)$ with the difference $\mathbf{w}(t+\Delta)-\mathbf{w}(t).$

Recall the increments of the standard Wiener process are independent and Gaussian; in particular, $\mathbf{w}(t+\Delta)-\mathbf{w}(t)$ is an $M$ -dimensional Gaussian random vector with mean $\mathbf{0}$ and covariance matrix $\Delta\mathbf{I},$ where $\mathbf{I}$ is the $M\times M$ identity matrix. So we generate a sequence of vectors $\left\{\boldsymbol{\xi}_{k}\right\}_{k=0}^{K-1}$ representing the Wiener increments, in which all vectors are randomly and independently chosen according to the $M$ -dimensional Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\Delta\mathbf{I}.$ We let $\mathbf{x}_{0}\in\mathbb{R}^{N}$ be an initial starting vector randomly chosen according to the Gaussian distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{S},$ where $\mathbf{S}$ solves the Lyapunov equation (3.2).

We now construct the time series $\left\{\mathbf{x}_{k}\right\}_{k=0}^{K-1},$ where $\mathbf{x}_{k}=\mathbf{x}(k\Delta)$ for each $0\leq k<K.$ Explicitly, we have

\displaystyle\mathbf{x}_{k+1}

\displaystyle=\mathbf{x}_{k}-\Delta\mathbf{B}\ \mathbf{x}_{k}+\boldsymbol{% \Sigma}\ \boldsymbol{\xi}_{k}

(5.1)

for each $0\leq k<K-1.$ This constructed time series $\left\{\mathbf{x}_{k}\right\}_{k=0}^{K-1}$ is what we use to approximate a realization of the OU process over the specific time interval $\left[0\ ,\ (K-1)\ \Delta\right].$

In Figure 5.1, we plot a sample time series approximating a realization of a $2$ -dimensional (cyclic) OU process over the finite time interval $[0,100]$ with model parameters $\mathbf{B}=\begin{bmatrix}1.2&-0.2\\ -0.2&1.2\end{bmatrix}$ and $\boldsymbol{\Sigma}$ equal to the $2\times 2$ identity matrix. To generate the time series, we used the scheme (5.1) with $K=10001$ iterations and time step size $\Delta=0.01.$

5.2 OU Process Empirical Lead Matrix Generation

Throughout, we fix $K\in\mathbb{N}$ and $\Delta>0.$ Consider a time series $\left\{\mathbf{x}_{k}\right\}_{k=0}^{K-1}$ produced according to our scheme in (5.1). We compute its empirical lead matrix, denoted $\mathbf{A}(K,\Delta),$ via the shoelace formula mentioned in (1.22).

As a consequence of the Strong Law of Large numbers identity stated in Theorem 3.7, for each fixed $\Delta>0,$ we have

\displaystyle\lim_{K\rightarrow\infty}\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}

\displaystyle=\mathbf{Q},

(5.2)

for almost all realizations of the OU process, where $\mathbf{Q}$ is the theoretical lead matrix.

We show numerical evidence of the discrete law of large numbers identity (5.2) via an example. Consider the $5$ -dimensional OU process with model parameters $\mathbf{B}=\text{Circ}(2.1,-0.2,-0.4,-0.6,-0.8)$ and $\boldsymbol{\Sigma}$ equal to the $5\times 5$ identity matrix. In Figure 5.2, for each $K\in\left\{10^{2}+1\ ,\ \dots\ ,\ 10^{6}+1\right\}$ and each $\Delta\in\left\{10^{-1}\ ,\ 10^{-2}\right\},$ we display a scatterplot, in which the first coordinate is the $(m,n)$ -th entry of the time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}$ and the second coordinate is the $(m,n)$ -th entry of the theoretical lead matrix $\mathbf{Q}.$ We observe as $K$ increases, the scatterplot approaches the straight line $y=x,$ which means the entries of the time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}$ approach the corresponding entries of the theoretical matrix $\mathbf{Q}$ . Furthermore, our example shows that even when decreasing the time step size from $10^{-1}$ to $10^{-2}$ , we still would need a large number of iterations $K$ if we want an accurate estimate of $\mathbf{Q}$ from a time-averaged empirical lead matrix.

5.3 Cyclicity Analysis of the Cyclic OU Process

We revisit the main problem statement of our thesis, in which we consider signal propagation model governed by the cyclic OU process with model parameters $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma}$ for each $\epsilon>0,$ where the circulant friction matrix $\mathbf{B}(\epsilon)$ is of the form

\displaystyle\mathbf{B}(\epsilon)

\displaystyle=\text{Circ}\left(\epsilon-\sum_{p=2}^{N}b_{p}\ ,\ b_{2}\ ,\ % \dots\ ,\ b_{N}\right).

for some fixed propagation coefficients $b_{2}\ ,\ \dots\ ,\ b_{N}\leq 0.$

Throughout this section, we fix the dimension $N=100$ assume the first propagation coefficient $b_{2}$ is equal to $-1$ while all other propagation coefficients are $0.$ This means the expected network structure is such that the $n$ -th sensor is receptive only to activity within the $(n+1)$ -th sensor, its immediate neighbor to the right, in which $n+1$ is indexed mod $N$ . We also assume $\boldsymbol{\Sigma}$ is a volatility matrix such that the diffusion matrix $\mathbf{D}$ is a diagonal matrix whose positive entries are equal. We let $d_{s}\in\{0,1\}$ be the $s$ -th diagonal entry of $\mathbf{D}.$ Recall $d_{s}=1$ corresponds to the situation where noise is injected into the $s$ -th sensor. For different choices of $\boldsymbol{\Sigma}$ and perturbation $\epsilon>0,$ in which $\epsilon$ is some power of $10,$ we produce a time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}$ corresponding to a realization of the $100$ -dimensional cyclic OU process with model parameters $\mathbf{B}(\epsilon)$ and $\boldsymbol{\Sigma}.$ Here, we fix $K=10^{6}+1$ iterations and $\Delta=0.01$ as the time step size. We present numerical results as to whether the leading eigenvector of the empirical lead matrix enables us to recover the structure of the cyclic network structure induced by $\mathbf{B}(\epsilon).$

5.3.1 Injection of Noise into Only One Sensor

Suppose we inject noise into only one sensor. Here, we consider the situation where we inject noise into the last sensor i.e. $d_{100}=1$ and $d_{s}=0$ for $1\leq s<100.$ We let $\epsilon\in\left\{10^{-11},10^{-10},10^{-1},10^{0},10^{3},10^{4}\right\}.$ We chose these particular values of $\epsilon$ for a specific reason. For $\epsilon<10^{-11},$ the friction matrix $\mathbf{B}(\epsilon)$ is considered unstable due to precision issues. For $\epsilon>10^{4},$ there was no change to the structure of the leading eigenvector. For the other listed values of $\epsilon,$ we noticed intermediary changes within the structure of the leading eigenvector.

In Figure 5.3, for each such value of $\epsilon,$ we plot the heatmap of a time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}(\epsilon)$ and the logarithms of its eigenvalue moduli in descending order. We also plot the component phases and moduli of the leading eigenvector. In Table 5.1, we record the ratio of the largest to third largest eigenvalue for each such empirical lead matrix.

$\epsilon$	$\left\|\lambda_{1}/\lambda_{3}\right\|$
$10^{-11}$	$\approx 2.00022$
$10^{-10}$	$\approx 2.00022$
$10^{-1}$	$\approx 8.00132$
$1$	$\approx 201.8365$
$10^{3}$	$\approx 16064344115142$
$10^{4}$	$\approx 2.24629\times 10^{16}$

Table 5.1: The ratio of the first to third largest eigenvalue of each time-averaged empirical lead matrix in Figure 5.3.

5.3.2 Observations

For small $\epsilon,$ we observe the components of the leading eigenvector spiral around the origin. We observe that the last few component phases are decreasing. This means if the signal in our propagation model is broadcast from the last sensor, then the leading eigenvector correctly detects the next few sensors that receive the signal as it propagates. However, we do not see overall decreasing behavior within the component phases corresponding to other indices. This means the global order of the component phases does not reflect the expected order of sensors receiving the propagating signal. Therefore, the phases alone suggest we cannot fully recover the overall network structure. Moreover, we observe all but the last few leading eigenvector component moduli are close to $0,$ which suggests as the signal is propagating, all but the last few sensors are receptive to the signal. Finally, we observe there is no clear separation between the largest eigenvalue and third largest eigenvalue of the empirical lead matrix, which suggests the largest eigenvalue does not actually dominate the spectrum of the lead matrix.

As $\epsilon$ gets larger, the spiraling pattern among the leading eigenvector components begins to disappear. The last leading eigenvector component stays put on the real axis. But the second to last leading eigenvector component gets closer to the imaginary axis, while all other leading eigenvector components get closer to the origin. We see chaotic behavior in all but the last two component phases. This is because while the last two leading eigenvector components approach the origin, they do so from different directions. This suggests that only the last two sensors are receptive to the propagating signal, which means we cannot fully recover the expected network structure. Moreover, there is a noticeable gap between the largest eigenvalue and the third largest eigenvalue of the empirical lead matrix, which means the empirical lead matrix can be well-approximated via a low rank matrix.

5.3.3 Independent, Identically Distributed Noise Injected into Few Sensors

Now, suppose we inject independent and identical noise into only a few sensors. Consider the situation where we inject noise into the last $11$ sensors i.e. $d_{s}=1$ for all $90\leq s\leq 100$ and $d_{s}=0$ for all other $s.$ We let $\epsilon\in\left\{10^{-10},10^{-9},10^{-3},10^{-2},1,10^{2},10^{3}\right\}.$ In Figure 5.4, we plot the heatmap of a time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}(\epsilon)$ and the logarithms of its eigenvalue moduli in descending order. We also plot the component phases and moduli of the leading eigenvector. In Table 5.2, we record the ratio of the largest to third largest eigenvalue of each empirical lead matrix.

$\epsilon$	$\left\|\lambda_{1}/\lambda_{3}\right\|$
$10^{-10}$	$\approx 1.2063$
$10^{-9}$	$\approx 1.2063$
$10^{-3}$	$\approx 1.3219$
$10^{-2}$	$\approx 1.5622$
$1$	$\approx 1.1248$
$10^{2}$	$\approx 1.0965$
$10^{3}$	$\approx 0.9119$

Table 5.2: We tabulate the ratios between the first and third largest eigenvalues corresponding to each time-averaged empirical lead matrix in Figure 5.4.

5.3.4 Observations

For small $\epsilon,$ the components of the leading eigenvector exhibit a visible spiral pattern. We notice, however, the last few leading eigenvector component phases are not decreasing. This suggests if the signal in our propagation model is broadcast from the $100$ -th sensor, then we are not able to recover the correct order of the next few sensors receiving the signal. Moreover, we observe that the $90$ -th leading eigenvector component modulus is the largest, which suggests the $90$ -th sensor is most receptive to the propagating signal. Recall that the index $90$ is the smallest index corresponding to the sensor we inject noise. Despite the fact we inject noise into the $s$ -th sensor for $90\leq s\leq 100,$ the corresponding leading eigenvector component moduli suggest the other sensors receiving noise are not as receptive to the propagating signal. Finally, we observe there is no clear separation between the largest eigenvalue of the empirical lead matrix and third largest eigenvalue.

As $\epsilon$ gets larger, we see the spiraling pattern within the leading eigenvector components disappears. The leading eigenvector components eventually lie on either the real axis, imaginary axis, or the origin.

5.3.5 Independent, Identically Distributed Noise Injected into a Large Number of Sensors

Suppose we inject noise into a larger number of sensors. In particular, suppose we inject noise into the last $51$ sensors i.e. $d_{s}=1$ for all $50\leq s\leq 100$ and $d_{s}=0$ for all other $s.$ We let $\epsilon\in\left\{10^{-9},10^{-8},10^{-3},10^{-2},1,10^{2},10^{3},10^{4}\right\}.$ In Figure 5.5, we plot the heatmap of a time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}(\epsilon)$ and the logarithms of its eigenvalue moduli in descending order. We also plot the component phases and moduli of the leading eigenvector. In Table 5.3, we record the ratio of the largest to third largest eigenvalue for each such empirical lead matrix.

$\epsilon$	$\left\|\lambda_{1}/\lambda_{3}\right\|$
$10^{-9}$	$\approx 1.8548$
$10^{-8}$	$\approx 1.8548$
$10^{-3}$	$\approx 1.387$
$10^{-2}$	$\approx 1.194$
$1$	$\approx 1.007$
$10^{2}$	$\approx 1.005$
$10^{3}$	$\approx 1.005$
$10^{4}$	$\approx 1.005$

Table 5.3: We tabulate the ratios between the first and third largest eigenvalues corresponding to each time-averaged empirical lead matrix in Figure 5.5.

5.3.6 Observations

For small $\epsilon,$ we observe the components of the leading eigenvector appear to form an elliptical pattern. The component principal arguments exhibit linear monotonically decreasing behavior. This suggests the leading eigenvector phases do approximately reflect the expected structure of the network. Moreover, the component moduli are nonzero and appear to exhibit a symmetric pattern. Finally, we observe there is no clear separation between the first and third largest eigenvalues of the lead matrix.

When $\epsilon$ gets larger, the elliptical pattern disappears, and the leading eigenvector components approach either the real axis, imaginary axis, or the origin. The last $50$ leading eigenvector component moduli are large, while the others are close to $0.$ Therefore, Cyclicity Analysis suggests we cannot fully recover the structure of the network.

5.3.7 Independent, Identically Distributed Noise Injected into All Sensors

Now, suppose we inject independent and identical noise into all sensors. This corresponds to the situation where $\boldsymbol{\Sigma}$ is a volatility matrix such that the diffusion matrix $\mathbf{D}$ is the identity matrix.

In Figure 5.6, we plot the heatmap of a time-averaged empirical lead matrix $\frac{\mathbf{A}(K,\Delta)}{(K-1)\Delta}(\epsilon)$ for each $\epsilon\in\left\{10^{-9},10^{-8},10^{-2},1,10,10^{3},10^{4}\right\}$ and the logarithms of its eigenvalue moduli in descending order. We also plot the component phases and moduli of the leading eigenvector. In Table 5.3, we record the ratio of the largest to third largest eigenvalue for each such empirical lead matrix.

$\epsilon$	$\left\|\lambda_{1}/\lambda_{3}\right\|$
$10^{-9}$	$\approx 2.001$
$10^{-8}$	$\approx 2.001$
$10^{-2}$	$\approx 1.036$
$1$	$\approx 1.0009$
$10$	$\approx 1.0002$
$10^{3}$	$\approx 1.001$
$10^{4}$	$\approx 1.001$

Table 5.4: We tabulate the ratios between the first and third largest eigenvalues corresponding to each time-averaged empirical lead matrix in Figure 5.6.

5.3.8 Observations

For small $\epsilon,$ the components of the leading eigenvector lie on a circle, as indicated by their moduli being equal. Observe the cyclic order of the leading eigenvector component phases correctly matches the order of sensors receiving the signal in our model. Therefore, Cyclicity Analysis enables us to fully recover the network structure.

For large $\epsilon,$ the components of the leading eigenvector still lie on a circle, but they lie on either the positive real, positive imaginary, negative real, or the negative imaginary axis. Therefore, the leading eigenvector component phases do not enable us to recover the network structure.

Chapter 6 Conclusion

In this thesis, we posed very general questions pertaining to the lead lag dynamics amongst the components of an $N$ -dimensional signal. Firstly, if two signals evolve similarly throughout time, which signal leads and which signal follows ? Secondly, if all signals evolve similarly throughout time, what is the order in which such signals evolve ?

We discussed Cyclicity Analysis, an established way of answering the aforementioned two questions for cyclic signals, which are repeating, yet aperiodic, signals. Cyclicity Analysis is different from many standard tools in that it is time reparameterization invariant; the results of the analysis do not depend on how one observes the signals throughout time. It involves two parts. The first part is to construct an $N\times N$ lead matrix, in which the sign of the $(m,n)$ -th entry heuristically captures the pairwise leader follower relationship between the $m$ -th and $n$ -th component signals. It is a specific construction based on iterated path integrals, which are the functionals of the signal $\mathbf{x}$ that are time reparameterization invariant. The second part of Cyclicity Analysis is to assume a model in which the component signals of $\mathbf{x}$ trace some periodic signal up to time reparameterization, scaling constants, and time shifts. Under this model, one determines the order in which the component signals evolve throughout time, which is the cyclic order of those time shifts. Under certain baseline situations, the leading eigenvector of the lead matrix can approximately recover the order in which the component signals of $\mathbf{x}$ evolve throughout time. More specifically, the cyclic order of the eigenvector component phases reflects the cyclic order of the offsets.

In this thesis, we investigated Cyclicity Analysis in a novel setting. We considered an $N$ -dimensional Ornstein-Uhlenbeck (OU) stochastic process with model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ in which $\mathbf{B}$ is a friction matrix and $\boldsymbol{\Sigma}$ is a volatility matrix. We considered a signal propagation model governed by the cyclic OU process, in which the model parameter $\mathbf{B}$ is a circulant friction matrix. In this model, there is a signal propagating throughout a network of $N$ sensors. The matrix $\mathbf{B}$ represents the cyclic network structure of the sensors, while the diffusion matrix $\mathbf{D}=\boldsymbol{\Sigma}\boldsymbol{\Sigma}^{\text{T}}/2$ is the covariance matrix whose entries are pairwise covariances of noises injected into any two sensors.

To make Cyclicity Analysis suitable for analyzing the OU process, we defined an auxiliary $N\times N$ lead process, a matrix-valued stochastic process in which each random matrix in the collection represents a lead matrix corresponding to a realization of the OU process. The first main result of the thesis is that the lead process obeys a strong law of large numbers identity: its time average converges almost surely to a skew-symmetric matrix $\mathbf{Q}$ .

Then, with different choices of cyclic OU model parameters $\mathbf{B}$ and $\boldsymbol{\Sigma},$ we investigated whether Cyclicity Analysis could enable us to recover the network structure induced by $\mathbf{B}$ under our signal propagation model. More specifically, we analyzed whether the structure of the leading eigenvector of $\mathbf{Q}$ could enable us to recover the network structure induced by $\mathbf{B}.$

For our analysis, we assumed the ground truth network structure is such that every sensor receives the propagating signal from exactly one sensor to its right. In addition, we assumed independent and identical noise is injected into multiple sensors, which corresponds to the situation where $\mathbf{D}$ is a diagonal diffusion matrix whose nonzero diagonal entries are equal.

With the imposed assumptions on $\mathbf{B}$ and $\boldsymbol{\Sigma},$ our experimental results showed that Cyclicity Analysis was in fact able to partially recover the structure of the network. For example, in the situation where we inject noise only into the last node, the cyclic order of the component phases of the leading eigenvector of $\mathbf{Q}$ partially reflected the expected cyclic order of the sensors receiving the signal in our model. As we increase the number of sensors receiving noise, we observed the cyclic order of the component phases approximately match the expected cyclic order of the sensors. In the ultimate case identical noise is injected into all sensors, the cyclic orders perfectly match.

There is future work that can be done. For example, what if we remove the assumption that $\mathbf{D}$ has equal positive diagonal entries ? This would correspond to a situation where independent, but not identical, noise is injected into multiple sensors. What if $\mathbf{B}$ represents a network structure such that each sensor does not receive the signal from only one other sensor, but multiple sensors ? Does the leading eigenvector of the lead matrix $\mathbf{Q}$ capture those phenomena ?

References

[1] Elias M Stein and Rami Shakarchi. Fourier analysis: an introduction, volume 1. Princeton University Press, 2011.
[2] Athanasios Papoulis. The Fourier integral and its applications. McGraw-Hill New York, 1967.
[3] Yuliy Baryshnikov and Emily Schlafly. Cyclicity in multivariate time series and applications to functional mri data. In 2016 IEEE 55th conference on decision and control (CDC), pages 1625–1630. IEEE, 2016.
[4] Ivan Abraham, Somayeh Shahsavarani, Benjamin Zimmerman, Fatima Husain, and Yuliy Baryshnikov. Slow cortical waves through cyclicity analysis. bioRxiv, 2021.
[5] Kuo-Tsai Chen. Integration of paths–a faithful representation of paths by noncommutative formal power series. Transactions of the American Mathematical Society, 89(2):395–407, 1958.
[6] Kuo-Tsai Chen. Algebras of iterated path integrals and fundamental groups. Transactions of the American Mathematical Society, 156:359–379, 1971.
[7] Kuo-Tsai Chen. Iterated integrals of differential forms and loop space homology. Annals of Mathematics, 97(2):217–246, 1973.
[8] Kuo-Tsai Chen. Iterated path integrals. Bulletin of the American Mathematical Society, 83(5):831–879, 1977.
[9] Terry J Lyons. Differential equations driven by rough signals. Revista Matemática Iberoamericana, 14(2):215–310, 1998.
[10] Terry J Lyons, Michael Caruana, and Thierry Lévy. Differential equations driven by rough paths. Springer, 2007.
[11] Sheldon Axler. Linear algebra done right. Springer Nature, 2024.
[12] Eugene L Allgower and Phillip H Schmidt. Computing volumes of polyhedra. mathematics of computation, 46(173):171–174, 1986.
[13] George E Uhlenbeck and Leonard S Ornstein. On the theory of the brownian motion. Physical review, 36(5):823, 1930.
[14] Emmanuel Gobet and Qihao She. Perturbation of ornstein-uhlenbeck stationary distributions: expansion and simulation. 2016.
[15] Vicky Fasen. Statistical estimation of multivariate ornstein–uhlenbeck processes and applications to co-integration. Journal of Econometrics, 172(2):325–337, 2013.
[16] Attilio Meucci. Review of statistical arbitrage, cointegration, and multivariate ornstein-uhlenbeck. 2009.
[17] Krzysztof Bartoszek, Jesualdo Fuentes-González, Venelin Mitov, Jason Pienaar, Marcin Piwczyński, Radosław Puchałka, Krzysztof Spalik, and Kjetil Lysne Voje. Model Selection Performance in Phylogenetic Comparative Methods Under Multivariate Ornstein–Uhlenbeck Models of Trait Evolution. Systematic Biology, 72(2):275–293, 12 2022.
[18] Claude Godrèche and Jean-Marc Luck. Characterising the nonequilibrium stationary states of ornstein–uhlenbeck processes. Journal of Physics A: Mathematical and Theoretical, 52(3):035002, 2018.
[19] Valentin Courgeau and Almut ED Veraart. Likelihood theory for the graph ornstein-uhlenbeck process. Statistical Inference for Stochastic Processes, pages 1–34, 2022.
[20] Alexander S. Poznyak. Chapter 9 - stable matrices and polynomials. In Alexander S. Poznyak, editor, Advanced Mathematical Tools for Automatic Control Engineers: Deterministic Techniques, pages 139–174. Elsevier, Oxford, 2008.
[21] Bernt Øksendal. Stochastic Differential Equations, pages 65–84. Springer Berlin Heidelberg, Berlin, Heidelberg, 2003.
[22] M. Gilson, N. E. Kouvaris, G. Deco, and G. Zamora-López. Framework based on communicability and flow to analyze complex network dynamics. Phys. Rev. E, 97:052301, May 2018.
[23] Ivan T Abraham. On geometric & topological methods for analysis of biophysical time series data. PhD thesis, 2022.
[24] Yongli Li, Tianchen Wang, Baiqing Sun, and Chao Liu. Detecting the lead–lag effect in stock markets: definition, patterns, and investment strategies. Financial Innovation, 8(1):51, 2022.
[25] Anish Mitra, Abraham Z Snyder, Carl D Hacker, and Marcus E Raichle. Lag structure in resting-state fmri. Journal of neurophysiology, 111(11):2374–2391, 2014.
[26] Hiroaki Sakoe and Seibi Chiba. Dynamic programming algorithm optimization for spoken word recognition. IEEE transactions on acoustics, speech, and signal processing, 26(1):43–49, 1978.
[27] Rajesh Singh, Dipanjan Ghosh, and R Adhikari. Fast bayesian inference of the multivariate ornstein-uhlenbeck process. Physical Review E, 98(1):012136, 2018.
[28] Matthieu Gilson, Enzo Tagliazucchi, and Rodrigo Cofré. Entropy production of multivariate ornstein-uhlenbeck processes correlates with consciousness levels in the human brain. Physical Review E, 107(2):024121, 2023.
[29] Hong Qian. Mathematical formalism for isothermal linear irreversibility. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 457(2011):1645–1655, 2001.
[30] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[31] Robert M Gray et al. Toeplitz and circulant matrices: A review. Foundations and Trends® in Communications and Information Theory, 2(3):155–239, 2006.
[32] Nikola Sandrić. A note on the birkhoff ergodic theorem. Results in mathematics, 72:715–730, 2017.
[33] Vincenzo Capasso and Devaid Bakstein. Introduction to Continuous-Time Stochastic Processes. Springer, 2021.
[34] Rafail Khasminskii. Stochastic stability of differential equations. Springer, 2012.
[35] Rick Durrett. Probability: theory and examples, volume 49. Cambridge university press, 2019.
[36] D Martin and LV Ahlfors. Complex analysis. New York: McGraw-Hill, 1966.
[37] L-K Hua. Introduction to number theory. Springer Science & Business Media, 1982.
[38] Ronald L. Graham, Donald E. Knuth, and Oren Patashnik. Concrete Mathematics: A Foundation for Computer Science. Addison-Wesley, Reading, 1989.
[39] Richard S Varga. Geršgorin and his circles, volume 36. Springer Science & Business Media, 2011.
[40] Peter Denton, Stephen Parke, Terence Tao, and Xining Zhang. Eigenvectors from eigenvalues: A survey of a basic identity in linear algebra. Bulletin of the American Mathematical Society, 59(1):31–58, 2022.
[41] Silvia Noschese, Lionello Pasquini, and Lothar Reichel. Tridiagonal toeplitz matrices: properties and novel applications. Numerical linear algebra with applications, 20(2):302–326, 2013.
[42] Peter E. Kloeden and Eckhard Platen. Introduction to Stochastic Time Discrete Approximation. Springer Berlin Heidelberg, Berlin, Heidelberg, 1992.

	$\displaystyle\\|\mathbf{b}(t,\mathbf{y})\\|+\\|\boldsymbol{\sigma}(t,\mathbf{y})% \\|_{F}$	$\displaystyle\leq C(1+\\|\mathbf{y}\\|)$		(2.22)
	$\displaystyle\\|\mathbf{b}(t,\mathbf{z})-\mathbf{b}(t,\mathbf{y})\\|$	$\displaystyle\leq D(\\|\mathbf{z}-\mathbf{y}\\|),$		(2.23)

	$\displaystyle\\|\mathbf{b}(\mathbf{y},t)\\|+\\|\boldsymbol{\sigma}(\mathbf{y},t)% \\|_{F}$	$\displaystyle=\\|\mathbf{b}(\mathbf{y},t)\\|+\\|\boldsymbol{\Sigma}\\|_{F}$
		$\displaystyle\leq D\\|\mathbf{y}\\|+\\|\boldsymbol{\Sigma}\\|_{F}$
		$\displaystyle\leq C\\|\mathbf{y}\\|+C=C(1+\\|\mathbf{y}\\|).$

$\displaystyle\left\|\frac{R_{n,m}(t)}{t}\right\|$	$\displaystyle\leq\left\|\frac{R_{n,m}(t)-R_{n,m}(K(t))}{t}\right\|+\left\|\frac{R% _{n,m}(K(t))}{t}\right\|$
	$\displaystyle=\left\|\frac{R_{n,m}(t)-R_{n,m}(K(t))}{t}\right\|+\left\|\frac{R_{n% ,m}(K(t))}{K(t)}\right\|$
	$\displaystyle\leq\frac{K(t)^{\frac{2}{3}}}{K(t)}+\left\|\frac{R_{n,m}(K(t))}{K(% t)}\right\|$
	$\displaystyle=\frac{1}{K(t)^{\frac{1}{3}}}+\left\|\frac{R_{n,m}(K(t))}{K(t)}% \right\|,$	(3.39)