1

A Widely Linear Complex Autoregressive Process
of Order One

arXiv:1511.04128v3 [stat.ME] 15 Mar 2017

Adam M. Sykulski, Member, IEEE, Sofia C. Olhede, Member, IEEE, and Jonathan M. Lilly, Senior Member, IEEE

Abstract—We propose a simple stochastic process for modeling
improper or noncircular complex-valued signals. The process is
a natural extension of a complex-valued autoregressive process,
extended to include a widely linear autoregressive term. This
process can then capture elliptical, as opposed to circular,
stochastic oscillations in a bivariate signal. The process is order one and is more parsimonious than alternative stochastic
modeling approaches in the literature. We provide conditions for
stationarity, and derive the form of the covariance and relation
sequence of this model. We describe how parameter estimation
can be efficiently performed both in the time and frequency
domain. We demonstrate the practical utility of the process
in capturing elliptical oscillations that are naturally present in
seismic signals.

c 2016 IEEE. Personal use of this material is permitted.
Permission from IEEE must be obtained for all other uses, in
any current or future media, including reprinting/republishing
this material for advertising or promotional purposes, creating
new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in
other works.
I. I NTRODUCTION
OMPLEX-valued stochastic processes are useful models
for parameterizing bivariate signals. Such models are in
widespread use in applications including oceanography [1]
and functional Magnetic Resonance Imaging [2]. The theory
for complex-valued representations has been developed both
in the context of stochastic processes [3], as well as for
signal processing [4], with notable recent developments in [5]
and [6]. The complex-valued representation is sometimes
preferred to the bivariate representation, due to its compactness
and interpretability [7]. For example, complex-valued signals
can be naturally decomposed into analytic and anti-analytic
signals [8], and provide a practical framework for assessing
impropriety or noncircularity in a complex signal [9]. On the
other hand, the bivariate vector representation can provide
better physical understanding of the generating mechanism of
the modeled process, and its usage is commonplace in the time
series community, see e.g. [10, Ch. 10–11].

C

The work of A. M. Sykulski was supported by a Marie Curie International Outgoing Fellowship within the 7th European Community Framework
Programme. S. C. Olhede acknowledges funding from EPSRC Responsive
mode (EP/L025744/1), EPSRC Leadership Fellowship (EP/I005250/1), and
ERC CoG 2015 - 682172 NETS. The work of J. M. Lilly was supported by
award 1031002 from the Physical Oceanography program of the United States
National Science Foundation.
A. M. Sykulski and S. C. Olhede are with the Department of Statistical
Science, University College London, Gower Street, London WC1E 6BT, UK
(emails: a.sykulski@ucl.ac.uk, s.olhede@ucl.ac.uk).
J. M. Lilly is with NorthWest Research Associates, PO Box 3027, Bellevue,
WA, USA (email: lilly@nwra.com)

The complex autoregressive process [11], [12] is a generalization of a real-valued autoregressive process, in which
the autoregressive coefficients and noise increments are both
complex-valued. Typically, the real and imaginary components
of the noise are assumed to be independent and identically
distributed [12], thus creating statistically isotropic (i.e. circular) oscillations in the signal. The process is therefore said
to be statistically circular or proper, defined subsequently, as
opposed to one that is noncircular or improper.
In many real-world observations of complex-valued signals,
noncircular or improper structure is expected to be present;
examples include seismic traces [13], oceanographic velocity
measurements [14], and wind measurements [5]. In Fig. 1, we
display a bivariate signal from a seismic trace of the 1991
Solomon Islands Earthquake, previously studied in [15]–[17].
In Fig. 1(c), the elliptical orbital shape of the oscillations
become apparent when viewed in the complex plane, and we
will refer to such motion as “elliptical oscillations.” In cases
such as these, in which the signal appears improper, a proper
process would be a poor choice of model and would fail to
summarize important characteristics of the data.
Motivated by this, in this paper we generalize the complex
autoregressive processes to a widely linear complex autoregressive process that is statistically improper or noncircular.
The notion of wide linearity was introduced in [18], and
we use this concept to relate the complex-valued process Zt
∗
to its previous value Zt−1 and its complex conjugate Zt−1
.
Our process is order one and hence Markovian, such that
Zt is not dependent on Zt−k (given Zt−1 ) for k > 1, and
is a special case of the widely linear autoregressive moving
average (ARMA) model of [19]. The widely linear ARMA is a
general and flexible framework for improper processes that is
well understood in the context of signal prediction [19], [20],
and signal estimation [21], in settings where parameter values
are assumed to be known. In practice, these parameters would
need to be estimated when modeling real-world signals.
In this paper, we provide time- and frequency-domain
techniques for estimating the parameters of our order one
widely linear process, where we will also derive the form of
the covariance and relation sequences. A key innovation will
be to relate the process to a real-valued bivariate vector autoregressive process, which will provide intuitive understanding
of the process generating mechanism, and will simplify the
problem to the estimation of five real-valued parameters.
We will demonstrate that despite its simplicity, the model
we propose can effectively capture the elliptical oscillatory
structure present in the seismic signal displayed in Fig. 1.
We contrast our widely linear process with the improper

2

10

0

5

-10
(b)

Yt × 104

(c)

15

Yt × 104

Xt × 104

(a)
10

0

10

-5

0

-10

-10

-15
16:25 16:30 16:35 16:40 16:45

-15 -10

-5

0

5

10

15

Xt × 104

time (UTC)

Fig. 1. Seismic traces from the Feb. 9th, 1991 Solomon Islands earthquake as
measured from the Pasadena recording station in California. The radial component, Xt , is displayed in (a) and the vertical component, Yt , is displayed
in (b). In (c) we display the complex-valued signal, Zt = Xt + iYt , from
16:31:55 to 16:44:35 (UTC), as indicated by the vertical dashed boundaries
in (a) and (b).

complex autoregressive process of [12], [22], where the impropriety is created using improper noise. We will demonstrate
how the inclusion of a widely linear autoregressive component
allows our process to reproduce stochastic elliptical oscillations despite being order one, whereas to generate elliptical
oscillations in the framework of [12], an order two process is
required, as investigated in [23]. We propose our order one
process as a simpler and more intuitive model for generating
elliptical oscillations in a complex-valued signal.
II. BACKGROUND
Consider a zero-mean complex-valued stationary stochastic
process Zt where t ∈ Z. The covariance sequence, sτ , and the
relation sequence, rτ , are defined at lag τ ∈ Z by
∗
sτ = E{Zt Zt+τ
},

rτ = E{Zt Zt+τ },

where E{·} denotes expectation and Zt∗ is the complex conjugate of Zt . The process Zt is said to be proper if
rτ = 0,

p
X
j=1

∀τ ∈ Z,

gj Zt−j + νt ,

gj , νt ∈ C,

Zt =

p
X

gj Zt−j +

(1)

where p is the order of the process, and where νt is a noise
process that is permitted to be noncircular or improper. The
coefficients gj are in general assumed to be complex-valued.
An important special case of this process is the proper complex
autoregressive process of order one [11], which we denote
by Zt′ and has three real-valued parameters {a, θ, σǫ2 }, and is
given by
′
Zt′ = aeiθ Zt−1
+ ǫt , a ≥ 0,
(2)

p
X

∗
hj Zt−j
+

q
X

