dp08 36 PDF

Download as pdf or txt
Download as pdf or txt
You are on page 1of 40

DEPARTMENT OF ECONOMICS

STATISTICAL FOURIER ANALYSIS:


CLARIFICATIONS AND INTERPRETATIONS









D.S.G. Pollock, University of Leicester, UK











Working Paper No. 08/36
October 2008


STATISTICAL FOURIER ANALYSIS:
CLARIFICATIONS AND INTERPRETATIONS
by D.S.G. Pollock (University of Leicester)
Email: stephen [email protected]
This paper expounds some of the results of Fourier theory that are es-
sential to the statistical analysis of time series. It employs the algebra of
circulant matrices to expose the structure of the discrete Fourier transform
and to elucidate the ltering operations that may be applied to nite data
sequences.
An ideal lter with a gain of unity throughout the pass band and a
gain of zero throughout the stop band is commonly regarded as incapable
of being realised in nite samples. It is shown here that, to the contrary,
such a lter can be realised both in the time domain and in the frequency
domain.
The algebra of circulant matrices is also helpful in revealing the nature
of statistical processes that are band limited in the frequency domain. In
order to apply the conventional techniques of autoregressive moving-average
modelling, the data generated by such processes must be subjected to anti-
aliasing ltering and sub sampling. These techniques are also described.
It is argued that band-limited processes are more prevalent in statis-
tical and econometric time series than is commonly recognised.
1. Introduction
Statistical Fourier analysis is an important part of modern time-series analysis,
yet it frequently poses an impediment that prevents a full understanding of
temporal stochastic processes and of the manipulations to which their data are
amenable. This paper provides a survey of the theory that is not overburdened
by inessential complications, and it addresses some enduring misapprehensions.
Amongst these is a misunderstanding of the eects in the frequency domain
of linear ltering operations. It is commonly maintained, for example, that an
ideal frequency-selective lter that preserves some of the elements of a time
series whilst nullifying all others is incapable of being realised in nite samples.
This paper shows that such nite-sample lters are readily available and
that they possess time-domain and frequency-domain representations that are
both tractable. Their representations are directly related to the classical
WienerKolmogorov theory of ltering, which presupposes that the data can
be treated as if they form a doubly-innite or a semi-innite sequence.
A related issue, which the paper aims to clarify, concerns the nature of
continuous stochastic processes that are band-limited in frequency. The paper
provides a model of such processes that challenges the common supposition
that the natural primum mobile of any continuous-time stochastic process is
Wiener process, which consists of a steam of innitesimal impulses.
The sampling of a Wiener process is inevitably accompanied by an aliasing
eect, whereby elements of high frequencies are confounded with elements of
1
D.S.G. POLLOCK: Statistical Fourier Analysis
frequencies that fall within the range that is observable in sampling. Here, it
is shown, to the contrary, that there exist simple models of continuous band-
limited processes for which no aliasing need occur in sampling.
The spectral theory of time series is a case of a non-canonical Fourier the-
ory. It has some peculiarities that originally caused considerable analytic di-
culties, which were rst overcome in a satisfactory manner by Norbert Wiener
(1930) via his theory of a generalised harmonic analysis. The statistical ltering
theory that is appropriate to such time series originated with Wiener (1941)
and Kolmogorov (1941a).
Many of the diculties of the spectral representation of a time series can
be avoided by concentrating exclusively on its correlation structure. Then,
one can exploit a simple canonical Fourier relationship that exists between the
autocovariance function of a discrete stationary process and its spectral density
function.
This relationship is by virtue of an inverted form of the classical Fourier
series relationship that exists between a continuous, or piecewise continuous,
periodic function and its transform, which is a sequence of Fourier coecients.
(In the inverted form of the relationship, which is described as the discrete-
time Fourier transform, it is the sequence that is the primary function and the
continuous periodic function that is its transform.) However, it is important to
have a mathematical model of the process itself, and this is where some of the
complications arise.
Figure 1 depicts what may be described as the canonical Fourier trans-
forms. For conciseness, they have been expressed using complex exponentials,
whereas they might be more intelligible when expressed in terms of ordinary
trigonometrical functions. When dealing with the latter, only positive frequen-
cies are considered, whereas complex exponentials are dened for both positive
and negative frequencies.
In the diagrams, the real-valued time-domain functions are symmetric,
with the consequence that their frequency-domain transforms are also real-
valued and symmetric. Under other conditions, they would be complex-valued,
giving rise to a pair of functions.
2. Canonical Fourier Analysis
The classical Fourier series, illustrated in Figure 1.(ii), expresses a continuous,
or piecewise continuous, period function x(t) = x(t+T) in terms of a sum of sine
and cosine functions of harmonically increasing frequencies {
j
= 2j/T; j =
1, 2, 3, . . .}:
x(t) =
0
+

j=1

j
cos(
j
t) +

j=1

j
sin(
j
t)
=
0
+

j=1

j
cos(
j
t
j
).
(1)
Here,
1
= 2/T is the fundamental frequency or angular velocity, which corre-
sponds to a trigonometric function that completes a single cycle in the interval
2
D.S.G. POLLOCK: Statistical Fourier Analysis
(i) The Fourier integral:
x(t) =
1
2
_

()e
it
d () =
_

x(t)e
it
dt
(ii) The classical Fourier series:
x(t) =

j=

j
e
i
j
t

j
=
1
T
_
T
0
x(t)e
i
j
t
dt
(iii) The discrete-time Fourier transform:
x
t
=
1
2
_

()e
it
d () =

t=
x
t
e
it
(iv) The discrete Fourier transform:
x
t
=
T1

j=0

j
e
i
j
t

j
=
1
T
T1

t=0
x
t
e
i
j
t
Figure 1. The classes of the Fourier transforms.
3
D.S.G. POLLOCK: Statistical Fourier Analysis
[0, T], or in any other interval of length T such as [T/2, T/2]. The second
expression depends upon the denitions

2
j
=
2
j
+
2
j
and
j
= tan
1
(
j
/
j
). (2)
The equality follows from the trigonometrical identity
cos(AB) = cos(A) cos(B) + sin(A) sin(B). (3)
In representing the periodic nature of x(t), it is helpful to consider mapping
the interval [0, T], over which the function is dened, onto the circumference
of a circle, such that the end points of the interval coincide. Then, successive
laps or circuits of the circle will generate successive cycles of the function.
According to Eulers equations, there are
cos(
j
t) =
1
2
(e
i
j
t
+ e
i
j
t
) and sin(
j
t) =
i
2
(e
i
j
t
e
ijt
). (4)
Therefore, equation (1) can be expressed as
x(t) =
0
+

j=1

j
+ i
j
2
e
i
j
t
+

j=1

j
i
j
2
e
i
j
t
, (5)
which can be written concisely as
x(t) =

j=

j
e
i
j
t
, (6)
where

0
=
0
,
j
=

j
i
j
2
and
j
=

j
=

j
+ i
j
2
. (7)
The inverse of the classical transform of (6) is demonstrated by writing

j
=
1
T
_
T
0
x(t)e
i
j
t
dt =
1
T
_
T
0
_

k=

k
e
i
k
t
_
e
i
j
t
dt
=
1
T

k=

k
_
_
T
0
e
i(
k

j
)t
dt
_
=
j
,
(8)
where the nal equality follows from an orthogonality condition in the form of
_
T
0
e
i(
k

j
)t
dt =
_
0, if j = k;
T, if j = k.
(9)
The relationship between the continuous periodic function and its Fourier trans-
form can be summarised by writing
x(t) =

j=

j
e
i
j
t

j
=
1
T
_
T
0
x(t)e
i
j
t
dt. (10)
4
D.S.G. POLLOCK: Statistical Fourier Analysis
The question of the conditions that are sucient for the existence of such
a Fourier relationship is an essential one; and there are a variety of Fourier the-
ories. The theory of Fourier series is concerned with establishing the conditions
under which the partial sums converge to the function, in some specied sense,
as the number of the terms increases.
At the simplest level, it is sucient for convergence that x(t) should be
continuous and bounded in the interval [0, T]. However, the classical theory
of Fourier series is concerned with the existence of the relationship in the case
where x(t) is bounded but is also permitted to have a nite number of maxima
and minima and a nite number of jump discontinuities. It can be shown that,
in that case, as successive terms are added, the Fourier series of (1) converges
to
1
2
{x(t + 0) + x(t 0)}, (11)
where x(t + 0) is the value as t is approached from the right and x(t 0) is
the value as it is approached from the left. If x(t) is continuous at the point in
question, then the Fourier series converges to x(t).
Here, the relevant criterion is that of overall convergence in mean-square
rather than of pointwise convergence. The question of the mode of converge was
the issue at stake in the paradox known as Gibbs phenomenon, which concerns
the absence of pointwise convergence in the presence jump discontinuities.
This phenomenon is illustrated in Figure 2, where it is apparent that not
all of the oscillations in the partial sums that approximate a periodic square
wave (i.e. a time-domain rectangle) are decreasing at a uniform rate as the
number of terms n increases. Instead, the oscillations that are adjacent to the
point of discontinuity are tending to a limiting amplitude, which is about 9% of
the jump. However, as n increases, the width of these end-oscillations becomes
vanishingly small; and thus the mean-square convergence of the Fourier series
is assured. Gibbs phenomenon is analysed in detail by Carslaw (1930); and
Carslaw (1925) has also recounted the history of its discovery.
Parsevals theorem asserts that, under the stated conditions, which guar-
antee mean-square convergence, there is
1
T
_
T
0
|x(t)|
2
dt =
2
0
+
1
4

j=1
(
2
j
+
2
j
) =

j=
|
j
|
2
, (12)
where |
j
|
2
=
j

j
=
2
j
/4. This indicates that the energy of the function x(t),
which becomes its average power if we divide by T, is equal to the sum of
the energies of the sinusoidal components. Here, it is essential that the square
integral converges, which is to say that the function possesses a nite energy.
The fullment of such an energy or power condition characterises what we have
chosen to describe as the canonical Fourier transforms.
There are some powerful symmetries between the two domains of the
Fourier transforms. The discrete-time Fourier transform, which is fundamental
to time-series analysis, is obtained by interchanging the two domains of the
classical Fourier series transform. It eects the transformation of a sequence
x(t) = {x
t
; t = 0, 1, 2, . . .} within the time domain, that is square summable,
5
D.S.G. POLLOCK: Statistical Fourier Analysis
0.00
0.25
0.50
0.75
1.00
1.25
0.00
0 T/8 T/4 3T/8 T/2
n = 30
n = 6
n = 3
Figure 2. The Fourier approximation of a square wave. The diagram shows
approximations to the positive half of a rectangle dened on the interval
[T/2, T/2] and symmetric with respect to the vertical axis through zero.
into a continuous periodic function in the frequency domain via an integral ex-
pression that is mean-square convergent. It is illustrated in Figure 1.(iii). We
may express the relationship in question by writing
x
t
=
1
2
_