kj ǫt−j +

q
X

lj ǫ∗t−j . (3)

j=0

j=0

j=1

j=1

This more general framework has larger flexibility in modeling improper signals. A challenge however is that as
{gj , hj , kj , lj } are complex-valued, then the number of parameters that need to be estimated is large, even for moderate
values of p and q. In this paper we shall focus on the special
case of p = 1 and q = 0 in (3), thus creating a simple orderone widely linear improper process.
Finally, [12] more generally show that any second-order
stationary process can be expressed as a widely linear filter
of a complex-proper white noise process ǫt such that
Zt =

∞
X

j=−∞

and is improper otherwise. A proper process therefore has a
relation sequence equal to zero, and is also commonly referred
to as a circular process [4]. It can be shown that the secondorder statistical properties of a proper process are isotropic,
such that its distribution is invariant to rotation. Second-order
proper processes can be treated much like second-order realvalued processes, in that the second-order properties are fully
specified by the covariance sequence.
The improper complex autoregressive process of [12], [22]
is given by
Zt =

where the complex autoregressive coefficient has been expressed in terms of an amplitude a and phase θ, both realvalued. Here {ǫt } is a sequence of i.i.d. complex-valued Gaussian noise where the real and imaginary parts are independent,
and each has zero mean and variance given by σǫ2 > 0. The
process is stationary if and only if a < 1, with variance given
by 2σǫ2 /(1 − a2 ). In such cases, a is commonly referred to as
the damping parameter. For identifiability a is not permitted to
be negative, as a negative autoregression is achieved instead
when a > 0 and π/2 + 2kπ < θ < 3π/2 + 2kπ, where
k ∈ Z. The parameter θ ∈ R is the angle of a rotation of
the process at each time step and is usually referred to as
the spin parameter [24]. The process Zt′ is an example of a
proper process, and we will refer back to this process when
we construct our widely linear improper process.
The improper complex autoregressive model (1) was generalized to a class of autoregressive moving average (ARMA)
models in [19], which uses finite-length widely linear filters
on both the autoregressive and noise components such that

kj ǫt−j +

∞
X

lj ǫ∗t−j ,

j=−∞

kj , lj , ǫt ∈ C.

(4)

This is in essence the complex-analogue to the well-known
Wold decomposition for real-valued processes [25], and in
general this will be an infinite-order process.
III. T HE W IDELY L INEAR C OMPLEX

AUTOREGRESSIVE

ORDER ONE PROCESS

In this section we introduce the widely linear complex
autoregressive process of order one. We do this by extending (2) to a widely linear form for the autoregressive and noise
components of the process, such that it is parameterized by
seven real-valued parameters. We will subsequently reduce this
model to five free parameters, by constraining two parameters,
for practical reasons that we shall discuss shortly.
We call our process the widely linear complex autoregressive process of order one, denoted by Zt , which has parameters
{λ, α, γ, φ, σν2 , cν }, and is given by
∗
Zt = λeiα Zt−1 + γeiφ Zt−1
+ νt ,

λ, γ ≥ 0,

(5)

where {νt } is a sequence of i.i.d. complex-valued Gaussian
noise with variance σν2 = E{νt νt∗ } > 0 and relation at lag
zero specified by cν = E{νt νt } ∈ C. The noise process
νt is specified by three real-valued parameters as cν is in

3

(a)

8

8

4

4
0

Y

0
-4

-4

-8

-8

-12
-12

-12
-12

-8

-4

0

X

8

12

3

2

2

1

1

0

-1

-2

-2

-3

-3
-1

0

X

0

1

2

3

4

4

8

12

t

(d)

0

-1

-2

-4

4

3

-3

-8

X

Yt

t

4

t

(c)

4

-4
-4

(b)

12

t

Yt

12

Y

general complex-valued, therefore our process has seven realvalued parameters. The remaining four parameters {λ, α, γ, φ}
define an iterative relationship between the process Zt and the
∗
previous value Zt−1 as well as its conjugate Zt−1
, where γ
and φ are respective analogues of the damping parameter, λ,
and the spin parameter, α, in this widely linear setting. Further
intuition for each of the parameters is gained when we derive
covariance sequences in Sections IV.
The widely linear complex autoregressive process of order
one is recovered from the ARMA model (3) of [19] by taking
p = 1 and q = 0. Furthermore, we combine the widely linear
noise term k0 ǫt + l0 ǫ∗t into the noise process νt , which has
variance σν2 = |k0 |2 +|l0 |2 and relation at lag zero cν = 2k0 l0 .
For this reason νt is commonly referred to as doubly white
noise [12], as it is the pointwise superposition of two complexproper white noise processes. As we have implicitly chosen
the initial phase angle of the noise, we have for parsimony
reduced the number of real-valued parameters from the eight
used in (3) (with p = 1 and q = 0), to the seven we have
defined in (5), without any loss of generality.
The model (5) introduces wide linearity in the autoregressive component of an order one process, and allows the process
to map out stochastic elliptical oscillations, as we now demonstrate. In Fig. 2, we contrast realizations from our process with
an improper order one process (1) from the framework of [12]
using doubly white noise. The model (1) only has impropriety
in the noise component—equivalent to setting γ = 0 in (5). It
is clear from the figure that when the autoregressive coefficient
λ is close to unity (in panels (a) and (b)), then the widely
linear complex autoregressive process (5) has a tendency to
generate elliptical oscillations (panel (a)), whereas an improper
order one process (1) generates oscillations that appear to be
circular (panel (b)). This is because the motion (1) is largely
determined by Zt = λeiα Zt−1 , a deterministic component
which specifies a circular oscillation. When λ is low (panels
(c) and (d)), then the noise term νt dominates both signals.
As νt is doubly white noise, then both processes are improper
in their distribution, but neither generate elliptical oscillations
that resemble the seismic traces seen in Fig 1(c)—this can only
be generated using (5) with a larger autoregressive coefficient,
and with γ > 0.
In general, the process (5) has seven real-valued parameters,
whereas (2) has three. The additional four parameters in (5)
are present because each component in (2)—the deterministic
autoregressive component and the noise component—has been
given its own improper elliptical structure with two parameters
respectively, one to stretch (γ and |cν |) and one to rotate (φ and
arg{cν }). Here we have defined |cν | as the complex modulus
of cν , and arg{cν } as the complex argument or phase. In
practical problems the four parameters {γ, φ, |cν |, arg{cν }}
will be difficult to identify simultaneously from observed
signals. For this reason we reduce the widely linear process (5)
to five free parameters, effectively making two parameters
redundant, which in our case will be |cν | and arg{cν }. This
is achieved by “aligning” the elliptical structure of the deterministic autoregression and the stochastic noise. The appeal of
reducing to five parameters is that the process will then have
the same number of free parameters as an improper order one

-4
-4

-3

-2

-1

0

X

t

1

2

3

4

t

Fig. 2. Simulated signals of length 512 (—-), from the widely linear complex
autoregressive process (5) in (a) and (c), and from an order one improper
process (1) in (b) and (d). In (a) and (b), the parameters are λ = 0.99,
α = π/6, σν2 = 1, and cν = exp(3πi/4)/5. In (c) and (d), the parameters
are λ = 0.297, α = π/6, σν2 = 1, and cν = 2 exp(3πi/4)/3. For the widely
linear complex autoregressive process (5) in (a) and (c), we additionally set
γ = 0.099 and φ = −π/4, which are set to zero in (b) and (d) for the
model (1). Expected second moment ellipses for the standard deviation are
also given (- - -), which are calculated from Section IV.