()e
it
d () =

t=
x
t
e
it
. (13)
The periodicity is now in the frequency domain such that () = (+2)
for all . The periodic function completes a single cycle in any interval of
length 2; and, sometimes, it may be appropriate to dene the function over
the interval [0, 2] instead of the interval [, ]. In that case, the relationship
of (13) will be unaected; and its correspondence with (10), which denes the
classical Fourier series, is claried.
There remain two other Fourier transforms, which embody a complete
symmetry between the two domains. The rst of these is the Fourier integral
transform of Figure 1.(i) which is denoted by
x(t) =
1
2
_

()e
it
d () =
_

x(t)e
it
dt. (14)
Here, apart from the scaling, which can be revised to achieve symmetry, there
is symmetry between the transform and its inverse. (The factor 1/2 can be
eliminated from the rst integral by changing the variable from angular velocity
, measured in radians per period, to frequency f = /2, measured in cycles
per period.)
6
D.S.G. POLLOCK: Statistical Fourier Analysis
One of the surprises on rst encountering the Fourier integral is the ape-
riodic nature of the transform, which is in spite of the cyclical nature of the
constituent sinusoidal elements. However, this is readily explained. Consider
the sum of two elements within a Fourier synthesis that are at the frequencies

m
and
n
respectively. Then, their sum will constitute a periodic function of
period = 2/

if and only if
m
= m

and
n
= n

are integer multiples


of a common frequency

.
For a sum of many sinusoidal elements to constitute a periodic function, it
is necessary and sucient that all of the corresponding frequencies should be
integer multiples of a fundamental frequency. Whereas this condition is fullled
by the classical Fourier series of (10), it cannot be fullled by a Fourier integral
comprising frequencies that are arbitrary irrational numbers.
A sequence that has been sampled at the integer time points from a con-
tinuous aperiodic function that is square-integrable will have a transform that
is a periodic function. This outcome is manifest in the discrete-time Fourier
transform, where it can be seen that the period in question has the length of
2 radians. Let () be the transform of the continuous aperiodic function,
and let
S
() be the transform of the sampled sequence {x
t
; t = 0, 1, 2, . . .}.
Then
x
t
=
1
2
_

()e
it
d =
1
2
_

S
()e
it
d. (15)
The equality of the two integrals implies that

S
() =

j=
( + 2j). (16)
Thus, the periodic function
S
() is obtained by wrapping () around a circle
of circumference of 2 and adding the coincident ordinates.
The relationship between the Fourier series and the discrete-time Fourier
transform enables us to infer that a similar eect will arise from the regular
sampling of a continuous function of frequency. Thus, the corresponding time-
domain function will be wrapped around a circle of circumferences T = 2/
1
,
where
1
is both the fundamental frequency and the sampling interval in the
frequency domain.
The Fourier integral is employed extensively in mathematical physics. In
quantum mechanics, for example, it is used in one domain to describe a localised
wave train and, in the other domain, to describe its frequency composition,
which is a resolution of its energy amongst a set of frequencies. There is an
inverse relationship between the dispersion of the wave train in space and the
dispersion of its frequency components. The product of the two dispersions is
bounded below according to the famous Heisenberg uncertainty principle.
A simple and important example of the Fourier integral is aorded by
the sinc function wave packet and its transform, which is a frequency-domain
rectangle. Figure 3 depicts a continuous sinc function of which the Fourier
transform is a rectangle on the frequency interval [, ]:

0
(t) =
1
2
_

e
it
d(t) =
_
e
it
i2t
_

=
sin t
t
. (17)
7
D.S.G. POLLOCK: Statistical Fourier Analysis
0
0.25
0.5
0.75
1
0
0.25
0 2 4 6 8 0 2 4 6 8
Figure 3. The sinc function wave-packet
0
(t) = sin(t)/t comprising frequencies
in the interval [0, ].
Here, the nal equality is by virtue of (4), which expresses a sine function as a
combination of complex exponential functions.
The function
0
(t) can also be construed as a wave packet centred on
time t = 0. The gure also represents a sampled version the sinc function.
This would be obtained, in the manner of a classical Fourier series, from the
rectangle on the interval [, ], if this were regarded as a single cycle of a
periodic function.
The function
0
(t) with t I = {0, 1, 2, . . .} is nothing but the unit
impulse sequence. Therefore, the set of all sequences {
0
(t k); t, k I}, ob-
tained by integer displacements k of
0
(t), constitutes the ordinary orthogonal
Cartesian basis in the time domain for the set of all real-valued time series.
When t R is a real-valued index of continuous time, the set of displaced
sinc functions {
0
(t k); t R, k I} constitute a basis for the set of contin-
uous functions of which the frequency content is bounded by the Nyquist value
of radians per unit time interval. In common with their discretely sampled
counterparts, the sequence of continuous sinc functions at integer displacements
constitutes an orthogonal basis.
To demonstrate the orthogonality, consider the fact that the correspond-
ing frequency-domain rectangle is an idempotent function. When multiplied
by itself it does not change. The time-domain operation corresponding to this
frequency-domain multiplication is an autoconvolution. The time-domain func-
tions are real and symmetric with
0
(kt) =
0
(tk), so their autoconvolution
is the same as their autocorrelation:
(k) =
_

0
(t)
0
(k t)dt =
_

0
(t)
0
(t k)dt. (18)
Therefore, the sinc function is its own autocovariance function: (k) =
0
(k).
The zeros of the sinc function that are found at integer displacements from the
centre correspond the orthogonality of sinc functions separated from each other
by these distances.
8
D.S.G. POLLOCK: Statistical Fourier Analysis
Table 1. The classes of Fourier transforms*
Periodic Aperiodic
Continuous Discrete aperiodic Continuous aperiodic
Fourier series Fourier integral
Discrete Discrete periodic Continuous periodic
Discrete FT Discrete-time FT
* The class of the Fourier transform depends upon the nature of the function
which is transformed. This function may be discrete or continuous and it
may be periodic or aperiodic. The names and the natures of corresponding
transforms are shown in the cells of the table.
A sinc function wave packet that is limited to the frequency band [, ]
[0, ], if we are talking only of positive frequencies, has the functional form of
(t) =
1
t
{sin(t) sin(t)}
=
2
t
cos{( + )t/2} sin{( )t/2}
=
2
t
cos(t) sin(t),
(19)
where = (+)/2 is the centre of the band and = ()/2 is half its width.
The equality follows from the identity sin(A+B) sin(AB) = 2 cos Asin B.
Finally, we must consider the Fourier transform that is the most important
one from the point of view of this paper. This is the discrete Fourier transform
of Figure 1.(iv) that maps from a sequence of T data {x
t
; t = 0, 1, . . . , T 1}
in the time domain to a set of T ordinates {
j
; j = 0, 1, . . . , T 1} in the
frequency domain at the points {
j
= 2j/T; j = 0, 1, . . . , T 1}, which are
equally spaced within the interval [0, 2], and vice versa. The relationship can
be expressed in terms of sines and cosines; but it is more conveniently expressed
in terms of complex exponentials, even in the case where the sequence {x
t
} is
a set of real data values:
x
t
=
T1

j=0

j
e
i
j
t

j
=
1
T
T1

t=0
x
t
e
i
j
t
. (20)
Since all sequences in data analysis are nite, the discrete Fourier transform is
used, in practice, to represent all manner of Fourier transforms.
Although the sequences in both domains are nite and of an equal number
of elements, it is often convenient to regard both of them as representing single
cycles of periodic functions dened over the entire set of positive and negative
integers. These innite sequences are described and the periodic extensions of
their nite counterparts.
The various Fourier transforms may be summarised in a table that records
the essential characteristics of the function that is to be transformed, namely
9
D.S.G. POLLOCK: Statistical Fourier Analysis
whether it is discrete or continuous and whether it is periodic or aperiodic.
This is table 1.
There are numerous accounts of Fourier analysis that can be accessed in
pursuit of a more thorough treatment. Much of the text of Carslaw (1930)
is devoted to the issues in real analysis that were raised in the course of the
development of Fourier analysis. The book was rst published in 1906, which
is close to the dates of some of the essential discoveries that claried such mat-
ters as Gibbs phenomenon. It continues to be a serviceable and authoritative
text. Titchmarsh (1937), from the same era, provides a classical account of the
Fourier integral.
The more recent text of Kreider, Kuller, Ostberg and Perkins (1966) deals
with some of the analytic complications of Fourier analysis in an accessible
manner, as does the brief text of Churchill and Brown. The recent text of
Stade (2005) is an extensive resource. That of Nahin (2006) is both engaging
and insightful.
Many of the texts that are intended primarily for electrical engineers are
also useful to statisticians and econometricians. Baher (1990) gives much of
the relevant analytic details, whereas Brigham (1988) provides a practical text
that avoids analytic complications.
Following the discoveries of Wiener (1930) and Kolmogorov (1941b), the
mathematical exposition of the spectral analysis of time series developed
rapidly. Detailed analyses of the analytic issues are to be found in some of the
classical texts of time series analysis, amongst which are those of Wold (1938),
Doob (1953), Grenander and Rosenblat (1957), Hannan (1960), Yaglom (1962)
and Rozanov (1967). Yaglom and Rozanov were heirs to a Russian tradition
that began with Kolmogorov (1941b)
It was not until the 1960s that spectral analysis began to nd widespread
application; and the book of Jenkins and Watts (1968) was symptomatic of
the increasing practicality of the subject. Econometricians also began to take
note of spectral analysis in the 1960s; and two inuential books that brought
it to their attention were those of Granger and Hatanaka (1964) and Fishman
(1969). Nerlove, Grether and Carvalho (1979) continued the pursuit of the
spectral analysis of economic data.
Latterly, Brillinger (1975), Priestly (1981), Rosenblatt (1985) and Brock-
well and Davis (1987) have treated Fourier analysis with a view to the statistical
analysis of time series, as has Pollock (1999).
3. Representations of the Discrete Fourier Transform
It is often helpful to express the discrete Fourier transform in terms of the nota-
tion of the z-transform. The z-transforms of the sequences {x
t
; t = 0, . . . , T 1}
and {
j
; j = 0, . . . , T 1} are the polynomials
(z) =
T1

j=0

j
z
j
and x(z) =
T1

t=0
x
t
z
t
, (21)
wherein z is an indeterminate algebraic symbol that is commonly taken to be a
10
D.S.G. POLLOCK: Statistical Fourier Analysis
complex number, in accordance with the fact that the equations (z) = 0 and
x(z) = 0 have solutions within the complex plane.
The polynomial equation z
T
= 1 is important in establishing the connec-
tion between the z-transform of a sequence of T elements and the corresponding
discrete Fourier transform. The solutions of the equation are the complex num-
bers
W
j
T
= exp
_
i2j
T
_
= cos
_
2j
T
_
i sin
_
2j
T
_
; j = 0, 1, . . . , T 1. (22)
These constitute a set of points, described as the T roots of unity, that are
equally spaced around the circumference of the unit circle at angles of
j
=
2j/T radians. Multiplying any of these angles by T will generate a multiple
of 2 radians, which coincides with an angle of zero radians, which is the angle,
or argument, attributable to unity within the complex plane.
Using this conventional notation for the roots of unity, we can write the
discrete Fourier transform of (20) and its inverse as

j
=
1
T
T1

t=0
x
t
W
jt
T
and x
t
=
T1

j=0

j
W
jt
T
. (23)
The advantage of the z-transforms is that they enable us to deploy the
ordinary algebra of polynomials and power series to manipulate the objects of
Fourier analysis. In econometrics, it is common to replace z by the so-called
lag operator L, which operates on doubly-innite sequences and which has the
eect that Lx(t) = x(t 1).
It is also useful to replace z by a variety of matrix operators. The matrix
lag operator L
T
= [e
1
, e
2
, . . . , e
T1
, 0] is derived from the identity matrix I
T
=
[e
0
, e
1
, . . . , e
T1
] of order T by deleting its leading column and appending a
column of zeros to the end of the array.
Setting z = L
T
within the polynomial x(z) = x
0
+ x
1
z + + x
T1
z
T1
gives rise to a banded lower-triangular Toeplitz matrix of order T. The matrices
I
T
= L
0
T
, L
T
, L
2
T
, . . . , L
T1
T
form a basis for the vector space comprising all such
matrices. The ordinary polynomial algebra can be applied to these matrices
with the provision that the argument L
T
is nilpotent of degree T, which is to
say that L
q
T
= 0 for all q T.
Circulant Matrices
A matrix argument that has properties that more closely resemble those
of the complex exponentials is the circulant matrix K
T
= [e
1
, . . . , e
T1
, e
0
],
which is formed by carrying the leading column of the identity matrix I
T
to
the back of the array. This is an orthonormal matrix of which the transpose is
the inverse, such that K

T
K
T
= K
T
K

T
= I
T
.
The powers of the matrix form a T-periodic sequence such that K
T+q
T
=
K
q
T
for all q. The periodicity of these powers is analogous to the periodicity
of the powers of the argument z = exp{i2/T}, which is to be found in the
Fourier transform of a sequence of T elements.
11
D.S.G. POLLOCK: Statistical Fourier Analysis
The matrices K
0
T
= I
T
, K
T
, . . . , K
T1
T
form a basis for the set of all circu-
lant matrices of order Ta circulant matrix X = [x
ij
] of order T being dened
as a matrix in which the value of the generic element x
ij
is determined by the
index {(i j) mod T}. This implies that each column of X is equal to the pre-
vious column rotated downwards by one element. The generic circulant matrix
has a form that may be illustrated follows:
X =
_

_
x
0
x
T1
x
T2
. . . x
1
x
1
x
0
x
T1
. . . x
2
x
2
x
1
x
0
. . . x
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
x
T1
x
T2
x
T3
. . . x
0
_

_
. (24)
There exists a one-to-one correspondence between the set of all polynomials
of degree less than T and the set of all circulant matrices of order T. Thus,
if x(z) is a polynomial of degree less that T, then there exits a corresponding
circulant matrix
X = x(K
T
) = x
0
I
T
+ x
1
K
T
+ + x
T1
K
T1
T
. (25)
A convergent sequence of an indenite length can also be mapped into
a circulant matrix. If {c
i
} is an absolutely summable sequence obeying the
condition that

|c
i
| < , then the z-transform of the sequence, which is
dened by c(z) =

c
j
z
j
, is an analytic function on the unit circle. In that
case, replacing z by K
T
gives rise to a circulant matrix C

= c(K
T
) with
nite-valued elements. In consequence of the periodicity of the powers of K
T
,
it follows that
C

=
_

j=0
c
jT
_
I
T
+
_

j=0
c
(jT+1)
_
K
T
+ +
_

j=0
c
(jT+T1)
_
K
T1
T
= c

0
I
T
+ c

1
K
T
+ + c

T1
K
T1
T
.
(26)
Given that {c
i
} is a convergent sequence, it follows that the sequence of
the matrix coecients {c

0
, c

1
, . . . , c

T1
} converges to {c
0
, c
1
, . . . , c
T1
} as T
increases. Notice that the matrix c

(K
T
) = c

0
I
T
+ c

1
K
T
+ + c

T1
K
T1
T
,
which is derived from a polynomial c

(z) of degree T 1, is a synonym for the


matrix c(K
T
), which is derived from the z-transform of an innite convergent
sequence.
The polynomial representation is enough to establish that circulant matri-
ces commute in multiplication and that their product is also a polynomial in
K. That is to say
If X = x(K
T
) and Y = y(K
T
) are circulant matrices,
then XY = Y X is also a circulant matrix.
(27)
The matrix operator K
T
has a spectral factorisation that is particularly
useful in analysing the properties of the discrete Fourier transform. To demon-
strate this factorisation, we must rst dene the so-called Fourier matrix. This
12
D.S.G. POLLOCK: Statistical Fourier Analysis
is a symmetric matrix U = T
1/2
[W
jt
; t, j = 0, . . . , T 1], of which the generic
element in the jth row and tth column is W
jt
, where W = exp{i2/T} is the
rst root of unity, from which we are now omitting the subscripted T for ease
of notation. On taking account of the T-periodicity of W
q
, the matrix can be
written explicitly as
U =
1

T
_

_
1 1 1 . . . 1
1 W W
2
. . . W
T1
1 W
2
W
4
. . . W
T2
.
.
.
.
.
.
.
.
.
.
.
.
1 W
T1
W
T2
. . . W
_

_
. (28)
The second row and the second column of this matrix contain the T roots of
unity. The conjugate matrix is dened as

U = T
1/2
[W
jt
; t, j = 0, . . . , T 1];
and, by using W
q
= W
Tq
, this can be written explicitly as

U =
1

T
_

_
1 1 1 . . . 1
1 W
T1
W
T2
. . . W
1 W
T2
W
T4
. . . W
2
.
.
.
.
.
.
.
.
.
.
.
.
1 W W
2
. . . W
T1
_

_
. (29)
The matrix U is a unitary, which is to say that it fulls the condition

UU = U

U = I. (30)
The operator K, from which we now omit the subscript, can be factorised
as
K =

UDU = U

D

U, (31)
where
D = diag{1, W, W
2
, . . . , W
T1
} (32)
is a diagonal matrix whose elements are the T roots of unity, which are found
on the circumference of the unit circle in the complex plane. Observe also that
D is T-periodic, such that D
q+T
= D
q
, and that K
q
=

UD
q
U = U

D
q

U for any
integer q. Since the powers of K form the basis for the set of circulant matrices,
it follows that all circulant matrices are amenable to a spectral factorisation
based on (31).
Circulant matrices have represented a mathematical curiosity ever since
their rst appearance in the literature in a paper by Catalan (1846). The
literature on circulant matrices, from their introduction until 1920, was sum-
marised in four papers by Muir (1911)(1923). A recent treatise on the subject,
which contains a useful bibliography, has been provided by Davis (1979); but
his book does not deal with problems in time-series analysis. An up-to-date
account, orientated towards statistical signal processing, has been provided by
Gray (2002).
13
D.S.G. POLLOCK: Statistical Fourier Analysis
The Matrix Discrete Fourier Transform
The matrices U and

U are entailed in the discrete Fourier transform. Thus,
the equations of (20) or (23) can be written as
= T
1/2
Ux = T
1
Wx and x = T
1/2

U =

W, (33)
where x = [x
0
, x
1
, . . . x
T1
]

and = [
0
,
1
, . . . ,
T1
]

, and where W = [W
jt
]
and

W = [W
jt
] are the matrices of (28) and (29) respectively, freed from their
scaling factors. We may note, in particular, that
= T
1/2
Ue
0
= T
1/2

Ue
0
and e
0
= T
1/2
U = T
1/2

U, (34)
where e
0
= [1, 0, . . . 0]

is the leading column of the identity matrix and =


[1, 1, . . . , 1]

is the summation vector. Thus, a spike or an impulse located at


t = 0 in the time domain or at j = 0 in the frequency domain is transformed
into a constant sequence in the other domain and vice versa. These are extreme
examples of the inverse relationship between the dispersion of a sequence and
that of its transform.
Given the circulant the matrix X = x(K) dened in (25), the discrete
Fourier transform and its inverse may be represented by the following spectral
factorisations:
X =

Ux(D)U and x(D) = UX

U. (35)
Here, x(D) = Tdiag{
0
,
1
, . . . ,
T1
} is a diagonal matrix containing the spec-
tral ordinates of the data. The equations of (33) are recovered by post multi-
plying the X by e
0
and x(D) by and using the results of (34).
Assuming that the data are mean adjusted, i.e. that

X = 0, the cross
product of the matrix X gives rise to the matrix C

of the circular autocovari-


ances of the data. Thus
T
1
X

X = C

=

Uc

(D)U
= T
1

Ux(

D)x(D)U,
(36)
where we have used X

= Ux(D)

U =

Ux(

D)U.
The values of the elements of [c

0
, c

1
, . . . , c

T1
]

= c

= C

e
0
are given by
the formula
c

=
1
T
T1

t=0
x
t
x
t+
; where x
t
= x
(t mod T)
or, equivalently,
c

=
1
T
T1

t=0
x
t
x
t+
+
1
T
1