linear model (1), but will produce elliptical oscillations using
widely linear forcing (cf. Fig. 2). The order two process for
elliptical oscillations of [23] in general requires seven realvalued parameters.
To reduce our model to five free parameters, we equate
the process to a bivariate process with elliptical covariance
structure. We start from the complex-proper autoregressive
process of order one (2). We rewrite this process as a bivariate
process in terms of (Xt′ Yt′ )T , where Zt′ = Xt′ + iYt′ ,
 ′
 
 ′ 
ǫ1,t
Xt
Xt−1
,
(6)
=
aR
+
σ
ǫ
′
ǫ2,t
Yt′
Yt−1
where,
R=



cos θ − sin θ
.
sin θ cos θ

It follows that (6) is an isotropic real-valued bivariate vector
process. Here ǫ1,t and ǫ2,t are i.i.d. Gaussian random variables
with mean 0 and variance 1. The matrix R accomplishes a
rotation by the angle −π < θ ≤ π in the Cartesian plane. For
identifiability we set a ≥ 0 and for stationarity we require a <
1, as we show in Appendix A. We then construct a new process
(Xt Yt )T , which we call the elliptical bivariate autoregressive
process (of order one), using an elliptical transformation
 ′
 
Xt
Xt
= QP
,
(7)
Yt′
Yt
where



1 
0
cos ψ − sin ψ
,
Q=
, P = ρ
sin ψ cos ψ
0 ρ

in which 0 ≤ ψ < π defines the orientation of an ellipse,
and 0 < ρ ≤ 1 defines the eccentricity, 0 ≤ ε < 1, where

4

p
ε = 1 − ρ4 . We call this transformation “elliptical” as it first
stretches the X-axis by a factor of 1/ρ, and compresses the
Y -axis by a factor of ρ, before then rotating these axes through
angle ψ. This transformation leads to elliptical statistical
properties, and the process remains stationary for a < 1.
This can be seen by examining the 2 × 2 covariance matrices
associated with (6) and (7) which are respectively given by
 
 ′2

σǫ2
Xt Xt′ Yt′
=
E
I,
2
1 − a2
Xt′ Yt′ Yt′

and

Xt2 Xt Yt
E
Xt Yt Yt2




=



σǫ2
1 − a2



QP 2 QT ,

where I is the identity matrix. The covariance matrix of
(Xt′ Yt′ )T is isotropic or circular, whereas in general the
covariance matrix of (Xt Yt )T is elliptical with orientation
ψ and ratio of semi-minor to semi-major axis ρ2 .
The elliptical bivariate autoregressive process of order one is
defined by five parameters, namely {a, θ, ρ, ψ, σǫ2 }. We now
express the widely linear complex autoregressive process of
order one as Zt = Xt + iYt , and in Proposition 1 we relate
the parameters of this process to those of the elliptical bivariate
autoregressive process to form a five-parameter process.
Proposition 1: Suppose the process (Xt Yt )T is an elliptical
bivariate autoregressive process of order one, as defined in (6)
and (7) by the parameters {a, θ, ρ, ψ, σǫ2 }. This process is
equivalent to a widely linear complex autoregressive process
of order one (5), specified by Zt = Xt + iYt where



i sin θ 1
2
Zt−1
+ρ
Zt = a cos θ +
2
ρ2


π
a sin θ 1
∗
+
− ρ2 ei(2ψ− 2 ) Zt−1
2
ρ2
 iψ

π
e
+ σǫ
ǫ1,t + ρei(ψ+ 2 ) ǫ2,t .
ρ

The relationship with the parameters {λ, α, γ, φ, σν2 } in the
specification of (5) is given in Table I, with the final redundant
parameter, the relation at lag zero cν , specified by


 γ 
π
1
2
i2ψ
2
2
cν = σǫ
ei(φ+ 2 ) .
e
=
σ
−
ρ
(8)
ν
ρ2
λ sin α
The proof of Proposition 1 is given in Appendix B. Therefore in order to represent the elliptical bivariate process (7)
in terms of the widely linear complex process (5), we have
required to use only five parameters. We can therefore simply
fix cν using either {a, θ, ρ, ψ, σǫ2 } or {λ, α, γ, φ, σν2 }, as given
in (8). This equivalence can be verified using the transformations in Table I. We note that this is how the value of cν was
set earlier in Fig 2.
The eccentricity parameter, ε, can also be found in terms of
the widely linear complex process parameters using Table I
s
p
2γ
.
(9)
ε = 1 − ρ4 =
λ| sin α| + γ

Therefore for the widely linear process to return valid values
for the eccentricity, 0 ≤ ε < 1, we require
γ ≤ λ| sin α|.

(10)

Additionally, p
for the widely linear process to be stationary we
require a = λ2 − γ 2 < 1. Therefore from (10) it follows
that a > 0. Then we see that stationarity is guaranteed when
λ < 1, and otherwise for stationarity we require that
γ>

p
λ2 − 1,

when λ ≥ 1.

(11)

These inequalities also ensure that θ and σǫ2 return valid values
when mapping parameters from the complex to bivariate
specifications. Increasing γ increases the eccentricity, until
eventually larger values of γ are not valid. This means that
there is a range of values for γ, which depends on both λ and
α, for our five-parameter process to be valid and stationary.
Combining the inequalities in (10) and (11) we see that we
require λ2 < 1/ cos2 α for our process to be valid and
stationary. When α = 0, the case of no spin, then γ = 0
from (10) and we require λ < 1 for stationarity. However as
α increases then interestingly, there are parameter values for
which λ > 1 and the process can still be stationary, unlike the
proper case, although this then requires a non-zero γ as can
be seen from (11).
We gain further insight by analyzing the relationships in
Table I. The first three parameters on each side of the table,
{λ, α, γ} and {a, θ, ρ}, have a direct one-to-one mapping,
where {λ, α} and {a, θ} become identical as ρ → 1 or
γ → 0. The parameters a and λ are monotonic functions
of each other, as are α and θ, which have the same sign
in the range (−π, π]. The bivariate ellipse orientation ψ and
the complex-conjugate spin parameter φ are directly related,
but are adjusted depending on the sign of θ and α. The ratio
of the variance parameters, σν2 and σǫ2 , is determined by the
eccentricity. The effect of each parameter can therefore be
closely related row-by-row in Table I.
By connecting to a bivariate process, we have gained the
advantages of both specifications: we benefit from the compactness and applicability of a complex representation, and we
benefit from the interpretability and physical understanding
gained from a bivariate representation. A particularly useful
feature of complex signals is that we can perform hypothesis
tests for impropriety [9], and we will demonstrate the insight
gained from such an analysis in our seismic data example in
Section VI.

IV. C OVARIANCE

AND

R ELATION SEQUENCE