t=0
x
t
x
t+T
= c

+ c
T
,
(37)
where c

is the ordinary empirical autocovariance of lag .


The vector c

of circular autocovariances is obtained via a circular corre-


lation of the data vector x. One can envisage the data disposed at T equally-
spaced points around the circumferences of two discs on a common axis. The
14
D.S.G. POLLOCK: Statistical Fourier Analysis
1
2
3
4
5
6
7
8
0
0
0
0
1
2
3
4
5
6
7
8
0
0
0
0
Figure 4. A device for nding the circular correlation of two sequences.
The upper disc is rotated clockwise through successive angles of 30 degrees.
Adjacent numbers on the two discs are multiplied and the products are
summed and divided by the sample size T to obtain the circular covariances.
same numbers on the two discs are aligned and products are formed and added.
This generates the sum of squares. Then, the upper disc is rotated by an angle
of 2/T radians and the displaced numbers are multiplied together and added.
The process is repeated. The T values formed via a complete rotation are each
divided by T, and the result is the sequence of circular autocovariances. This
device is illustrated in Figure 4, where T = 12 and where identical sequences
on the two discs end with four zeros.
The core of the factorisation of the matrix of C

of (36) is the real-valued


diagonal matrix
c

(D) = T
1
x(

D)x(D) = TDiag{|
0
|
2
, |
1
|
2
, . . . , |
T1
|
2
}, (38)
of which the graph of the elements is described as the periodogram. (The
frequencies indexed by j = 0, . . . , T 1 extend from 0 to 2(T 1)/T. However,
for real-valued data, there is
j
=
Tj
, and so it is customary to plot the
periodogram ordinates only over a frequency range from 0 to , as in Figure
14.)
Since e

0
C

e
0
= c

0
= c
0
and Ue
0
=

Ue
0
= T
1/2
, we get the following
expression for the variance:
c
0
= T
1
T1

j=0
x
2
t
= T
2

x(

D)x(D) =
T1

j=0
|
j
|
2
.
(39)
This is the discrete version of Parsevals theorem. In statistical terms, it rep-
resents a frequency-specic analysis of variance.
15
D.S.G. POLLOCK: Statistical Fourier Analysis
4. Fourier Analysis of Temporal Sequences
It is clear that the discrete-time Fourier transform is inappropriate to a doubly-
innite stationary stochastic sequence y(t) = {y
t
; t = 0, 1, 2, . . .} dened of
over the set of positive and negative integers. Such a sequence, which has
innite energy, is not summable.
In the appendix, it is demonstrated there there exists a spectral represen-
tation of y(t) in the form of
y(t) =
_

e
it
dZ(), (40)
where Z() is continuous non-dierentiable complex-valued stochastic process
dened on the interval [, ]. The increments dZ(), dZ() of the process
are uncorrelated for all = . However, the statistical expectation
E{dZ()dZ

()} = dF() = f()d (41)


comprises a function F(), known as the spectral distribution function, which
is continuous and everywhere dierentiable in the interval, provided that y(t)
contains no perfectly regular components of nite power. If there are regular
components within y(t), then F() will manifest jump discontinuities at the
corresponding frequencies.
(The emissions of light produced when electrons jump between the adjacent
atomic orbits produce discrete spectra, with distribution functions resembling
staircases. Koopmans (1974, p. 49) contends that continuous spectra with
superimposed discrete components (i.e. mixed spectra) are amongst the most
common types. However, this is hardly the case for the spectra of econometric
data. Apart from the stylised representation of economic index with a perfectly
regular seasonal component, the spectra of economic data are of the purely
continuous variety.)
The discrete-time transform comes into eect when it is applied not to the
data sequence itself but, instead, to its autocovariance function
() = {
t
; = 0, 1, 2, . . .} where

= E(y
t
y
t
). (42)
Such an autocovariance sequence is absolutely summablei.e. there is

| <
and, therefore, a fortiori, its transform is mean-square convergent. The
Fourier transform of the autocovariance function is the spectral density function
or power spectrum, which, in view of the symmetry of (), can be expressed
as
f() =
1
2
_

0
+

=1

(e
i
+ e
i
)
_
=
1
2
_

0
+ 2

=1

cos()
_
.
(43)
16
D.S.G. POLLOCK: Statistical Fourier Analysis
1 2 3 4
1.0
0.5
0.5
1.0
Figure 5. The values of the function cos{(11/8)t} coincide with those
of its alias cos{(5/8)t} at the integer points {t = 0, 1, 2, . . .}.
Notice, however, that compared with equation (14), the scaling factor 1/2 has
migrated from the Fourier integral to the series. Now, the inverse transform
takes the form of

=
_

e
i
f()d. (44)
Setting = 0 gives
0
=
_
f()d, which indicates that the power of the
process, which is expressed by its variance, is equal to the integral of the spectral
density function over the interval [, ].
This a form of Parsevals theorem. The sequence y(t) is not square sum-
mable and, therefore, its power is cannot be expressed in terms of a sum of
squares. Instead it is expressed, in the time domain, via the statistical moment

0
= E(y
2
t
).
In the frequency domain, this power is attributed to an innite set of sinu-
soidal components, of which the measure of power is provided by the continuous
spectral density function f(). One can see, by a direct analogy with a contin-
uous probability density function, that the power associated with a component
at a particular frequency is vanishingly small, i.e. of measure zero.
Here, we are seeing some of the perplexities that are associated with
Wieners generalised harmonic analysis and with the spectral representation
of a stochastic process. These matters are dealt with in the appendix.
The Problem of Aliasing
Notice that the (positive) frequencies in the spectral analysis of a discrete-
time process are limited to the interval [0, ], bounded by the so-call Nyquist
frequency . The process of sampling imposes a limit of the observability of
high-frequency components. To be observable in sampled data, a sinusoidal
motion must take no less than two sample periods to complete its cycle, which
limits its frequency to no more than radians per period.
A motion of higher angular velocity will be mistaken for one of a velocity
that lies within the observable interval, and this is described as the problem of
17
D.S.G. POLLOCK: Statistical Fourier Analysis
0
0.25
0.5
0.75
1
0 2 2
A B
Figure 6. The function f(/2 + ), represented by the semi-continuous line, su-
perimposed upon the function f(/2). The sum of these functions is a 2-periodic
function f
H
(), which represents the spectral density of a subsampled process. The
segments A and B of f(/2 + ) that fall within the interval [, ] may be con-
strued as the eects of wrapping f(/2), dened over the interval [2, 2], around
a circle of circumference 2.
aliasing. The problem is illustrated by Figure 5 which shows that, by sampling
it at the integer points {t = 0, 1, 2, . . .}, a sinusoid with frequency of (11/8)
radians per period, which is in excess of the limiting Nyquist value of , will
be mistaken for one with a frequency of (5/8).
The extent to which aliasing is a problem depends upon the structure of
the particular time series under analysis. It will be suggested later that, for
many econometric time series, the problem does not arise, for the reason that
their (positive) frequencies are band-limited to a subinterval of the range [0, ].
However, this fact, which has its own problems, is not commonly recognised in
econometric time-series analysis.
To understand the statistical aspects of aliasing, we may consider a sta-
tionary stochastic process y(t) = {y
t
; t = 0, 1, 2, . . .} with an autocovariance
function () = {

; = 0, 1, 2, . . .}. The relationship between the spectral


density function and the autocovariance function is a follows:
() =
_

f()e
i
d f() =
1
2

e
i
. (45)
When alternate values are selected from the data sequence, the autocovariance
function is likewise subsampled, and there is
(2) =
_

f()e
i(2)
d
=
1
2
_
2
2
f(/2)e
i
d (46)
=
1
2
__

2
f(/2)e
i
d +
_

f(/2)e
i
d +
_
2

f(/2)e
it
d
_
.
18
D.S.G. POLLOCK: Statistical Fourier Analysis
Here, we have dened = 2, and we have used the change of variable tech-
nique to obtain the second equality.
Within the integrand of (46), there is exp{i} = exp{i( 2)}. Also,
f({/2} ) = f({/2} + ), by virtue of the 2-periodicity of the function.
Therefore,
_

2
f(/2)e
i
d =
_

0
f({/2} )e
i(2)
d
=
_

0
f({/2} + )e
i
,
(47)
where the rst equality is an identity and the second exploits the results above.
Likewise,
_
2

f(/2)e
i
d =
_
0

f({/2} + )e
i(+2)
d
=
_
0

f({/2} + )e
i
.
(48)
It follows that within the expression of (46), the rst integral may be translated
to the interval [0, ], whereas the third integral may be translated to the interval
[, 0]. After their translations, the rst and the third integrands can be
combined to form the segment of the function f( + /2) that falls in the
interval [, ]. The consequence is that
(2) =
1
2
_

{f(/2) + f( + /2)}e
i
d. (49)
It follows that
(2) f
H
() =
1
2
{f(/2) + f( + /2)}, (50)
where f
H
() = f
H
( + 2) is the spectral density function of the subsampled
process. Figure 6 shows the relationship between f(/2) and f( + /2).
The superimposition of the shifted frequency function f(+/2) upon the
function f(/2) corresponds to a process of aliasing. To visualise the process,
one can imagine a single cycle of the original function f() over the interval
[, ]. The dilated function f(/2) manifests a single cycle over the inter-
val [2, 2]. By wrapping this segment twice around the circle dened by
exp{i} with [, ], we obtain the aliased function {f(/2)+f(+/2)}.
the
Observe that, if f() = 0 for all || > /2, then the support of the dilated
function would be the interval [, ], and the rst and third integrals would
be missing from the nal expression of (46). In that case, there would be no
aliasing, and the function f
H
() would represent the spectrum of the original
process accurately.
In the case of a nite sample, the frequency-domain representation and the
time-domain representation are connected via the discrete Fourier transform.
19
D.S.G. POLLOCK: Statistical Fourier Analysis
We may assume that the size T of the original sample is an even number.
Then, the elements of the sequence {x
0
, x
2
, x
4
, . . . , x
T1
}, obtained by selecting
alternate points, are described by
x
2t
=
T1

j=0

j
exp{i
1
2tj}, where
1
=
2
T
=
(T/2)1

j=0
{
j
+
j+(T/2)
} exp{i
2
tj}, where
2
= 2
1
=
4
T
.
(51)
Here, the second equality follow from the fact that
exp{i
2
j} = exp{i
2
[j mod (T/2)]}.
To envisage the eect of equation (51), we can imagine wrapping the se-
quence
j
; j = 0, 1, . . . , T 1 twice around a circle of circumference T/2 so that
its elements fall on T/2 points at angles of
2j
= 4j/T; j = 1, 2, . . . , T/2 from
the horizontal.
Undoubtedly, the easiest way to recognise the eect of aliasing in the case
of a nite sample is to examine the eect upon the matrix representation of the
discrete Fourier transform. For conciseness, we adopt the expressions under
(33), and we may take the case were T = 8. Subsampling by a factor of 2 is
a matter of removing from the matrix

U alternate rows corresponding to the
odd-valued indices. This gives
_

_
x
0
x
2
x
4
x
6
_

_
=
_

_
1 1 1 1 1 1 1 1
1 W
6
8
W
4
8
W
2
8
1 W
6
8
W
4
8
W
2
8
1 W
4
8
1 W
4
8
1 W
4
8
1 W
4
8
1 W
2
8
W
4
8
W
6
8
1 W
2
8
W
4
8
W
6
8
_

_
_

7
_

_
=
_

_
1 1 1 1
1 W
3
4
W
2
4
W
4
1 W
2
4
1 W
2
4
1 W
4
W
2
4
W
3
4
_

_
_

0
+
4

1
+
5

2
+
6

3
+
7
_

_
.
(52)
We see that, in the nal vector on the RHS, the elements
0
,
1
,
2
,
3
are
conjoined with the elements
4
,
5
,
6
,
7
, which are associated with frequencies
that exceed the Nyquist rate.
5. Linear Filters in Time and Frequency
In order to avoid the eects of aliasing in the process of subsampling the data,
it is appropriate to apply an anti-aliasing lter to remove from the data those
frequency components that are in excess of the Nyquist frequency of the sub-
sampled date. In the case of downsampling by a factor of two, we should need
20
D.S.G. POLLOCK: Statistical Fourier Analysis
to remove, via a preliminary operation, those components of frequencies in
excess of /2.
A linear lter is dened by a set of coecients {
k
; k = p, . . . , q}, which
are applied to the data sequence y(t) via a process of linear convolution to give
a processed sequence
x(t) =
q

k=p

k
y(t k). (53)
On dening the z-transforms (z) =

k

k
z
k
, y(z) =

t
y
t
z
t
and x(z) =

t
x
t
z
t
, we may represent the ltering process by the following polynomial or
series equation:
x(z) = (z)y(z). (54)
Here, (z) is described as the transfer function of the lter.
To ensure that the sequence x(t) is bounded if y(t) is bounded, it is nec-
essary and sucient for the coecient sequence to be absolutely summable,
such that

k
|
k
| < . This is known by engineers as the bounded input
bounded output (BIBO) condition. The condition guarantees the existence of
the discrete-time Fourier transform of the sequence. However, a meaning must
be sought for this transform.
Consider, therefore, the mapping of a (doubly-innite) cosine sequence
y(t) = cos(t) through a lter dened by the coecients {
k
}. This produces
the output
x(t) =

k
cos([t k])
=

k
cos(k) cos(t) +

k
sin(k) sin(t)
= cos(t) + sin(t) = cos(t ),
(55)
where =

k

k
cos(k), =

k

k
sin(k), =

(
2
+
2
) and =
tan
1
(/). These results follow in view of the trigonometrical identity of (3).
The equation indicates that the lter serves to alter the amplitude of the
cosine signal by a factor of and to displace it in time by a delay of = /
periods, corresponding to an angular displacement of radians. A (positive)
phase delay is inevitable if the lter is operating in real time, with j 0, which
is to say that only the present and past values of y(t) are comprised by the
lter. However, if the lter is symmetric, with
k
=
k
, which necessitates
working o-line, then = 0 and, therefore, = 0.
Observe that the values of and are specic to the frequency of
the cosine signal. The discrete-time Fourier transform of the lter coecients
allows the gain and the phase eects to be assessed over the entire range of
frequencies in the interval [0, ]. The transform may be obtained by setting
z = exp{i} = cos() i sin() within (z), which constrains this argument
to lie on the unit circle in the complex plane. The resulting function
(exp{i}) =

k
cos(k) i

k
sin(k)
= () i()
(56)
21
D.S.G. POLLOCK: Statistical Fourier Analysis
is the frequency response function, which is, in general, a periodic complex-
valued function of with a period of 2. In the case of a symmetric lter, it
becomes a real-valued and even function, which is symmetric about = 0.
For a more concise notation, we may write () in place of (exp{i}).
In that case, the relationship between the coecients and their transform can
be conveyed by writing {
k
} ().
The eect of a linear lter on a stationary stochastic process, which is
not summable, can be represented in terms of the autocovariance generating
functions f
y
(z) and f
x
(z) of the sequence and its transform, respectively. Thus
f
x
(z) = |(z)|
2
f
y
(z), (57)
where |(z)|
2
= (z)(z
1
). Setting z = exp{i} gives the relationship
between the spectral density functions, which entails the squared gain |()|
2
=
|(e
i
)|
2
of the lter.
The Ideal Filter
Imagine that it is required to remove from a data sequence the components
of frequencies in excess of
d
within the interval [0, ]. This requirement may
arise from an anti-aliasing operation, preliminary to sub sampling, or it may
arise from a desire to isolate a component that lies in the frequency range
[0,
d
]. We shall discuss the motivation in more detail later.
The ideal lter for the purpose would be one that preserves all components
with (positive) frequencies in the subinterval [0,
d
) and nullies all those with
frequencies in the subinterval (
d
, ]. The specication of the frequency re-
sponse function of such a lter over the interval [, ] is
() =
_

_
1, if || (0,
d
),
1/2, if =
d
,
0, otherwise.
(58)
The coecients of the lter are given by the sinc function:

k
=
1
2
_

d

d
e
ik
d =
_

, if k = 0;
sin(
d
k)
k
, if k = 0 .
(59)
These coecients form a doubly-innite sequence and, in order to apply
such a lter to a data sample of size T, it is customary to truncate the sequence,
retaining only a limited number of its central elements (Figure 7).
The truncation gives rise to a partial-sum approximation of the ideal fre-
quency response that has undesirable characteristics. In particular, there is
a ripple eect whereby the gain of the lter uctuates within the pass band,
where it should be constant with a unit value, and within the stop band, where
it should be zero-valued. Within the stop band, there is a corresponding prob-
lem of leakage whereby the truncated lter transmits elements that ought to be
22
D.S.G. POLLOCK: Statistical Fourier Analysis
0
0.1
0.2
0.3
0.4
0.5
0
0.1
0 5 10 15 0 5 10 15
Figure 7. The central coecients of the Fourier transform of the frequency response
of an ideal low pass lter with a cut-o point at = /2. The sequence of coecients
extends indenitely in both directions.
0
0.25
0.5
0.75
1
1.25
0
0 /2 /2
Figure 8. The result of applying a 17-point rectangular window to the coecients
of an ideal lowpass lter with a cut-o point at = /2.
0
0.25
0.5
0.75
1
0
0 /2 /2
Figure 9. The result of applying a 17-point Blackman window to the coecients of
an ideal lowpass lter with a cut-o point at = /2.
23
D.S.G. POLLOCK: Statistical Fourier Analysis
blocked (Figure 8). This is nothing but Gibbs eect, which has been discussed
already in section 3 and illustrated in Figure 2.
The classical approach to these problems, which has been pursued by elec-
trical engineers, has been to modulate the truncated lter sequence with a
so-called window sequence, which applies a gradual taper to the higher-order
lter coecients. (A full account of this has been given by Pollock, 1999.)
The eect is to suppress the leakage that would otherwise occur in regions
of the stop band that are remote from the regions where the transitions occur
between stop band and pass band. The detriment of this approach is that it
exacerbates the extent of the leakage within the transition regions (Figure 9).
The implication here seems to be that, if we wish to make sharp cuts in
the frequency domain to separate closely adjacent spectral structures, then we
require a lter of many coecients and a very long data sequence to match it.
However, Figures 8 and 9 are based upon the discrete-time Fourier transform,
which entails an innite sequence in the time domain and a continuous function
in the frequency domain; and their appearances are liable be misleading.
In practice, one must work with a nite number of data points that trans-
form into an equal number of frequency ordinates. If one is prepared to work in
the frequency domain, then one can impose an ideal frequency selection upon
these ordinates by setting to zero those that lie within the stop band of the
ideal lter and conserving those that lie within the pass band. One can ignore
whatever values a continuous frequency response function might take in the in-
terstices between these frequency ordinates, because they have no implications
for the nite data sequence.
The same eects can be achieved by operations performed within the time
domain. These entail a circular lter, which is applied to the data by a pro-
cess of circular convolution that is akin to the process of circular correlation
described and illustrated in section 3. In the case of circular convolution, the
data sequence and the sequence of lter coecients are disposed around the
circle in opposite directionsclockwise and anti clockwise. The set of T coe-
cients of a circular lter are derived by applying the (inverse) discrete Fourier
transform to a sample of the ordinates of the ideal frequency response, taken
at the T Fourier frequencies
j
= 2j/T; j = 0, 1, . . . , T 1.
To understand the nature of the lter, consider evaluating the frequency
response of the ideal lter at one of the points
j
. The response, expressed in
terms of the doubly-innite coecient sequence of the ideal lter, is
(
j
) =

k=

k
W
jk
, (60)
where W
j
= exp{i
j
} = cos(
j
) i sin(
j
). However, since W
q
is a T-
periodic function, there is W q = W (q mod T), where the upward arrow
signies exponentiation. Therefore,
24
D.S.G. POLLOCK: Statistical Fourier Analysis
(
j
) =
_

q=

qT
_
+
_

q=

qT+1
_
W
j
+
+
_

q=

qT+T1
_
W
j(T1)
(61)
=

0
+

1
W
j
+ +

T1
W
j(T1)
, for j = 0, 1, 2, . . . , T 1.
Thus, the coecients

0
,

1
, . . . ,

T1
of the circular lter would be ob-
tained by wrapping the innite sequence of sinc function coecients of the ideal
lter around a circle of circumference T and adding the overlying coecients.
By applying a sampling process to the continuous frequency response function,
we have created an aliasing eect in the time domain. However, summing the
sinc function coecients is not a practical way of obtaining the coecients of
the nite sample lter. They should be obtained, instead, by transforming the
frequency-response sequence or else by one of the available analytic formulae.
In the case where
d
= d
1
= 2d/T, the coecients of the lter that
fulls the specication of (58) within a nite sample of T points are given by