In this section we compute the covariance and relation
sequences for the widely linear complex autoregressive process
of order one. It follows directly from (5) that the process
is Gaussian, as it is a linear combination of complex-valued
Gaussian random variables. Therefore the covariance and
relations sequences fully specify the process. These sequences
would have complicated expansions if expressed analytically,
so instead we find recurrence relationships between the lags.
2
) and relation at lag zero
First we find that the variance (σZ

5

TABLE I
T HIS TABLE PROVIDES A MAPPING BETWEEN THE PARAMETERS OF THE ELLIPTICAL BIVARIATE PROCESS OF ORDER ONE (7), AND THE WIDELY LINEAR
2 > 0 FOR THE WIDELY LINEAR COMPLEX
COMPLEX AUTOREGRESSIVE PROCESSES OF ORDER ONE (5). W E REQUIRE λ ≥ 0, −π < α ≤ π, γ ≥ 0 AND σν
PROCESS , AND a ≥ 0, −π < θ ≤ π, 0 < ρ ≤ 1 AND σǫ2 > 0 FOR THE ELLIPTICAL BIVARIATE PROCESS . T HE PARAMETERS φ AND ψ ARE UNRESTRICTED
IN THIS MAPPING . T HE FUNCTION atan2 IS THE FOUR QUADRANT INVERSE TANGENT, acos IS THE INVERSE COSINE FUNCTION AND sgn IS THE SIGNUM
FUNCTION . T HESE FUNCTIONS ARE CHOSEN TO ENSURE THAT α AND θ HAVE A ONE - TO - ONE MAPPING IN THE RANGE (−π, π].

Bivariate elliptical to widely linear complex
r

2
2
λ = a cos2 θ + sin4 θ ρ12 + ρ2
α = atan2
γ=

n

a
| sin θ|
2

sin θ
2





− ρ2

1
ρ2

1
ρ2


o
+ ρ2 , cos θ

θ = sgn(α)acos



ρ=

φ = 2ψ − sgn(θ) π2


σν2 = σǫ2 ρ12 + ρ2

ψ=

2
σZ
= E{Zt Zt∗ }


∗
= E λeiα Zt−1 + γeiφ Zt−1
+ νt ×

∗
λe−iα Zt−1
+ γe−iφ Zt−1 + νt∗



2
= (λ2 + γ 2 )σZ
+ λγei(α−φ) cZ + λγei(φ−α) c∗Z + σν2 ,
(12)

cZ = E{Zt Zt }


∗
= E λeiα Zt−1 + γeiφ Zt−1
+ νt ×

∗
λeiα Zt−1 + γeiφ Zt−1
+ νt

=

2 2iα

+λ e

cZ +

γ 2 e2iφ c∗Z

+ cν .



(13)

2
σZ

We now combine (12) and (13) to solve for
and cZ ,
yielding
 2 
 2 
λ2 + γ 2
λγei(α−φ) λγei(φ−α)
σZ
σZ
 cZ  =  2λγei(α+φ)
λ2 e2iα
γ 2 e2iφ   cZ 
c∗Z
c∗Z
2λγe−i(α+φ) γ 2 e−2iφ λ2 e−2iα
 2
σν
+  cν  (14)
c∗ν

and hence

where



λ| sin α|−γ
λ| sin α|+γ

φ
2

 2
 2
σZ
σν
 cZ  = M −1  cν  ,
c∗Z
c∗ν


1 − λ2 − γ 2 −λγei(α−φ) −λγei(φ−α)
M =  −2λγei(α+φ) 1 − λ2 e2iα
−γ 2 e2iφ  .
−i(α+φ)
2 −2iφ
−2λγe
−γ e
1 − λ2 e−2iα

(15)



The analytic form for the inverse matrix M −1 in (15) is
provided as part of the online software available at http://
ucl.ac.uk/statistics/research/spg/software
and is not included here for space considerations. When
simulating signals from the process, in addition to satisfying
the inequalities specified in (10), then for the process to be
stationary we require that the first observation is generated

q

λ2
λ2 −γ 2


cos α

1/4

+ sgn(α) π4
√

σǫ2 = σν2

(cZ ) are given by

2
2λγei(α+φ) σZ

Widely linear complex to bivariate elliptical
p
a = λ2 − γ 2

λ2 sin2 α−γ 2
2λ| sin α|

from the complex-valued normal distribution with mean zero,
2
variance σZ
, and relation at lag zero cZ . See Section V-A for
more detail on the complex-valued normal distribution. This
is how the signals in Fig. 2 have been simulated, and more
details on this can be found in the supporting online code.
2
After computing σZ
and cZ from (15), the covariance
sequence, sτ , and the relation sequence, rτ , can be found for
any τ > 0 using the following recurrence relationships

∗
sτ = E Zt Zt+τ


∗
−iφ
∗
Zt+τ −1 + νt+τ
= E Zt λe−iφ Zt+τ
−1 + γe
−1
= λe−iα sτ −1 + γe−iφ rτ −1 ,
rτ = E {Zt Zt+τ }


∗
= E Zt λeiφ Zt+τ −1 + γeiφ Zt+τ
−1 + νt+τ −1
= λeiα rτ −1 + γeiφ sτ −1 .

(16)

(17)

2
Therefore, after solving for s0 = σZ
and r0 = cZ , we iterate
to find sτ and rτ using (16) and (17). To find the sequences
for negative lags we use the simple relationship s−τ = s∗τ and
r−τ = rτ .
From (16) and (17) we see that λ and α contribute to the
exponential decay of the autocovariance—this is expected as
autoregressive processes are short memory. Conversely, γ and
φ have a “flipping” effect on sτ and rτ , where the covariance
is dependent on the relation at previous lags and vice-versa.
This is a consequence of the widely linear form in (5), where
the conjugate action creates an iterative ‘flip’ of the process
about the real axis in the complex plane.

V. M ETHODS FOR E STIMATING PARAMETERS
In this section we detail how the parameters of the widely
linear complex autoregressive process of order one can be
estimated from an observed signal using maximum likelihood. We first resolve the exact form of the likelihood in
Section V-A, and then provide an approximate method in
the frequency domain in Section V-B. The latter method has
practical advantages in real-world applications, as we shall
discuss.

6

A. Maximum Likelihood
Suppose that Zt follows a stationary widely linear complex
autoregressive process or order one, as specified in (5) with
parameters satisfying (10) and (11), then the probability distribution of Zt follows a complex-valued normal distribution.
For a general complex-valued normally distributed random
variable z, we denote its distribution by NC (µ, σ 2 , c), with
mean µ, variance σ 2 , and relation at lag zero c. The probability
density function of z is then given by
1
×
p(z) = p
π (σ 2 )2 + |c|2
(

 
)
 σ 2 c −1 z − µ
1 ∗
∗
exp − z − µ z − µ
.
c∗ σ 2
z ∗ − µ∗
2

It then follows that the probability distribution of the widely
linear complex autoregressive process or order one, denoted
Zt , is given by

2
, cZ ,
Zt ∼ NC 0, σZ

2
where σZ
and cZ are found using (15). For a given observed
signal z0 , . . . , zN −1 from the process Zt , the probability of observing the first value z0 directly follows from the probability
distribution of Zt ,

1
p(z0 ; θ) = p 2
×
π (σZ )2 + |cZ |2
(
−1  )
 2
 σZ
1 ∗
z0
cZ
,
exp − z0 z0
2
z0∗
c∗Z σZ
2

2
where θ = {σZ
, cZ }. Next we make use of the Markovian
property of the process to find the conditional distribution of
Zt given Zt−1 = zt−1 , for 1 ≤ t ≤ N − 1

∗
, σν2 , cν ,
(Zt |Zt−1 = zt−1 ) ∼ NC λeiα zt−1 + γeiφ zt−1

such that the condition probability of observing zt given zt−1
is
1
p(zt |zt−1 ; θ) = p
×
π (σν2 )2 + |cν |2
(
)
 

 σν2 cν −1 zt − µzt
1 ∗
∗
,
exp − zt − µzt zt − µzt
c∗ν σν2
zt∗ − µ∗zt
2

where θ = {µzt , σν2 , cν } and

∗
.
µzt = λeiα zt−1 + γeiφ zt−1

(18)

The likelihood of observing the signal z0 , . . . , zN −1 is found
by evaluating
p(z0 , . . . , zN −1 ; θ) = p(z0 ; θ)

N
−1
Y
t=1

p(zt |zt−1 ; θ).

The log-likelihood (denoted ℓt (θ)) is therefore
ℓt (θ) = log (p(z0 ; θ)) +

N
−1
X
t=1

log (p(zt |zt−1 ; θ)) ,

which for a widely linear complex autoregressive process of
order one is found to be

N −1
log (σν2 )2 + |cν |2
ℓt (θ) = −N log π −
2


 
N
−1
X
 σν2 cν −1 zt − µzt
1 ∗
zt − µ∗zt zt − µzt
−
c∗ν σν2
zt∗ − µ∗zt
2
t=1
 2
−1  
 σZ
 1
1
cZ
z0
2 2
− log (σZ
) + |cZ |2 − z0∗ z0
,
2
c∗Z σZ
z0∗
2
2
(19)
with µzt given in (18). The optimal parameter choice θ̂ is then
found by maximizing the log-likelihood (19)
θ̂ = arg max ℓt (θ),
θ∈Θ

(20)

where Θ is the permitted parameter range for θ, recalling that
for the five-parameter process, cν is specified by (8), and that
the inequalities (10) and (11) should also be satisfied.
B. Frequency Domain “Whittle” Likelihood
The parameters of the widely linear complex autoregressive
process of order one can also be computed in the frequency
domain using Whittle’s approximation to maximum likelihood
[26], known as the ‘Whittle likelihood.’ This approximation
of the time-domain likelihood is in the frequency domain,
and relies solely on applying Fourier Transforms which can
be computed in O(N log N ) operations. We use a biascorrected form of the Whittle likelihood, which was extended
to complex-valued signals in [27] and is given by
X
ℓ(θ) = −
log (|f (ω; θ)|) + J H (ω)f −1 (ω; θ)J(ω) ,
ω∈Ω

(21)
where Ω is the set of Fourier frequencies used in the estimation, θ is the unknown parameter vector, and J(ω) and f (ω; θ)
are given by

N −1 
1 X Zt
−iωt
,
(22)
J(ω) = √
∗ e
N t=0 Zt


S̄(ω; θ) R̄(ω; θ)
f (ω; θ) =
.
(23)
R̄∗ (ω; θ) S̄(−ω; θ)

The vector J(ω) is the Discrete Fourier Transform for a
signal Z0 , . . . , ZN −1 , evaluated at the Fourier frequencies. The
matrix f (ω; θ) contains the expected periodogram, S̄(ω; θ) =
E{|J(ω)|2 }, and the expected complementary periodogram,
R̄(ω; θ) = E{J(ω)J(ω)}, and is dependent on both the signal
length N , and the parameter vector θ. Expected periodograms
are used in (21), rather than theoretical spectral densities,
as this removes the known bias effects from using Discrete
Fourier Transforms in (21) for finite N . We then compute
S̄(ω; θ) and R̄(ω; θ) directly from the covariance and relation
sequences using the relationships


N
−1
X
|τ | −iωτ
sτ (θ) 1 −
e
,
(24)
S̄(ω; θ) =
N
τ =−(N −1)

R̄(ω; θ) =

N
−1
X

τ =−(N −1)



|τ | −iωτ
rτ (θ) 1 −
e
.
N

(25)

(a)

(b)
15
10

Yt × 104

130
120
110
100
90
80
70
60
50
40
-π

5
0
-5
-10
-15

-π/2 -π/4 0

π/4 π/2

π

-20

130

0
t

(c)

130
120

110

110

power (dB)

120

100
90
80
70
60
-π

-10

10

20

π/4 π/2

π

X × 104

frequency (ω, rad/s)

power (dB)

The expected complementary periodogram R̄(ω; θ) will in
general be a complex-valued quantity. The usefulness of (24)
and (25) is that they can be computed in O(N log N ) time
as they are Discrete Fourier Transforms, ensuring the Whittle
likelihood remains an O(N log N ) procedure in this bias
corrected version. This method is proven to be statistically
consistent in [28], and is shown to significantly remove bias
effects for sample sizes as large as 1,000 data points.
For the widely linear complex autoregressive process of
order one, sτ and rτ in (24) and (25) are computed using
the recurrence relationships given in Section IV, which is an
O(N ) computation. Both f (ω; θ) and J(ω) are then computed
using Fast Fourier Transforms, thus giving O(N log N ) efficiency overall. The optimal parameter choice θ̂ is then found
by maximizing the Whittle likelihood in the same way as (20).
The advantage of performing maximum likelihood in the
frequency domain is that we can restrict the range of frequencies used in the summation in (21). This allows the parameter
estimation to be performed semi-parametrically, by ignoring
frequencies that are known to be contaminated or not specified
well by the model. For example, this was used in [29] to
remove the effect of eddies when estimating the parameters
of a turbulent flow model for the ocean surface. We will also
employ such semi-parametric procedures in our data example
in Section VI.
We note that the Whittle likelihood can be alternatively used
with tapered spectral estimates in (22), where the triangle
kernel 1 − |τ |/N in (24) and (25) is then replaced with a
modified kernel of smaller width, as documented in [28].
Tapering the likelihood helps reduce mean square error in
parameter estimates when the process is long memory or has
steep spectral slopes. As the widely linear complex autoregressive process of order one is short memory, then tapering
is unnecessary and we use the periodogram approach defined
in (21)–(25).

power (dB)

7

(d)

100
90
80
70

-π/2 -π/4 0

π/4 π/2

frequency (ω, rad/s)

π

60
-π

-π/2 -π/4 0

frequency (ω, rad/s)

Fig. 3. Panel (a) displays the periodogram (—) and model fit (- - -) of the
signal Zt displayed in Fig. 1(c). Panel (b) shows Zt in the interval (UTC)
16:33:46 to 16:36:26 (—) and 16:38:19 to 16:40:59 (—). The periodogram
(—) and model fit (- - -) are displayed for Zt in the intervals (c) 16:33:46
to 16:36:26 and (d) 16:38:19 to 16:40:59. In (a), (c), and (d), parameter
estimation is performed for ω ∈ [−π/4, π/4], as indicated by the vertical
dashed boundaries, and we extend the fitted lines to all frequencies (-·-).

complex plane in Fig. 1(c). The signal has evident improper
structure, as can be seen from the elliptical paths of the signal
in the complex plane.
We first fit the widely linear complex autoregressive process
of order one to the entire signal displayed in Fig. 1(c) using the
Whittle likelihood, as detailed in Section V. The periodogram
of the data, and the resulting model fit of the periodogram,
are displayed in Fig. 3(a). For all parameter estimates in
this section, we perform the Whittle likelihood estimation
semi-parametrically over a reduced range of frequencies, ω ∈
[−π/4, π/4], as the signal energy is strongly concentrated
within this frequency range. For complex-valued signals, the
spectrum is defined at both negative and positive frequencies,
and will in general be asymmetric. The two peaks of different
VI. A PPLICATION TO S EISMIC DATA
magnitude on each side of the spectra, at approximately the
In this section we investigate using the widely linear same frequency, indicate elliptical oscillatory motion. Our
complex autoregressive process of order one as a model for fitted process has located the frequency of these peaks, but
seismic data. We analyze the seismic trace from the Feb 9th overall is a poor fit to the periodogram. This is due to the non1991 Solomon Islands earthquake, as presented in Fig. 1 in stationarity of the signal. Inspecting Fig. 1(c) in more detail we
Section I. The data is sampled every second and is freely can see that the amplitude, eccentricity and orientation of the
available from http://ds.iris.edu/wilber3. All the elliptical oscillations are changing in time. Our model, which
results in this section (and all figures in this paper) are exactly is stationary, is not able to capture this variable structure.
reproducible with MATLAB code available from http://
To investigate these nonstationary effects we separately
ucl.ac.uk/statistics/research/spg/software. analyze two segments of the data, each 161 seconds (or data
The seismic trace consists of three components: a radial, points) long, spanning the periods 16:33:46 to 16:36:26 and
vertical, and transverse signal [30]. We model the radial and 16:38:19 to 16:40:59 (UTC) respectively. The complex-valued
vertical components as a complex-valued signal, as they are signals corresponding to these time periods are displayed in
strongly coupled due to the presence of a Rayleigh wave. We Fig. 3(b). The choice of window length is motivated by the
do not analyze the transverse component as the expression of example signal itself, and has been selected such that it is
the Rayleigh wave does not exist in the transverse signal [17]. as short as possible (to capture as much time variability as
The radial, Xt , and vertical, Yt , signal are displayed in possible), while still being able to robustly estimate all five
Figs. 1(a) and 1(b) respectively. Our analysis will first focus free parameters. For automated windowing procedures we
on the segment between the dashed vertical lines in the figure. refer the reader to [31], which is outside the scope of this
We combine these signals within this partition to form a single paper. The periodograms and model fits of each segment are
complex-valued signal, Zt = Xt + iYt , as displayed on the displayed in Figs. 3(c) and 3(d). The process is now seen to

8

(a)

0.9

(b)
5π/8

orientation (ψ)

0.8
0.7
0.6
0.5

3π/8
π/4

0.4
0.3

π/2

16:34

16:37

16:40

time (UTC)

16:43

16:34

16:37

16:40

16:43

time (UTC)

Fig. 5. Estimates of (a) the eccentricity, ε, and (b) the orientation, ψ, in
radians, of the seismic signal of Fig. 1. (—-) are the estimates from fitting a
widely linear complex autoregressive process of order one across a 161-second
sliding window, with the 95% confidence intervals given in gray. (- · -) are the
estimates from the method of [33], and (- - -) are the estimates obtained from
the Fourier Transform evaluated at the two frequency peaks. All estimates
and confidence intervals, including the nonparametric techniques, have been
smoothed in time using a moving average window of width 11.

Fig. 4. Spectrograms using a 161-second sliding window of the seismic
signal of Fig. 1. Panel (a) is the evolving periodogram of Zt , (b) is the
evolving model fit to the periodogram using our widely linear process, (c) is
the magnitude of the evolving complementary periodogram of Zt and (d) is
the magnitude of the evolving model fit to the complementary periodogram.
The color scale is given in decibels.

be a good fit to these shorter signals. The optimal parameters
are significantly different for each segment. For example, the
eccentricity estimate in the first segment is 0.29 whereas in
the second segment it is 0.56. These differing eccentricity
estimates are related to the different ratios of the amplitudes
of the two peaks in each respective periodogram.
The appropriateness of our process to modeling shorter segments of this signal suggests that a locally stationary modeling
assumption should be used, see e.g. [32]. To investigate this in
more detail we perform the model fit to a rolling 161-second
window over the entire signal.
In Fig. 4(a) and 4(b) we plot the spectrogram of the data—
that is, a moving window of the periodogram—together with
the spectrogram of the expected periodogram from our model
fit. These spectrograms are only shown for the frequencies that
have been used in the fit. Note that the zero frequency is not
included in the fit as we have removed the sample mean for
each segment. From the figure it can be seen that the widely
linear complex autoregressive process of order one captures
the overall shape of the spectrum at each time slice, as well as
its evolution over time. In particular, the process has captured
the gradually changing frequency of the oscillations.
We also display, in Fig. 4(c) and 4(d), the time-frequency
plots for the magnitude of the complementary periodogram,
and the resulting model fit (respectively). The complementary
periodogram forms a Fourier pair with the sample relation
sequence, and as a consequence complementary periodograms
from observed improper processes are expected to exhibit
noticeable structure. This structure has been captured well in
the model fit, which is important, as the complementary periodogram is where information regarding impropriety—such
as the expected orientation of elliptical motion—is contained.
The complementary periodogram is complex-valued, but we
have not included plots for its phase here, for space saving

considerations.
In Fig. 5(a) and 5(b) we display the time-varying eccentricity and orientation parameters from the model fit, calculated using Table I and (9). Other informative time-varying
summaries can also be found, such as of the noise variance
or damping terms. These can be generated as part of the
online software. In Fig. 5 we compare with two alternative
methods. First we compare against results obtained from the
nonparametric deterministic approach of [33], which models
the signal as a time varying ellipse, thus providing a good
comparison to our results despite being a complementary approach. Secondly, we compare against a simple nonparametric
approach of comparing the Fourier transform at the positive
and negative frequency peaks in the power spectral density,
which we denote as ±ωmax , where it can be readily shown that
eccentricity and orientation estimates can be obtained from
p
2 |JZ (ωmax )JZ (−ωmax )|
ε̂ =
|JZ (ωmax )| + |JZ (−ωmax )|
1
ψ̂ = [arg{JZ (ωmax )} + arg{JZ (−ωmax )}] ,
2
with JZ (ω) defined as the top row of (22). In Fig. 5,
the estimated eccentricities and orientations broadly agree
across the different methods. The values obtained from our
parametric model are in general smoother, as the method
smooths over information across frequencies when estimating
parameters. The usefulness of our stochastic process is that it
prescribes a generating mechanism providing physical insight
and the ability to replicate signals, which the alternative purely
diagnostic metrics do not provide.
Another useful feature of a stochastic modeling approach is
that we can calculate confidence intervals for parameter estimates, by numerically computing the Hessian of the Whittle
likelihood, as detailed in [29]. In Fig. 5 we include the 95%
confidence intervals for our parameter estimates, where care
must be taken when assessing significance across time as these
are not pointwise simultaneous confidence intervals.
Finally, another advantage of the stochastic modeling approach is that we can perform a parametric hypothesis test
for impropriety, to test for statistical significance for when

9

200

wave. Propriety is rejected within our main analysis window of
Figs. 4–6, and also afterwards where there is still some seismic
activity (as can be seen in Fig. 1). We can see that the rejection
of propriety in favor of our model is most significant at time
points where the signal is most eccentric, around 16:40 (cf.
Fig. 5(a)), which is intuitive as here a circular/proper model
is the least appropriate.

W

150

100

50

0

16:26

16:30

16:34

16:38

16:42

16:46

16:50

16:54

time (UTC)

Fig. 6. The likelihood ratio statistic W (—-) over time, smoothed with a
moving average window of width 11, from fitting the proper (2) and widely
linear improper (5) complex autoregressive order one processes to Zt across a
161-second sliding window. The signal is then divided into 11 non-overlapping
windows, as indicated by the vertical gray-dashed lines, where the vertical
black solid lines indicate the analysis window of Figs. 4–6. The likelihood
ratio statistic for each of these windows is then given by (—-). We also
display the 95th percentile of a χ22 distribution (- - -), and after applying a
False Discover Rate (FDR) procedure to control the rate of false positives,
we reject propriety in each segment except the first two.

an improper model should be used. This approach is simpler
than performing the test of [9], under the assumption that the
parametric model is appropriate. We perform the parametric
test by also fitting the proper complex autoregressive order
one process (2) to rolling windows of observations, in exactly
the same manner as performed for the widely linear improper
process. We then perform a likelihood ratio test, as proposed in
[27], to see if there is significant statistical evidence to suggest
the null hypothesis of a proper process should be rejected
in favor of an improper process. To do this we compute the
likelihood ratio statistic, given by W = 2{ℓ(θalt) − ℓ(θnull )},
where “alt” and “null” denote the alternative and null models
respectively. We compare the likelihood ratio statistic with the
95th percentile of a χ22 distribution. The χ22 distribution is used
because there are two additional parameters in the alternate
than in the null.
The results of the test are displayed in Fig. 6, where we
have extended the analysis and computed W over a longer
period of time. Similarly to Fig. 5, care must be taken here
when performing such an analysis over time, as a correction
must be made for multiple testing, to control the rate of
false positives. As a result, we have divided the analysis
into 11 non-overlapping windows, as indicated, and reported
the likelihood ratio test statistic within each window. Then
to control the rate of false positives, rather than rejecting
all segments with p-values less than .05 (found using the
χ22 distribution), a False Discovery Rate (FDR) procedure is
applied using the Benjamini-Hochberg procedure [34]. This
procedure ranks the p-values in ascending order (denoted
p1 , . . . , p11 ) and finds the largest j such that pj ≤ .05j/11,
and then rejects all segments corresponding to p1 , . . . , pj . This
procedure formally requires data segments to be independent,
and while mild correlations do exist, these can only result in
positive correlations between the statistics (as χ22 distributions
can only be positively correlated). Therefore we may still
employ this procedure, but the rejection rate is conservative.
The FDR analysis, which is included in the online code,
rejects all but the first two segments, which have the highest
associated p-values and are before the arrival of the Rayleigh

VII. C ONCLUSIONS
In this paper we have proposed a widely linear complex
autoregressive process of order one. The key novelty of the
stochastic process is that impropriety is constructed by relating
the process to its conjugate at the previous timestep using
a widely linear representation, building on ideas developed
in [19] for higher order ARMA processes. Our approach is in
contrast to alternative approaches to modeling improper complex autoregressive processes, where only the noise component
is improper. Our stochastic process can generate improper
structure in the form of elliptical oscillations, which is not
possible using alternative order one processes in the literature.
We reduced the specification of our process from seven free
parameters to five free parameters, by relating the process to
a bivariate elliptical process with interchangeable parameters,
and “aligning” the ellipticity of the signal and noise. Reducing
to five free parameters has the advantage that parameters of
our process are easier to identify and estimate in real-world
problems. Furthermore, linking to a bivariate process provided
the benefits of using both representations, where we were
then able to find conditions for stationarity, and describe the
structure of the elliptical oscillations. In general the parameter
connections between the representations are non-trivial, but
transforming between representations provided useful insight,
shedding light on the behavior of the more compact complex
widely linear representation.
A promising direction for gaining insight into the full
unconstrained seven-parameter widely linear specification, is
to relate the process to a bivariate process with two separate
elliptical transforms—one each for the autoregressive and
noise components in (6)—thus also now being specified by
seven parameters. Such an analysis should be performed if the
problem is known to have elliptical signal and noise structure
that is not aligned.
Another important innovation of this paper is that we have
provided time- and frequency-domain techniques to parameter estimation, and then applied them to demonstrate how
our widely linear improper process can effectively capture
elliptical oscillations in observed seismic traces. An advantage
of the complex-valued approach is that frequency-domain
understanding then becomes natural, as the asymmetry in
the power spectra defines preferred direction of rotation.
The process we propose has the potential to be applied to
other improper signals in numerous applications, including for
example modeling elliptical eddies as documented in [14], or
modeling phase-shifted stochastic cycles in econometric time
series [35].

10

A PPENDIX A
S TATIONARITY OF

THE BIVARIATE ELLIPTICAL
AUTOREGRESSIVE PROCESS

From [36, Ch. 11] we have that for (6) to be stationary we
require that the eigenvalues of the matrix


a cos θ −a sin θ
aR =
a sin θ a cos θ
have modulus less than one. The eigenvalues of this matrix
are found to be
λ2 = a cos θ − ia sin θ.

λ1 = a cos θ + ia sin θ,

Therefore as |λ1 | = |λ2 | = a, it follows that requiring both
|λ1 | < 1 and |λ2 | < 1 for stationarity is equivalent to requiring
that a < 1 (as we already have that a ≥ 0).
P ROOF

A PPENDIX B
OF P ROPOSITION 1

Combining (6) and (7) we have the relationship
 
  ′ 
 
ǫ
Xt−1
Xt
,
= QP aR
+ σǫ 1,t
′
ǫ2,t
Yt−1
Yt
′
′
Yt−1
)T for (Xt−1 Yt−1 )T , and
and then substituting (Xt−1
−1
T
using that Q = Q , it follows that
 



 
ǫ
Xt−1
Xt
. (26)
+ σǫ 1,t
= QP aRP −1 QT
ǫ2,t
Yt−1
Yt

To simplify (26) we first define




1 0
0 −1
,
I=
, J=
1 0
0 1

K=




0 1
.
1 0

It then follows that
P RP

−1

sin θ
= cos θI+
2






sin θ 1
1
2
2
+ ρ J−
+ ρ K.
ρ2
2
ρ2

We then use the properties that
QIQT = I,

QJQT = J,

QKQT =



− sin 2ψ cos 2ψ
,
cos 2ψ sin 2ψ

which allows for a simple expression for L = QP RP −1 QT
where


sin θ 1
2
J
+
ρ
L = cos θI +
2
ρ2



sin θ 1
− sin 2ψ cos 2ψ
2
.
−
ρ
−
cos 2ψ sin 2ψ
2
ρ2
Therefore (26) simplifies to
 
 


ǫ
Xt
Xt−1
+ QP σǫ 1,t .
= aL
ǫ2,t
Yt
Yt−1

(27)

To reformulate this in terms of Zt = Xt + iYt we use the
relationship
 

 

1
Xt
Zt
1 1
= T
.
(28)
,
where
T
=
Yt
−i i
Zt∗
2

Substituting (28) into (27) and using that T H T = 2I,
where the subscript H denotes the Hermitian transpose, after
rearranging we see that
 
 


a H
Zt
ǫ
Zt−1
H
= T LT
(29)
+ T QP σǫ 1,t .
∗
Zt∗
ǫ2,t
Zt−1
2
Expanding (29) and taking the top row we then get the
relationship given in the proposition.
R EFERENCES
[1] J. Gonella, “A rotary-component method for analysing meteorological
and oceanographic vector time series,” Deep-Sea Res., vol. 19, pp. 833–
846, 1972.
[2] D. B. Rowe, “Modeling both the magnitude and phase of complexvalued fMRI data,” NeuroImage, vol. 25, no. 4, pp. 1310–1324, 2005.
[3] K. S. Miller, Complex stochastic processes: an introduction to theory
and application. Addison-Wesley, 1974.
[4] P. J. Schreier and L. L. Scharf, Statistical signal processing of complexvalued data: the theory of improper and noncircular signals. Cambridge
University Press, 2010.
[5] T. Adali, P. J. Schreier, and L. L. Scharf, “Complex-valued signal
processing: The proper way to deal with impropriety,” IEEE T. Signal
Proces., vol. 59, no. 11, pp. 5101–5125, 2011.
[6] A. T. Walden, “Rotary components, random ellipses and polarization: a
statistical perspective,” Phil. Trans. R. Soc. A, vol. 371:20110554, 2013.
[7] S. C. Olhede and A. T. Walden, “Local directional denoising,” IEEE T.
Signal Proces., vol. 53, pp. 4725–4730, 2005.
[8] ——, “‘Analytic’ wavelet thresholding,” Biometrika, vol. 91, no. 4, pp.
955–973, 2004.
[9] P. J. Schreier, L. L. Scharf, and A. Hanssen, “A generalized likelihood
ratio test for impropriety of complex signals,” IEEE Signal Proc. Let.,
vol. 13, no. 7, pp. 433–436, 2006.
[10] J. D. Hamilton, Time series analysis. Princeton University Press, 1994.
[11] A. Le Breton, “A note on maximum likelihood estimation for the
complex-valued first-order autoregressive process,” Stat. Probabil. Lett.,
vol. 7, no. 2, pp. 171–173, 1988.
[12] B. Picinbono and P. Bondon, “Second-order statistics of complex signals,” IEEE T. Signal Proces., vol. 45, no. 2, pp. 411–420, 1997.
[13] J. C. Samson, “Pure states, polarized waves, and principal components
in the spectra of multiple, geophysical time-series,” Geophys. J. Int.,
vol. 72, no. 3, pp. 647–664, 1983.
[14] J. M. Lilly and J. C. Gascard, “Wavelet ridge diagnosis of time-varying
elliptical signals with application to an oceanic eddy,” Nonlinear Proc.
Geoph., vol. 13, no. 5, pp. 467–483, 2006.
[15] J. M. Lilly and J. Park, “Multiwavelet spectral and polarization analyses
of seismic records,” Geophys. J. Int., vol. 122, no. 3, pp. 1001–1021,
1995.
[16] S. C. Olhede and A. T. Walden, “Polarization phase relationships via
multiple Morse wavelets. II. Data analysis,” P. Roy. Soc. A-Math. Phy.,
vol. 459, no. 2031, pp. 641–657, 2003.
[17] J. M. Lilly, “Modulated oscillations in three dimensions,” IEEE T. Signal
Proces., vol. 59, no. 12, pp. 5930–5943, 2011.
[18] B. Picinbono and P. Chevalier, “Widely linear estimation with complex
data,” IEEE T. Signal Proces., vol. 43, no. 8, pp. 2030–2033, 1995.
[19] J. Navarro-Moreno, “ARMA prediction of widely linear systems by
using the innovations algorithm,” IEEE T. Signal Proces., vol. 56, no.
7-2, pp. 3061–3068, 2008.
[20] J. Navarro-Moreno, J. Moreno-Kaiser, R. M. Fernández-Alcalá, and J. C.
Ruiz-Molina, “Widely linear prediction for transfer function models
based on the infinite past,” Comput. Stat. Data An., vol. 58, pp. 139–146,
2013.
[21] J. Navarro-Moreno, M. D. Estudillo-Martı́nez, R. M. Fernández-Alcalá,
and J. C. Ruiz-Molina, “Estimation of improper complex-valued random
signals in colored noise by using the Hilbert space theory,” IEEE T.
Inform. Theory, vol. 55, no. 6, pp. 2859–2867, 2009.
[22] P. Rubin-Delanchy and A. T. Walden, “Simulation of improper complexvalued sequences,” IEEE T. Signal Proces., vol. 55, no. 11, pp. 5517–
5521, 2007.
[23] ——, “Kinematics of complex-valued time series,” IEEE T. Signal
Proces., vol. 56, no. 9, pp. 4189–4198, 2008.
[24] M. Arató, A. N. Kolmogorov, and Y. G. Sinai, “Estimation of the
parameters of a complex stationary Gaussian Markov process,” Dokl.
Akad. Nauk SSSR, vol. 146, pp. 747–750, 1962.

11

[25] H. Wold, A study in the analysis of stationary time series. Almqvist
& Wiksell, 1954.
[26] P. Whittle, “Estimation and information in stationary time series,” Ark.
Mat., vol. 2, no. 5, pp. 423–434, 1953.
[27] A. M. Sykulski, S. C. Olhede, J. M. Lilly, and J. J. Early, “Stochastic
modeling and estimation of stationary complex-valued signals,” arXiv
preprint: 1306.5993, 2016.
[28] A. M. Sykulski, S. C. Olhede, and J. M. Lilly, “The de-biased Whittle likelihood for second-order stationary stochastic processes,” arXiv
preprint: 1605.06718, 2016.
[29] A. M. Sykulski, S. C. Olhede, J. M. Lilly, and E. Danioux, “Lagrangian
time series models for ocean surface drifter trajectories,” J. R. Statist.
Soc. C, vol. 65, no. 1, pp. 29–50, 2016.
[30] A. T. Walden and R. E. White, “Estimating the statistical bandwidth of
a time series,” Biometrika, vol. 77, no. 4, pp. 699–707, 1990.
[31] S. Adak, “Time-dependent spectral analysis of nonstationary time series,” J. Am. Stat. Assoc., vol. 93, no. 444, pp. 1488–1501, 1998.
[32] M. Fuentes, “Spectral methods for nonstationary spatial processes,”
Biometrika, vol. 89, no. 1, pp. 197–210, 2002.
[33] J. M. Lilly and S. C. Olhede, “Bivariate instantaneous frequency and
bandwidth,” IEEE T. Signal Proces., vol. 58, no. 2, pp. 591–603, 2010.
[34] Y. Benjamini and Y. Hochberg, “Controlling the false discovery rate: a
practical and powerful approach to multiple testing,” J. R. Statist. Soc.
B, pp. 289–300, 1995.
[35] G. Rünstler, “Modelling phase shifts among stochastic cycles,”
Economet. J., vol. 7, no. 1, pp. 232–248, 2004.
[36] E. Zivot and J. Wang, Modeling financial time series with S-Plus.
Springer Science & Business Media, 2007.