d
(k) =
_

_
2d
T
, if k = 0
cos(
1
k/2) sin(d
1
k)
T sin(
1
k/2)
, for k = 1, . . . , [T/2],
(62)
where
1
= 2/T and where [T/2] is the integral part of T/2.
A more general specication for a bandpass lter is as follows:
() =
_

_
1, if || (
a
,
b
),
1/2, if =
a
,
b
,
0, for elsewhere in [, ),
(63)
where
a
= a
1
and
b
= b
1
. This can be fullled by subtracting one lowpass
lter from another to create

[a,b]
(k) =

b
(k)

a
(k)
=
cos(
1
k/2){sin(b
1
k) sin(a
1
k)}
T sin(
1
k/2)
= 2 cos(g
1
k)
cos(
1
k/2) sin(d
1
k)
T sin(
1
k/2)
.
(64)
Here, 2d = b a is the width of the pass band (measured in terms of a number
of sampled points) and g = (a + b)/2 is the index of its centre. These results
have been demonstrated by Pollock (2006).
The relationship between the time-domain and the frequency-domain l-
tering can be elucidated by considering the matrix representation of the lter,
25
D.S.G. POLLOCK: Statistical Fourier Analysis
0
0.25
0.5
0.75
1
0
0 /2 /2
Figure 10. The frequency response of the 16-point wrapped lter dened over the
interval [, ). The values at the Fourier frequencies are marked by circles and
dots.
0
0.25
0.5
0.75
1
1.25
0
0.25
0 /2 /2
Figure 11. The frequency response of the 17-point wrapped lter dened over the
interval [, ). The values at the Fourier frequencies are marked by circles.
which is obtained by replacing the generic argument W
j
= exp{i
j
} in (60)
and (61) by the matrix operator K. This gives

(K) =

=

U

(D)U. (65)
Replacing z in the z-transform y(z) of the data sequence by K gives the cir-
culant matrix Y =

Uy(D)U. Multiplying the data matrix by the lter matrix
gives
x(K) = X =

Y = {

(D)U}{

Uy(D)U}
=

U{

(D)y(D)}U =

U

(D)UY,
(66)
where the nal equality follows in view of the identity UY = y(D)U.
Here, there is a convolution in the time domain, represented by the product
X =

Y of two circulant matrices, and a modulation in the frequency domain,


represented by the product x(D) =

(D)y(D) of two diagonal matrices. Of


these matrices,

(D) comprises the zeros and units of the ideal frequency


response, whereas y(D) comprises the spectral ordinates of the data.
26
D.S.G. POLLOCK: Statistical Fourier Analysis
It should be emphasised that () and

(
j
) are distinct functions of
which, in general, the values coincide only at the roots of unity. Strictly speak-
ing, these are the only points for which the latter function is dened. Never-
theless, there may be some interest in discovering the values that

(
j
) would
deliver if its argument were free to travel around the unit circle.
Figure 10 is designed to satisfy such curiosity. Here, it can be seen that,
within the stop band, the function

() is zero-valued only at the Fourier


frequencies. Since the data are compounded from sinusoidal functions with
Fourier frequencies, this is enough to eliminate from the sample all elements
that fall within the stop band and to preserve all that fall within the pass band.
In Figure 10, the transition between the pass band and the stop band oc-
curs at a Fourier frequency. This feature is also appropriate to other lters that
we might design, that have a more gradual transition with a mid point that also
falls on a Fourier frequency. Figure 11 shows that it is possible, nevertheless, to
make an abrupt transition within the space between two adjacent frequencies.
An alternative way of generating the ltered sequence within the time
domain is to create a periodic extension of the data sequence via successive
replications of the data. Then, the ltered values can be generated via an
ordinary linear convolution of the wrapped lter coecients of
0
(z) with the
data. It can be seen that these products are synonymous with the product
of the innite sequence of sinc function coecients with an innite periodic
extension of the data.
A modulation in the frequency domain, which requires only T multipli-
cations, is a far less laborious operation than the corresponding circular con-
volution in the time domain, which requires T
2
multiplications and T(T 1)
additions. To be set against the relative eciency of the frequency-domain
operations is the need to carry the data to that domain via a Fourier trans-
form and to carry the products of the modulation back to the time domain via
an inverse transformation. However, given the eciency of the Fast Fourier
transform, the advantage remains with the frequency-domain operations. Also,
the time-domain lter may have to be formed by applying an inverse Fourier
transform to a set of ordinates sampled from the frequency response of the
lter.
The ideal lowpass lter has been implemented in a program that can be
downloaded from the following web address:
http://www.le.ac.uk/users/dsgp1/
Various devices are available within the program for dealing with trended
data sequences. A trend function can be interpolated through the data to
deliver a sequence of residual deviations to which the lter can be applied.
Thereafter, the ltered residuals and the trend can be recombined to produce
a trendcycle trajectory.
Alternatively, the trend can be eliminated from the data by a twofold ap-
plication of the dierence operator. The lter can be applied to the dierenced
data and the result can be reinated via the summation operator, which is the
inverse of the dierence operator.
The program also uses the method of frequency-domain ltering to create
27
D.S.G. POLLOCK: Statistical Fourier Analysis
an ideal bandstop lter that is able to eliminate the seasonal component of a
trended econometric data sequence.
WienerKolmogorov Theory
The classical WienerKolmogorov ltering theory concerns a doubly in-
nite data sequence y(t) = {y
t
; t = 0, 1, 2, . . .}, which is the sum of a signal
sequence (t) and a noise sequence (t):
y(t) = (t) + (t). (67)
The signal and noise components are assumed to be generated by mutually
independent stationary stochastic processes of zero mean. Therefore, the auto-
covariance generating function of the data, which may be represented by

y
(z) =

(z) +

(z), (68)
is the sum of the autocovariance generating functions of the components.
On this basis, a lter may be formulated that produces a minimum mean
square error estimate of the signal. Its transfer function is given by
(z) =

(z)

(z) +

(z)
. (69)
The innite-sample lter is a theoretical construct from which a practical im-
plementation must be derived.
A straightforward way of deriving a lter that is applicable to the nite
data sequence contained in a vector y = [y
0
, y
1
, . . . , y
T1
]

is to replace the ar-


gument z within the various autocovariance generating functions by the matrix
lag operator L
T
= [e
1
, . . . , e
T1
, 0] and to replace z
1
by F
T
= L

T
. Then, the
generating function
y
(z) =
0
+

=1

(z

+ z

) gives rise to matrix

y
=
0
I
T
+
T1

=1

(F

T
+ L

T
), (70)
where the limit of the summation is when = T 1, since L
q
T
= 0 for all q T.
This is an ordinary autocovariance matrix. The autocovariance matrices of the
components

and

are derived likewise, and then there is


y
=

.
Using the matrices in place of the autocovariance generating functions in
(68) gives rise to the following estimate x of the signal component :
x =

)
1
y (71)
It is necessary to preserve the order of the factors of this product since, unlike
the corresponding z-transforms,

and
1
y
= (

)
1
do not commute
in multiplication. Such lters have been discussed in greater detail by Pollock
(2000), (2001), where the matter of trended data sequences is also dealt with.
28
D.S.G. POLLOCK: Statistical Fourier Analysis
For an alternative nite-sample representation of the lter, we may replace
z in (68) and (69) by K
T
=

UDU and z
1
by K

T
= K
T1
T
= UD

U. Then, the
generating functions are replaced by the following matrices of circular autoco-
variances:

(K) =

U

(D)U,

(K) =

U

(D)U
and

y
=

=
y
(K) =

U
y
(D)U.
(72)
Putting these matrices in place of the autocovariance matrices of (71) gives rise
to the following estimate of the signal:
x =

U

(D){

(D) +

(D)}
1
Uy =

UJ

Uy. (73)
The estimate may be obtained in the following manner. First, a Fourier
transform is applied to the data vector y to give Uy. Then, the elements of the
transformed vector are multiplied by those of the diagonal weighting matrix
J

(D){

(D) +

(D)}
1
. Finally, the products are carried back to the
time domain by the inverse Fourier transform.
The same method is indicated by equation (66), which represent the ap-
plication of the ideal lter. When it is post multiplied by e
0
, the latter delivers
the expression x =

U

(D)Uy, which can also be given the foregoing interpre-


tation.
6. The Processes Underlying the Data
In this section, we shall discuss the relationship between the sampled data and
the processes that generate it, which are presumed to operate in continuous
time.
If the continuous process contains no sinusoidal elements with frequencies
in excess of the Nyquist value of radians, then a perfect reconstruction of the
underlying function can be obtained from its sampled values by using the sinc
function sin(t)/t of Figure 3 as an interpolator. To eect the interpolation,
each of the data points is associated with a sinc function scaled by the sampled
value. These functions, which are set at unit displacements along the real line,
are added to form a continuous envelope.
The eect of a sinc function at any integer point t is to disperse the mass
of the sampled value over the adjacent points of the real line. The function has
a value of unity at t and a value of zero on all other integer points. Therefore,
in the process of interpolation, it adds nothing to the values at the adjacent
integer points. However, it does attribute values to the non-integer points that
lie in the interstices, thereby producing a continuous function from discrete
data. Moreover, the data can be recovered precisely by sampling the continuous
function at the integers.
Only if the sampled data is free of aliasing, will this process of interpolation
succeed in recovering the original continuous function. The question of whether
aliasing has occurred is an empirical one that relates to the specic nature of the
process underlying the data. However, the models of continuous-time processes
29
D.S.G. POLLOCK: Statistical Fourier Analysis
that are commonly adopted by econometricians and other statisticians suggest
that severe aliasing is bound to occur.
The discrete-time white-noise process (t) that is the primum mobile of
linear stochastic processes of the ARMA and ARIMA varieties is commonly
thought to originate from sampling a Wiener process W(t) dened by the fol-
lowing conditions:
(a) W(0) = 0,
(b) E{W(t)} = 0, for all t,
(c) W(t) is normal,
(d) dW(s), dW(t) for all t = s are independent stationary increments,
(e) (t) =
1
h
{W(t + h) W(t)} for some h > 0.
The increments dW(s), dW(t) are impulses that have a uniform power spectrum
distributed over the entire real line, with innitesimal power in any nite inter-
val. Sampling W(t) at regular intervals entails a process of aliasing whereby the
spectral power of the cumulated increments gives rise to a uniform spectrum
of nite power over the interval [, ].
An alternative approach is to imagine that the discrete-time white noise
is derived by sampling a continuous-time process that is inherently limited in
frequency to the interval [, ]. The sinc function of Figure 3, which is limited
in frequency to this interval, can also be construed as a wave packet centered
on time t = 0. A continuous stochastic function that is band-limited by the
Nyquist frequency of will be generated by the superposition of such functions
arriving at regular or irregular intervals of time and having random amplitudes.
No aliasing would occur in the sampling of such a process.
(The wave packet is a concept from quantum mechanicssee, for example
Dirac (1958) and Schi (1981)that has become familiar to statisticians via
wavelet analysissee, for example, Daubechies (1992) and Percival and Walden
(2000).)
Band-Limited Processes
It is implausible to assume that a Wiener process could be the primum
mobile of a stochastic process that is band limited within a sub interval of the
(positive) frequency range [0, ]. Such processes appear to be quite common
amongst the components of econometric time-series.
Figure 12 displays a sequence of the logarithms of the quarterly series
of U.K. Gross Domestic Product (GDP) over the period from 1970 to 2005.
Interpolated through this sequence is a quadratic trend, which represents the
growth path of the economy.
The deviations from this growth path are a combination of the low-
frequency business cycle with the high-frequency uctuations that are due to
the seasonal nature of economic activity. These deviations are represented in
Figure 13, which also shows an interpolated continuous function that is designed
to represent the business cycle.
30
D.S.G. POLLOCK: Statistical Fourier Analysis
11.25
11.5
11.75
12
12.25
12.5
12.75
0 25 50 75 100 125
Figure 12. The quarterly sequence of the logarithms of the GDP in the U.K. for the
years 1970 to 2005, inclusive, together with a quadratic trend interpolated by least
squares regression.
0
0.025
0.05
0.075
0.1
0.025
0.05
0.075
0 25 50 75 100 125
Figure 13. The residual sequence from tting a quadratic trend to the income data
of Figure 12. The interpolated line represents the business cycle.
0
0.002
0.004
0.006
0 /4 /2 3/4
Figure 14. The periodogram of the residuals obtained by tting a quadratic trend
through the logarithmic sequence of U.K. GDP. The parametric spectrum of a tted
AR(2) model is superimposed upon the periodogram.
31
D.S.G. POLLOCK: Statistical Fourier Analysis
0
0.05
0.1
0.15
0
0 5 10 15 20 0 5 10 15 20
Figure 15. The sinc function wave-packet
3
(t) = sin(t/8)/t comprising fre-
quencies in the interval [0, /8].
The periodogram of the deviations is shown in Figure 14. This gives a clear
indication of the separability of the business cycle and the seasonal uctuations.
The spectral structure extending from zero frequency up to /8 belongs to the
business cycle. The prominent spikes located at the frequency /2 and at the
limiting Nyquist frequency of are the property of the seasonal uctuations.
Elsewhere in the periodogram, there are wide dead spaces, which are punctu-
ated by the spectral traces of minor elements of noise.
The slowly varying continuous function interpolated through the deviations
of Figure 13 has been created by combining a set of sine and cosine functions of
increasing frequencies in the manner of equation (1), but with the summation
extending no further than the limiting frequency of the business cycle, which
is /8.
A process that is band limited to the frequency interval [0, /8] would be
derived from the random arrival of wave packets of random amplitudes, likewise
band limited in frequency. The sinc function is but one example such a packet;
and it has the disadvantage of being widely dispersed in time. Its advantage is
that it provides an orthogonal basis for the set of band-limited functions. The
basis set is
{
3
(t 8k); k = 0, 1, 2, . . .} with
3
(t) =
sin(t/8)
t
, (74)
wherein the functions are displaced, from one to the next, by 8 sample points.
(Here the subscript on
3
is in reference to the fact that 8 = 2
3
and that the
interval [0, /8] arises from the 3rd dyadic subdivision of [0, ].) The function at
t = 0 is represented in Figure 15. Whereas a continuous band-limited stochastic
function would be constituted by adding the envelopes of the succession of
continuous sinc functions, its sampled version would be obtained by adding the
sampled ordinates that are also represented in Figure 15.
A further point concerns the autocovariances of a band-limited process.
The autocovariance function of an ordinary linear stochastic process of the
ARMA variety is positive denite, which corresponds to the nonzero property
32
D.S.G. POLLOCK: Statistical Fourier Analysis
of the spectral density function. Likewise, the matrix of autocovariances is pos-
itive denite, which implies that all of its eignenvalues are positive. Therefore,
it is not possible to represent a band-limited process by such means.
However, it is straightforward to represent the autocovariances of a band-
limited process via a circulant matrix of the form
(K) = =

U(D)U. (75)
Here, the diagonal matrix (D) represents the spectral density of the process,
sampled at T frequency points. The band-limited nature of the process will
be reected by the zero-valued elements of the diagonal of (D) and by the
singularity of .
Fitting ARMA Models to the Business Cycle
One might wish to t a parametric model to the detrended data of Figure
13. For this purpose, the multiplicative seasonal ARIMA model of Box and
Jenkins (1976) seems, at rst, to be a natural choice. Such a model embodies
a seasonal dierence operator of the form (1 z
4
) = (1 z)(1 +z +z
2
+z
3
), of
which the ordinary dierence operator (1 z) is a factor. However, the data of
Figure 13, have already been detrended. Therefore, the presence or the redun-
dant dierence operator renders the model inappropriate. A more appropriate
recourse is to remove the seasonal component from the data before tting an
autoregressive model aimed at representing the business cycle component.
The seasonal component can be removed from the data of Figure 13 by
nullifying the spectral ordinates at the seasonal frequencies and at some ad-
jacent frequencies. This operation, which can be performed by the program
mentioned in Section 5, is represented in Figure 14 by the highlighted bands
that correspond to the stop bands of the deseasonalising lter. Also represented
in Figure 14 is the parametric spectrum of a second-order autoregressive AR(2)
model that has been tted to the seasonally adjusted data,
One would expect to be able to characterise the low-frequency dynamics
of the business cycle component via a pair of conjugate complex autoregressive
roots to be found within the estimated model. The arguments of such roots
should characterise the angular velocity of the cycle; and their moduli should
represent the degree of its damping. However, it is found that the roots in
question are real-valued. It has been reported by Pagan (1997), amongst others,
that this is the almost invariable outcome of such exercises.
We might wish, nevertheless, to use an ARMA model to characterise the
business cycle component that is supported on the interval [0, /8] and which
is represented by the smooth continuous function that has been interpolated
through the data of Figure 13. In general, when we wish to t an ARMA
model to a component of the data that is supported on a frequency interval
[0,
d
], where
d
< , it is necessary to map the component onto the full
frequency interval [0, ]. For this purpose, it may be necessary to synthesise a
continuous function from the Fourier ordinates within the subinterval [0,
d
],
as has been done in Figure 13 for the case where
d
= /8. The function can
33
D.S.G. POLLOCK: Statistical Fourier Analysis
then be resampled at the lesser rate of one observation in every /
d
periods
to produce data that are appropriate to the purpose.
In the case of the business cycle component of Figure 13, which has a
spectral signature that is supported on the interval [0, /8], it is sucient to
take one in eight of the data points of the anti-aliased sequence that can been
obtained by subjecting the original data to an ideal lowpass lter with a cut-o
at /8. The periodogram of the resulting subsampled sequence, which has the
frequency range of [0, ], will have a prole identical to that of the spectral of
structure that occupies the interval [0, /8] in Figure 14.
The AR(2) model that ts these data has a pair of conjugate complex roots
with an argument of = 1.35 radians (77

), which corresponds to a cyclical


period of 37 quarters, and a modulus of = 0.695. However, the evidence
of the interpolated function of Figure 13 and of the periodogram of Figure 14
suggests that the estimated period is an average of the durations longer major
cycle and a shorter minor cycle. However, with only 18 points available in the
subsampled data, it seems that there is insucient information to justify the
presentation of the AR(4) model that reects these putative cycles in two pairs
of conjugate complex roots.
The Recognition of the Frequency Limitation
Some reasons should be given for why the prevalence of band limited pro-
cesses has been overlooked by econometricians. One evident reason is that, in
the past, econometricians have rarely examined the periodograms of their data
sequences.
However, the band limited nature of a process is liable to be obscured
whenever there are signicant disjunctions in the periodic extension of the data
at the points where the end of one replication of the sample joins the beginning
of the next. Such disjunctions will give rise to a slew of spectral power that
will serve to conceal the underlying structure of a band-limited process.
The presence of a barely perceptible trend in the data will be enough to
cause such a disjunction. The periodic extension of a nite linearly trending se-
quence will have the appearance of the serrated edge of a carpenters saw. The
spectrum of a saw tooth function has the one-over-f prole of a rectangular
hyperbola, which will mask the ner spectral details. (See Hamming (1989),
for example.) Granger (1966) has described such spectra as typical of econo-
metric time series. His characterisation has greatly inuenced the perceptions
of econometricians. It has suggested to many of them that the spectral analysis
of economic data is liable to be dicult if not fruitless.
A third reason for the failure to perceive the band-limited nature of some
of these process may be that econometricians have been unduly inuenced by
the prevalence of the Wiener process as a model of the forcing functions of
continuous-time stochastic processes. A Wiener process, which has a frequency
band of innite width, is the antithesis of a band-limited process.
Within the computer program that is an adjunct of this paper, there are
several provisions for ensuring that the data can be wrapped seamlessly around
the circumference of a circle in a manner that avoids disjunctions in the periodic
extension.
34
D.S.G. POLLOCK: Statistical Fourier Analysis
If the data are to be detrended by the interpolation of a polynomial func-
tion, then extra weight can be given to the points at the beginning and the
end of the sample to ensure that the polynomial passes through their midst.
This may have the desired eect upon the residual sequence, which will be used
thereafter in the analysis.
If the data have already been detrended, then it may be desirable to add
tapered extrapolations at both ends which converge on a common asymptote.
The extrapolations are liable to be discarded after the data have been ltered.
Alternatively, if the detrended data display prominent seasonal uctuations of
a pattern that is changing over the course of the sample, then a segment may
be interpolated into the circular data in which the pattern of the nal year of
the sample gradually morphs into the pattern of the rst year.
References
Baher, H., (1990), Analog and Digital Signal Processing, John Wiley and Sons,
Chichester.
Box, G.E.P., and G.M. Jenkins, (1976), Time Series Analysis: Forecasting and
Control: Revised Edition, Holden-Day, San Francisco.
Brigham, E.O., (1988), The Fast Fourier Transform and its Applications,
Prentice-Hall, Englewood Clis, New Jersey.
Brillinger, D.R., (1975), Time Series: Data Analysis and Theory, Holt, Rine-
hart and Winston, News York.
Brockwell, P.J., and R.A. Davis, (1987), Time Series: Theory and Methods,
Springer-Verlag, New York.
Carslaw, H.S., (1925), A Historical Note on Gibbs Phenomenon in Fourier
Series and Integrals, Bulletin of the American Mathematical Society, 31, 420
424.
Carslaw, H.S., (1930), Introduction to the Theory of Fouriers Series and Inte-
grals, Third Edition, Macmillan, London.
Catalan, E., (1846), Recherches sur les Determinants, Bulletin de
lAcademie Royale de Belgique, 13, 534555.
Churchill, R.V., and J.W. Brown, (1978), Fourier Series and Boundary Value
Problems, Third edition, McGraw-Hill, New York.
Daubechies, I., (1992), Ten Lectures on Wavelets, Society for Industrial and
Applied Mathematics, Philadelphia.
Davis, P.J., (1979), Circulant Matrices, John Wiley and Sons, New York.
Dirac, P.A.M., (1958), The Principles of Quantum Mechanics, Fourth Edition,
Oxford University Press, Oxford.
Fishman, G.S., (1969), Spectral Methods in Econometrics, Harvard University
Press, Cambridge, Massachusetts.
35
D.S.G. POLLOCK: Statistical Fourier Analysis
Fuller, W.A., (1976), Introduction to Statistical Time Series, John Wiley and
Sons, New York.
Granger, C.W.J, (1966), The Typical Spectral Shape of an Economic Variable,
Econometrica, 34, 150161.
Granger, C.W.J., and M. Hatanaka, (1964), Spectral Analysis of Economic
Time Series, Princeton University Press, Princeton, New Jersey.
Gray, R.M., (2002), Toeplitz and Circulant Matrices: A Review, Information
Systems Laboratory, Department of Electrical Engineering, Stanford Univer-
sity, Stanford, California, [email protected].
Grenander, U., and M. Rosenblatt, (1957), Statistical Analysis of Stationary
Time Series, John Wiley and Sons, New York.
Hamilton, J.D., (1994), Time Series Analysis, Princeton University Press,
Princeton, New Jersey.
Hamming, R.W., (1989), Digital Filters: Third Edition, Prentice-Hall, Engle-
wood Clis, New Jersey.
Hannan, E.J., (1960), Time-series Analysis, Methuen, London.
Jenkins, G.M., and D.G. Watts, (1968), Spectral Analysis and its Applications,
Holden-Day, Oakland, Calif.
Kolmogorov, A.N., (1941a), Interpolation und Extrapolation von station aren
zualigen Folgen, Bulletin de lacademie des sciences de U.S.S.R., Ser. Math.,
5, 314.
Kolmogorov, A.N., (1941b) Stationary Sequences in a Hilbert Space, (in Rus-
sian), Bulletin of the Moscow State University, 2, 140.
Koopmans, L., (1974), The Spectral Analysis of Time Series, Academic Press,
NewYork.
Kreider, D.L., R.G. Kuller, D.R. Osterberg and F.W. Perkins, (1966), An In-
troduction to Linear Analysis, Addison Wesley Publishing Co., Reading Mass.
Muir, T., (1906), The Theory of Circulants in the Historical Order of Develop-
ment up to 1860, Proceedings of the Royal Society of Edinburgh, 26, 390398.
Muir, T., (1911), The Theory of Circulants from 1861 to 1880, Proceedings of
the Royal Society of Edinburgh, 32, 136149.
Muir, T., (1915), The Theory of Circulants from 1880 to 1900, Proceedings of
the Royal Society of Edinburgh, 36, 151179.
Muir, T., (1923), The Theory of Circulants from 1900 to 1920, Proceedings of
the Royal Society of Edinburgh, 44, 218241.
Nerlove, M., D.M. Grether and J.L. Carvalho, (1979), Analysis of Economic
Time Series: a Synthesis, Academic Press, New York.
Pagan, A., (1997), Towards an Understanding of Some Business Cycle Charac-
teristics, The Australian Economic Review, 30, 115.
36
D.S.G. POLLOCK: Statistical Fourier Analysis
Percival, D.B., and A.T. Walden, (2000), Wavelet Methods for Time Series
Analysis, Cambridge University Press, Cambridge.
Pollock, D.S.G., (1999), A Handbook of Time-Series Analysis, Signal Processing
and Dynamics, Academic Press, London.
Pollock, D.S.G., (2000), Trend Estimation and De-Trending via Rational
Square-Wave Filters, Journal of Econometrics, 99, 317334.
Pollock, D.S.G., (2001), Filters for Short Non-stationary Sequences, Journal of
Forecasting, 20, 341355.
Pollock, D.S.G., (2006), The Realisation of Finite-Sample Frequency-Selective
Filters. Discussion paper of the Department of Economics, Queen Mary Col-
lege, University of London.
Priestley, M.B., (1981), Spectral Analysis and Time Series, Academic Press,
London.
Rosenblatt, M., (1985), Stationary Sequences and Random Fields, Birkh` auser,
Basel.
Rozanov, Y.A., (1967), Stationary Random Processes, Holden-Day, San Fran-
cisco.
Schi, L.I., (1981), Quantum Mechanics, McGraw-Hill Book Co., New York.
Stade, E., (2005), Fourier Analysis, John Wiley and Sons, Hoboken, New Jer-
sey.
Titchmarsh, E.C., (1937), Introduction to the Theory of Fourier Integrals,
Clarendon Press, Oxford.
Wiener, N., (1930), Generalised Harmonic Analysis, Acta Mathematica, 55,
117258.
Wiener, N., (1941), Extrapolation, Interpolation and Smoothing of Stationary
Time Series. Report on the Services Research Project DIC-6037. Published in
book form in 1949 by MIT Technology Press and John Wiley and Sons, New
York.
Yaglom, A.M., (1962), An Introduction to the Theory of Stationary Random
Functions, Prentice-Hall, Englewood Clis, New Jersey.
Appendix: The Spectral Representation of a Stationary Process
Consider a sample {y
0
, y
1
, . . . , y
T1
} of T observations on a stationary stochas-
tic process. These are amenable to a discrete Fourier transform, whereby they
can be expressed as
y
t
=
n

j=0
_

j
cos(
j
t) +
j
sin(
j
t)
_
, (A.1)
where
j
= 2j/T and n is the integral part of T/2.
37
D.S.G. POLLOCK: Statistical Fourier Analysis
By allowing n to tend to innity, it is possible to express a sequence of
indenite length in terms of a sum of sine and cosine functions. However, in
the limit as n , the coecients
j
,
j
tend to vanish; and, therefore, an
alternative representation in terms of dierentials is called for.
By writing
j
= dA(
j
),
j
= dB(
j
), where A(), B() are step func-
tions with discontinuities at the points {
j
; j = 0, . . . , n}, the expression (A.1)
can be rendered as
y
t
=

j
_
cos(
j
t)dA(
j
) + sin(
j
t)dB(
j
)
_
. (A.2)
In the limit, as n , the summation is replaced by an integral to give the
expression
y(t) =
_

0
_
cos(t)dA() + sin(t)dB()
_
. (A.3)
Here, cos(t) and sin(t), and therefore y(t), may be regarded as innite se-
quences dened over the entire set of positive and negative integers.
Since A() and B() are discontinuous functions for which no derivatives
exist, one must avoid using ()d and ()d in place of dA() and dB().
Moreover, the integral in (A.3) is a FourierStieltjes integral.
In order to derive a statistical theory for the process that generates y(t),
one must make some assumptions concerning the functions A() and B().
So far, the sequence y(t) has been interpreted as a realisation of a stochastic
process. If y(t) is regarded as the stochastic process itself, then the functions
A(), B() must, likewise, be regarded as stochastic processes dened over
the interval [0, ]. A single realisation of these processes now corresponds to a
single realisation of the process y(t).
The rst assumption to be made is that the functions dA() and dB()
represent a pair of stochastic processes of zero mean, which are indexed on the
continuous parameter . Thus
E
_
dA()
_
= E
_
dB()
_
= 0. (A.4)
The second and third assumptions are that the two processes are mutu-
ally uncorrelated and that non-overlapping increments within each process are
uncorrelated. Thus
E
_
dA()dB()
_
= 0 for all , ,
E
_
dA()dA()
_
= 0 if = ,
E
_
dB()dB()
_
= 0 if = .
(A.5)
The nal assumption is that the variance of the increments is given by
V
_
dA()
_
= V
_
dB()
_
= 2dF()
= 2f()d.
(A.6)
38
D.S.G. POLLOCK: Statistical Fourier Analysis
It can be seen that, unlike A() and B(), F() is a continuous dierentiable
function. The function F() and its derivative f() are the spectral distribu-
tion function and the spectral density function, respectively.
In order to express (A.3) in terms of complex exponentials, we may dene
a pair of conjugate complex stochastic processes:
dZ() =
1
2
_
dA() idB()
_
,
dZ

() =
1
2
_
dA() + idB()
_
.
(A.7)
Also, the domain of the functions A(), B() may be extended from [0, ] to
[, ] by regarding A() as an even function such that A() = A() and by
regarding B() as an odd function such that B() = B(). Then there is
dZ

() = dZ(). (A.8)
From the conditions of (A.5), it follows that
E
_
dZ()dZ

()
_
= 0 if = ,
E{dZ()dZ

()} = f()d.
(A.9)
These results may be used to reexpress (A.3) as
y(t) =
_

0
_
(e
it
+ e
it
)
2
dA() i
(e
it
e
it
)
2
dB()
_
=
_

0
_
e
it
{dA() idB()}
2
+ e
it
{dA() + idB()}
2
_
=
_

0
_
e
it
dZ() + e
it
dZ

()
_
.
(A.10)
When the integral is extended over the range [, ], this becomes
y(t) =
_

e
it
dZ(), (A.11)
which is commonly described as the spectral representation of the process y(t).
39

You might also like