arXiv:1710.07792v2 [math.ST] 17 Nov 2017

Cointegrated Density-Valued Linear Processes
Won-Ki Seo
Department of Economics, University of California, San Diego
November 20, 2017

Abstract
In data rich environments we may sometimes deal with time series
that are probability density-function valued, such as observations of
cross-sectional income distributions over time. To apply the methods of functional time series analysis to such observations, we should
first embed them in a linear space in which the essential properties
of densities are preserved under addition and scalar multiplication.
Bayes Hilbert spaces provide one way to achieve this embedding. In
this paper we investigate the use of Bayes Hilbert spaces to model
cointegrated density-valued linear processes. We develop an I(1) representation theory for cointegrated linear processes in a Bayes Hilbert
space, and adapt existing statistical procedures for estimating the corresponding attractor space to a Bayes Hilbert space setting. We revisit
empirical applications involving earnings and wage densities to illustrate the utility of our approach.
I thank Brendan K. Beare for a lot of discussion at all stages of this paper. I am
also grateful to Dakyung Seong and Lam Nguyen for useful discussion on section 6 of this
paper.

1

1

Introduction

While the subject of time series analysis has traditionally dealt with time
series taking values in finite dimensional Euclidean space, a recent literature
on functional time series analysis deals with time series taking values in an
infinite dimensional Banach or Hilbert space. Each observation of such a time
series is a functional object; for example, it could be a square-integrable function, continuous function, or probability density function. Bosq (2000) gives
a rigorous theoretical treatment of linear processes in Banach and Hilbert
spaces. Horváth and Kokoszka (2012) discuss statistical aspects of functional
data and time series analysis, and provide various empirical applications.
Granger (1981) introduced the notion of cointegration as a way to model
long-run equilibrium relationships between real-valued time series. A recent
paper by Chang et al. (2016) is the first to consider the possibility of cointegration in a functional time series setting. The authors consider a time
series of probability densities taking values in the space of square-integrable
real functions on a compact interval K, denoted by L2 (K), and provide a
notion of cointegration adapted to this space. They develop associated statistical methods based on functional principal component analysis (FPCA),
and provide empirical applications to time series of earnings and stock return
densities.
Beare (2017) notes that a technical complication arises in the framework
developed by Chang et al. (2016): the nonnegativity property of probability
densities is incompatible with the type of nonstationarity exhibited by integrated time series, and this incompatibility cannot be resolved by a simple
demeaning of densities. Nontrivial examples of cointegrated density-valued
processes in L2 (K) therefore do not exist. Even in models of stationary probability density-valued time series, Petersen and Müller (2016) have observed
that it is generally inadvisable to treat such time series as taking values in
the subset of L2 (K) consisting of probability densities, as this subset does
not form a linear subspace of L2 (K). It is clear that an arbitrary linear
combination of densities is not a proper density; only a convex combination
of densities is a proper one.
The primary purpose of this paper is to show that the technical complications just discussed can be resolved by viewing probability density-valued
observations as elements of a Bayes Hilbert space. Such spaces were introduced by Egozcue et al. (2006) and developed further by van den Boogaart
et al. (2014). They are constructed in such a way that different elements
2

of the space correspond to different probability densities, and the essential
properties of probability densities are preserved under addition and scalar
multiplication. The notion of cointegration may easily be adapted to this
space as in Beare et al. (2017), who extend the framework developed by
Chang et al. (2016) to arbitrary complex separable Hilbert spaces.
A secondary purpose of the paper is to study the behavior of densityvalued autoregressive processes of order p (AR(p) processes) taking values
in a Bayes Hilbert space. We provide conditions under which an AR(p)
law of motion has a stationary solution or I(d) solution for some positive
integer d, and a necessary and sufficient condition for such an I(d) solution
to in fact be I(1). These results are closely related to the so-called GrangerJohansen representation theorem and its generalization to a possibly infinite
dimensional Hilbert space setting by Beare et al. (2017) and Beare and Seo
(2017).
A third contribution of the paper is the provision of statistical methods
to estimate the attractor space (to be defined later) for an I(1) process in a
Bayes Hilbert space. These methods are based on the functional unit root
test proposed by Chang et al. (2016). We illustrate their usefulness with
empirical applications to time series of cross-sectional densities of individual
earnings and wages.
The remainder of the paper is organized as follows. In Section 2, we review some background material on Bayes Hilbert spaces and other essential
mathematical concepts. In Section 3, we explain how Bayes Hilbert spaces
can provide a useful setting for models of cointegrated density-valued time
series. Results on I(d) representations of cointegrated AR(p) processes in
Bayes Hilbert space are provided in Section 4. Statistical tools for studying cointegrated density-valued time series are introduced in Section 5, and
illustrated with empirical applications in Section 6.

2

Preliminaries

Here we briefly review essential background for the study of cointegrated
density-valued linear processes, and fix standard notation and terminology.

3

2.1

Bayes Hilbert spaces

Bayes Hilbert spaces provide the setting for our treatment of cointegrated
density-valued linear processes. The discussion provided here will omit some
details to conserve space. The reader is referred to Egozcue et al. (2006) and
van den Boogaart et al. (2014) for a rigorous introduction to the subject.
Let (M, A) be a measurable space and P be the set of σ-finite positive
real-valued measures defined on it. For any measure λ, called a reference
measure, define
M(λ) := {ν ∈ P : ν(A) = 0 ⇔ λ(A) = 0,

∀A ∈ A}.

That is, M(λ) is the collection of σ-finite positive measures absolutely continuous with respect to λ. Due to the Radon-Nikodym theorem, any element
ν in M(λ) may be identified with the λ-density dν/dλ, so hereafter we always
regard an element in M(λ) as its λ-density.
We define an equivalence relation ' on M(λ) as follows. For λ-densities
f, g ∈ M(λ), the equivalence f ' g holds if there is some a > 0 such that
Z
Z
gdλ for all A ∈ A.
f dλ = a ·
A

A

The above says that if two measures ν1 and ν2 are proportional in the sense
that ν1 (A) = a · ν2 (A) for all A ∈ A, their λ-densities are equivalent under
'. The collection of '-equivalence classes of λ-densities, denoted by B(λ),
are the elements of a Bayes space associated with λ. We define two vector
operations (⊕, ), called perturbation and powering, as follows. For a ∈ R
and f, g ∈ B(λ), define
Z
f ⊕g '

f gdλ,
Z

a

f'

f a dλ.

Negative perturbation is defined as f g = f ⊕ (−1 g). B(λ) equipped
with the two vector operations (⊕, ) is a vector space.
We hereafter assume that λ is a finite measure and let B 2 (λ) be the
collection of λ-densities that have square-integrable logarithms, i.e.


Z
2
2
B (λ) := f ∈ B(λ) :
|log f | dλ < ∞ .

4

It turns out that B 2 (λ) is a vector subspace of B(λ). Define the inner product
h·, ·iB 2 (λ) in the following way: for f, g ∈ B 2 (λ),


Z
Z
Z 
log f − log f dλ
log g − log f dλ dλ.
hf, giB 2 (λ) =
The following result is obtained in van den Boogaart et al. (2014).
Lemma 2.1. B 2 (λ) equipped with (⊕, ) and h·, ·iB 2 (λ) is a separable Hilbert
space.
Define L2 (λ) as the set of measurable square-integrable functions f : Ω →
R. An element of L2 (λ) is regarded as a class of square-integrable functions
that are λ-almost surely equivalent. L2 (λ) is assumed to be
R equipped with
the inner product h·, ·iL2 (λ) defined as follows: hf˜, g̃iL2 (λ) = f˜g̃dλ for f˜, g̃ ∈
L2 (λ). Define


Z
2
2
L (λ) = f˜ ∈ L (λ) : f˜ = 0 ,
which is a closed subspace of L2 (λ). One useful fact for the study of densityvalued functional objects is that there exists an isometric isomorphism between B 2 (λ) and L2 (λ).
Lemma 2.2. The following hold:
(a) B 2 (λ) is isometrically isomorphic to L2 (λ).
(b) clr : B 2 (λ) → L2 (λ) defined as
Z
clr(f ) = log f −

log f dλ,

f ∈ B 2 (λ),

is an isometrical isomorphism. Moreover, clr−1 is given by
clr−1 (f˜) ' exp(f˜),

f˜ ∈ L2 (λ).

See van den Boogaart et al. (2014) for more details. The map clr is called
the centered log-ratio (clr) transformation. Since clr is a unitary linear map,
the following are naturally implied:
clr(f ⊕ g) = clr(f ) + clr(g),
hf, giB 2 (λ)

clr(a f ) = a · clr(f ),
Z
= hclr(f ), clr(g)iL2 (λ) = clr(f ) clr(g)dλ.
5

Let λ be the uniform measure defined on a compact interval K. The Bayes
Hilbert space associated with λ, B 2 (λ), was considered in Egozcue et al.
(2006), and it may be understood as the space of probability densities with
support K. Note that B 2 (λ) is the space of equivalence classes of positive
functions defined on K with square-integrable logarithms with respect to
the
R uniform measure λ. Within an equivalence class, the integral constraint
f dλ = 1 singles out the representative element as a probability density
K
function. Then for f, g ∈ B 2 (λ), two operations ⊕ and are be defined as
f (x)g(x)
,
f
(y)g(y)dy
K

f ⊕ g(x) = R

a

f (x) = R

f (x)a
.
a dy
f
(y)
K

If we denote the indicator function on K by 1K , then the '-equivalence
class associated with 1K is the neutral element of perturbation. Naturally,
1 ∈ R is the neutral element of powering. Perturbation and powering are
constructed such that the results of f ⊕ g and a f are always densities even
if one or both of f and g are not proper densities. If one regards an element
in B 2 (λ) as a probability density, perturbation can be interpreted as a Bayes
updating. This is why B 2 (λ) is called a Bayes (Hilbert) space in the previous
literature.
The fact that the resulting density from any arbitrary linear combination
of two elements of B 2 (λ) is a proper density is crucial to the approach taken
in this paper. In fact, the point of Beare (2017)’s comments on Chang et al.
(2016) is rooted in the fact that linear combinations of probability densities
in L2 (K) are not necessarily probability densities. It will be shown in the
subsequent sections that the Bayes Hilbert space is not just a good candidate on which to define (cointegrated) density-valued processes, but also its
Hilbert space structure makes it possible to use well-developed statistical
techniques in Hilbert spaces without further technicality.
We will consider a time series of densities with common support K.
Therefore, λ will be assumed to be the uniform measure on K and densityvalued linear processes will be defined in the associated Bayes Hilbert space
B 2 (λ). Allowing the reference measure to be another measure with compact
support is trivial. However, generalization to a reference measure with unbounded support is nontrivial and we will not consider this case. When the
reference measure has unbounded support, the space of probability measures
as a subset of B 2 (λ) is not a subspace, as the following example demonstrates.

6

Example 2.1. Let λ be the standard normal measure defined on the real
line and B 2 (λ) be the associated Bayes Hilbert space. Let ν be the normal
2
Rmeasure 2with mean 0 and variance 2. Then f := dν/dλ ' exp(x /4), so
| log f | dλ < ∞ since the fourth moment of the standard normal measure
is finite. However, note that f itself represents an infinite measure; i.e. it is
not integrable. So in this case, an element of B 2 (λ) need not be a proper
density.
Aitchison (1982) introduced the space B 2 (λ) equipped with (⊕, ) and
h·, ·iB 2 (λ) in the special case where λ is the discrete uniform measure supported on a finite set K. The term Aitchison geometry is used to connote
the geometric properties of this Bayes Hilbert space. The term generalized
Aitchison geometry may be preferred when referring to Bayes Hilbert spaces
constructed with general reference measures.

2.2

Bounded linear operators on Bayes Hilbert space

Here we briefly summarize essential concepts on bounded linear operators.
The reader is referred to Bosq (2000) and Conway (1994) for a detailed
introduction.
Let H be B 2 (λ) or its complexification BC2 (λ), defined as
BC2 (λ) = {f ⊕ ig : f, g ∈ B 2 (λ)},
where i denotes the imaginary unit. BC2 (λ) is understood to be equipped
with the inner product hf ⊕ ig, f 0 ⊕ ig 0 iBC2 (λ) = hf, f 0 iB 2 (λ) − ihf, g 0 iB 2 (λ) +
ihg, f 0 iB 2 (λ) + hg, g 0 iB 2 (λ) .
Let B(H) be the space of bounded linear operators, equipped with the
usual operator norm
kAkB(H) = sup kAf kH ,

A ∈ B(H).

kf kH ≤ 1

A ∈ B(H) is said to be compact if for two orthonormal bases (uj , j ∈ N)
and (vj , j ∈ N) of H, it can be written as
Af =

∞
M

γj hf, uj iH

j=1

for some sequence (γj , j ∈ N) tending to zero.
7

vj ,

Given A ∈ B(H), we define two fundamental subspaces, the range and
kernel of A, as follows.
ker A := {f ∈ H : Af = λ},
ran A := {Af : f ∈ H}.
The dimension of ker A is called the nullity of A, and the dimension of ran A
is called the rank of A.
For A ∈ B(H), there exists an operator A∗ ∈ B(H), called the adjoint of
A, that is uniquely determined by the following property
hAf, giH = hf, A∗ giH ,

for all f, g ∈ H.

Given a subset V ⊂ H, the orthogonal complement of V is denoted by
V ⊥ = {g ∈ H : hf, giH = 0 for all f ∈ V }.
The closure of V , defined as the union of V and its limit points, is denoted
by cl(V ). It turns out that the following relationship holds (see e.g. Conway,
1994, pp. 35-36)
ker A = (ran A∗ )⊥

and (ker A)⊥ = cl(ran A∗ ).

(2.1)

A ∈ B(H) is said to be positive semidefinite if hAf, f iH ≥ 0 for all f ∈ H.
If the inequality is strict for all f 6= λ, A is said to be positive definite.
Given a subspace V ∈ H, A|V denotes the restriction of an operator
A ∈ B(H), i.e. A|V : V → H.
Let K = L2 (λ) if H = B 2 (λ), and let K = L2C (λ) if H = BC2 (λ), where
L2C (λ) is the complexification of L2 (λ) similarly defined as BC2 (λ). Let clrC
denote the complexification of clr, i.e. clrC (f ⊕ ig) = clr(f ) + i · clr(g). Then
clrC is an isometric isomorphism between BC2 (λ) and L2C (λ). Since clr or clrC
is an isometric isomorphism between H and K, any bounded linear operator
in B(H) can be understood as the corresponding element in B(K). This
property will be used repeatedly in subsequent sections.

2.3

Random elements of Bayes Hilbert space

Here we summarize some essential concepts relating to random elements of
a Bayes Hilbert space. See Bosq (2000) for a detailed discussion of random
elements of Banach and Hilbert spaces.
8

Let (Ω, F, P ) be the underlying probability triple and let H = B 2 (λ) or
which is understood to be equipped with its Borel σ-field. An Hvalued random variable is a measurable map X : Ω → H. Let k · kH be the
norm defined on H. We say that X is integrable if EkXkH < ∞. If X is
integrable, there exists a unique element, denoted by EX, such that
BC2 (λ),

EhX, f iH = hEX, f iH ,

for all f ∈ H.

EX is called the expectation of X.
If EkXk2H < ∞ then X is said to be square-integrable. Let L2H denote
the space of square-integrable random variables X with EX = λ, equipped
with the norm
1/2
kXkL2H = EkXk2H
.
It turns out that L2H is a Banach space. For X ∈ L2H , the Cauchy-Schwarz
inequality implies that Xhf, XiH is integrable for all f ∈ H. Define the
operator CX ∈ B(H) by
CX (f ) := E (Xhf, XiH ) .
CX is called the covariance operator of X.

2.4

Operator pencils

Let H be a complex separable Hilbert space and U be an open connected set
in the complex plane C. An operator-valued map Q : U → B(H) is called
an operator pencil. In this paper, H = BC2 (λ) or L2C (λ) is considered.
For notational convenience, we put H = L2C (λ). An operator pencil Q(z)
is holomorphic on an open connected set D if and only if
Q(1) (z0 ) := lim

z→z0

Q(z) − Q(z0 )
z − z0

exists in the norm of B(H) for each z0 ∈ D. A well known fact is that when
Q(z) is holomorphic at z = z0 , it allows the power series expansion
Q(z) =

∞
X

Qk (z − z0 )k

k=0

for some sequence (Qk , k ∈ N) in B(H).
The set σ(Q) := {z ∈ U : Q(z) is not invertible} is called the spectrum
of Q, which turns out to be closed (Markus, 2012, p. 56).
9

3

Cointegrated density-valued linear processes

Hereafter the following is always assumed without explicitly stating it.
Assumption (R). λ is the uniform measure on a compact interval K ⊂ R.
Under the above assumption, B 2 (λ) is the Bayes Hilbert space studied in
Egozcue et al. (2006). The assumption is not essential but convenient. In
fact, extending the subsequent results of this paper to other measures supported on K, such as the truncated normal measure, can be easily done with
trivial modifications.
First we define linear processes in B 2 (λ), called Bayes linear processes.
Due to the Hilbert space structure, linear processes in B 2 (λ) may be easily
defined and studied according to the work of Bosq (2000, 2007). As seen in
Section 2, a Bayes linear process may be understood as a linear process of
probability densities with compact support K.
Then, we introduce I(1) processes in B 2 (λ) based on the previous work of
Beare et al. (2017) who provided a notion of cointegrated linear processes in
an arbitrary Hilbert space. Of course, it can be understood as an I(1) process
of probability densities with support K. We further investigate cointegration
in B 2 (λ).

3.1

Linear processes in B 2 (λ)

Suppose that we have a time series of densities f = (ft , t ≥ t0 ) with compact
t
for each t, where νt is a probability
support K ⊂ R. We assume that ft := dν
dµ
measure for each t and µ is the Lebesgue measure defined on R. Since
ft has compact support K, it is absolutely continuous with respect to λ.
This implies that the corresponding λ-density of the underlying probability
measure νt is µ(K)−1 ft . Thus, the unit integral constraint singles out ft
itself from the equivalence class containing ft . The time series of µ-densities
f = (ft , t ≥ t0 ) may therefore be regarded as a stochastic process taking
values in B 2 (λ).
Let ε = (εt , t ∈ Z) be an independent and identically distributed (iid)
sequence in L2B 2 (λ) and (Ak , k ≥ 0) be a sequence in B(B 2 (λ)) satisfying
P∞
2
k=0 kAk kB(B 2 (λ)) < ∞. For fixed t0 ∈ Z ∪ {−∞}, the sequence f = (ft , t ≥

10

t0 ) defined as
ft =

∞
M

Ak εt−k

(3.1)

k=0

is convergent in L2B 2 (λ) (Bosq, 2000, Lemma 7.1). We call the sequence f =
(ft , t ≥ t0 ) a Bayes linear process. Bayes linear processes are necessarily
stationary.
P∞
A Bayes linear process,
L∞ (3.1), is said to be standard2 if k=0 kAk kB(B 2 (λ)) <
∞. In this case, A := k=0 Ak is convergent in B(B (λ)). As in Beare et al.
(2017), we define the long-run covariance operator Λ ∈ B(B 2 (λ)) for a standard linear process as follows
Λ := ACε0 A∗ .
Since clr is an isometric isomorphism, it follows that the clr image of a
(standard) linear process in B 2 (λ) is a (standard) linear process in L2 (λ).
Denote by g̃ ∈ L2 (λ) the clr image of g ∈ B 2 (λ). Then the L2 (λ)-linear
process paired with (3.1) is
f˜t =

∞
X

Ãk ε̃t−k ,

k=0

where f˜t = clr(ft ), Ãk = clr ◦Ak ◦ clr−1 and ε̃ = (ε˜t , t ∈ Z) is an iid sequence
in L2 (λ).

3.2

I(1) processes in B 2 (λ) and cointegration

Chang et al. (2016) were the first to study I(1) processes taking values in
an infinite dimensional Hilbert space, specifically the space L2 (λ). Subsequently, Beare et al. (2017) considered I(1) processes in an arbitrary separable complex Hilbert space. In this section we specialize the latter setting
to the Bayes Hilbert space B 2 (λ). Let the time series of random densities
of interest f = (ft , t ≥ 0) be a sequence in B 2 (λ). Denote the time series
of first differences ∆f = (∆ft , t ≥ 1) with ∆ft = ft ft−1 . We say that
f = (ft , t ≥ 0) is I(1) if ∆ft satisfies
∆ft =

∞
M
k=0

11

Nk εt−k

(3.2)

for all t ≥ 1, where ε = (εt , t ∈ Z) is an iid sequence in L2B 2 (λ) , and (Nk , k ≥ 0)
P∞
2
is
k=0 kkNk kB(B 2 (λ)) < ∞ and N :=
La∞ sequence in B(B (λ)) satisfying
k=0 Nk 6= 0.
From Chang et al. (2016) and Beare et al. (2017), it can be shown that
the above I(1) sequence allows the so-called Beveridge-Nelson decomposition,
ft = (f0

ν0 ) ⊕ N ξ t ⊕ νt ,

t≥0

Lt
L∞
L
where ξt =
( ∞
s=1 εs and νt =
k=0 Nk εt−k with Nk = −1
j=k+1 Nj ).
This means that ft is obtained by combining three different components: an
initial condition f0 ν0 , a random walk component N ξt , and a stationary
component νt .
Given an I(1) sequence f = (ft , t ≥ 0), we define the cointegrating space
associated with f , denoted by C(f, B 2 (λ)), to be the set


C(f, B 2 (λ)) := g ∈ B 2 (λ) : hft , giB 2 (λ) , t ≥ 0 is stationary for some f0 .
Moreover, we call A(f, B 2 (λ)) := C(f, B 2 (λ))⊥ the attractor space. These
terms to indicate such subspaces are commonly used in the literature on
cointegration, see e.g. Johansen (1996).
Define the long-run variance operator of ∆f as Λ := N Cε0 N ∗ . If Cε0 is
positive definite, it must be the case that
ker Λ = ker N ∗ ,
and further Beare et al. (2017) showed that the cointegrating space is equal
to ker N ∗ .
The cointegrating space (or attractor space) can be identified by analyzing
its clr-image in L2 (λ), which may be convenient in practice. Suppose that
the clr image of (3.2) is
∆f˜t =

∞
X

Ñk ε̃t−k ,

k=0

where ∆f˜t = f˜t − f˜t−1 and Ñk = clr ◦Nk ◦ clr−1 . The Beveridge-Nelson
decomposition is
f˜t = (f˜0 − ν̃0 ) + Ñ ξ˜t + ν˜t ,
12

t≥0

where Ñ = clr ◦N ◦ clr−1 . We also define the long-run variance operator
of ∆f˜ as Λ̃ = Ñ Cε̃0 Ñ ∗ , where Cε̃0 is the covariance operator of ε̃0 . The
following result shows that we can fully identify the cointegrating space in
B 2 (λ) associated with f from the cointegrating space in L2 (λ) associated
with f˜, the clr image of f .
Proposition 3.1. If Cε0 is positive definite, then C(f, B 2 (λ)) = clr−1 ker Ñ ∗
Proof. Proposition 3.1 in Beare et al. (2017) implies that C(f, B 2 (λ)) =
ker N ∗ . Note that clr∗ = clr−1 since for g̃ := clr(g),
hf, clr−1 g̃iB 2 (λ) = hclr(f ), g̃iL2 (λ) .
This holds for any arbitrary g ∈ B 2 (λ). By the uniqueness of the adjoint
map, clr∗ = clr−1 . From these results, it is clear that
Ñ ∗ = clr ◦N ∗ ◦ clr−1 .
Since clr and clr−1 are bijective, it follows that ker N ∗ = clr−1 ker Ñ ∗ .
Chang et al. (2016) considered a cointegrated density-valued process with
values in L2 (λ). However, Beare (2017) pointed out that the attractor space
under their assumptions is always trivial, i.e. A(f, L2 (λ)) = {0}. This problem occurs because an I(1) process of probability densities in L2 (λ) cannot
satisfy the nonnegativity constraints required for each realization to be a
proper density; the space of probability densities with compact support K
is strictly smaller than L2 (λ) and not a subspace. Due to the geometry of
the Bayes Hilbert space, a cointegrated linear process in B 2 (λ) does not suffer from this problem. Furthermore, it is clear from the earlier discussion
of Bayes Hilbert spaces in Section 2.1 that the space of probability densities supported on K is isomorphic to L2 (λ), which is just a slightly smaller
subspace than L2 (λ) itself.

4

Cointegrated AR(p) processes in B 2(λ)

In this section, we study AR(p) processes with values in B 2 (λ) and obtain a
version of the Granger-Johansen representation theorem.
Suppose that a sequence of λ-densities (ft , t ≥ −p + 1) in B 2 (λ) evolves
according to
ft = A1 ft−1 ⊕ · · · ⊕ Ap ft−p ⊕ εt
13

(4.1)

for t ≥ 1, where ε = (εt , t ∈ Z) is an iid sequence in L2B 2 (λ) and A1 , . . . , Ap ∈
B(B 2 (λ)). Equation (4.1) may be understood as an autoregressive law of
motion for density-valued processes in B 2 (λ). The corresponding clr image
of (4.1) with values in L2 (λ) is given by
f˜t = Ã1 f˜t−1 + · · · + Ãp f˜t−p + ε˜t

(4.2)

for t ≥ 1, where Ãk = clr ◦Ak ◦ clr−1 . It is clear that Ãk ∈ B(L2 (λ)), and
that Ãk is compact if and only if Ak is compact. Let idL2 (λ) be the identity
operator on L2 (λ), and for z ∈ C define the operator pencil
Ã(z) = idL2 (λ) −z Ã1 − · · · − z p Ãp .
For each z ∈ C, Ã(z) is a bounded linear operator on L2C (λ), the complexification of L2 (λ). The operator pencil Ã(z) is polynomial in z; see Markus
(2012) for a detailed discussion of polynomial operator pencils. Let idB 2 (λ)
be the identity operator on B 2 (λ), and define
A(z) = idB 2 (λ) (z

A1 )

···

(z p

Ap ).

(4.3)

For each z, A(z) is a bounded linear operator on BC2 (λ). One can easily show
that A(z) = clr−1
C ◦Ã(z) ◦ clrC , so σ(A) = σ(Ã) holds.
To further analyze (4.1), we employ the following mutually exclusive assumptions.
Assumption (S). If z ∈ σ(A) then |z| > 1.
Let Dr be the open disk centered at zero with radius r > 0. Since σ(A)
is closed, Assumption (S) implies that there exists η > 0 such that A(z) is
invertible on D1+η .
Assumption (N).
(i) A1 , . . . , Ap ∈ B(B 2 (λ)) are compact operators.
(ii) If z ∈ σ(A) then |z| > 1 or z = 1. Moreover, 1 ∈ σ(A).
(iii) If P ∈ B(B 2 (λ)) is a projection on ran A(1), then the map
(idB 2 (λ) P)A(1) (1)|ker A(1) : ker A(1) → ran(idB 2 (λ) P)
is invertible.
14

Under Assumption (N)-(i), A(z) is an index-zero Fredholm operator for each
z ∈ C. This implies that ran A(1) is closed. (N)-(i) is not restrictive in practice since it is in fact common to assume the compactness of autoregressive
operators. (N)-(ii) implies the existence of η > 0 such that A(z) is invertible
everywhere on D1+η except at z = 1. The role of (N)-(iii) will be discussed
later.
The AR(p) law of motion (4.1) behaves very differently depending on
whether A(z) satisfies Assumption (S) or (N). In the former case it generates
a stationary process which may be represented as a standard Bayes linear
process.
Proposition 4.1. Suppose that Assumption (S) is satisfied. Then the AR(p)
law of motion (4.1) admits a unique stationary solution in L2B 2 (λ) given by
ft =

∞
M

Nk εt−k

(4.4)

k=0

for all t ≥ 1. Moreover, Nk is given by
Nk =

1
k!

N (k) (0)

and N (z) := A(z)−1 is holomorphic on D1+η for some η > 0.
Proof of Proposition 4.1. Consider the clr image of (4.2) in L2 (λ) and let
H = L2 (λ). Using the Markovian representation in Bosq (2000, p. 128), the
AR(p) law of motion in H can be re-expressed as an AR(1) law of motion in
Hp by writing
Ft = AFt−1 + Et ,
where

f˜t
 f˜ 
 t−1 
Ft =  .  ,
 .. 
f˜t−p+1


 
ε˜t
0
 
Et =  ..  ,
.
0




Ã1 Ã2 · · · Ãp−1 Ãp
id
0
0
 H 0 ···

A= .
.. . .
..
..  .
.
 .
.
.
.
. 
0
0 · · · idH 0

15

(4.5)

Define A(z) := idHp −Az, i.e.

idH −Ã1 z −Ã2 z −Ã3 z −Ã4 z
···
−Ãp z
 − id z
idH
0
0
···
0
H


0
− idH z idH
0
···
0
A(z) = 

.
.
.
.
..
.
..
..
..
..
..

.
0
0
0
0
− idH z idH


A11 (z) A12 (z)
=:
,
A21 (z) A22 (z)









where one can easily verify that A22 (z) is invertible. Define the Schur com−1
plement of A22 (z) as A+
11 (z) := A11 (z) − A12 (z)A22 A21 (z). From a little
algebra, it can be easily verified that A+
11 (z) = Ã(z). Since A22 (z) is invert+
ible, A(z) is invertible if and only if A11 is invertible. Therefore, it follows
that σ(A) = σ(Ã) ⊂ {z : |z| > 1}. Define
s(A) := {z ∈ C : A − z idHp is invertible}
and
ρ(A) := sup{|z| : z ∈ σ(A)},
the spectrum and spectral radius of the bounded linear operator A. From
the construction of the operator pencil A(z) and Assumption (S), it is clear
that
ρ(A) := sup{|z| : z ∈ s(A)} < 1.

(4.6)

In view of Gelfand’s formula (Conway, 1994, Proposition VII.3.8), (4.6) implies that there exists j such that
kAj kB(Hp ) < aj ,

a ∈ (0, 1).

(4.7)

From a well known result on the inverse of Banach-valued functions , A(z)−1
is given as
A(z)−1 = idH +Az + A2 z 2 + · · ·
which is convergent in D1+η , and indeed holomorphic on D1+η , for some η > 0
because of (4.7). This shows that the AR(1) law of motion in (4.5) may be
represented as
Ft =

∞
X

Ak Et−k ,

k=0

16

t ∈ Z.

Let Π : Hp → H be the coordinate projection given as Π(x1 , . . . , xp ) = x1
and let Π∗ : H → Hp be its adjoint. Then one can easily verify that
f˜t =

∞
X

ΠAk Π∗ ε̃t−k ,

t ∈ Z.

k=0
−1
Simple algebra yields ΠA(z)−1 Π∗ = A+
= Ã(z)−1 , and it is clear that
11 (z)
ΠAk Π∗ is the kth coefficient in the Taylor series of Ã(z)−1 around z = 0.
Holomorphicity of A(z)−1 on D1+η is inherited from the holomorphicity of
Ã(z)−1 on D1+η , which is implied by the holomorphicity of A(z)−1 established
above and the relation Ã(z)−1 = ΠA(z)−1 Π∗ .
P
Remark 4.1. Norm-summability, ∞
k=0 kÑk kB(B 2 (λ)) < ∞, is a natural consequence of holomorphicity of A(z)−1 on D1+η for some η > 0. This shows
that (4.4) is a standard Bayes linear process.

Under Assumption (N), a process in B 2 (λ) satisfying the AR(p) law of
motion (4.1) is no longer stationary. Loosely speaking, it corresponds to an
I(d) process for d ≥ 1.
Proposition 4.2. Under Assumption (N)-(i, ii), the AR(p) law of motion
does not admit a stationary solution. Instead, given a suitable initialization
f−d+1 , . . . , f0 , it is satisfied by a sequence in L2B 2 (λ) whose dth difference is
given by
∞
M
d
Nk εt−k
∆ ft =
k=0

for all t ≥ 1. Nk is defined by
Nk =

1
k!

N (k) (0)

and N (z) := (z − 1)d A(z)−1 is holomorphic on D1+η for some η > 0.
If Assumption (N)-(iii) also holds, then d = 1 and
lim (z − 1)

z→1

A(z)−1 =




−1

idB 2 (λ) P A(1) (1)|ker A(1)
idB 2 (λ) P , (4.8)

which does not depend on the choice of the projection P on ran A(1).

17

Proof. For notational convenience, we work with Ã(z) and set H = L2 (λ).
Assumption (N)-(i,ii) implies that 1 is an isolated singularity of Ã(z). From
Lemma 8.1, we know that there exist finite dimensional projections Q1 , . . . , Qd
such that, for z in a punctured neighborhood of 1,
Ã(z) = (P1 + (1 − z)Q1 ) · · · (Pd + (1 − z)Qd )G(z),
where Pi = idH −Qi for i = 1, . . . , d, and G(z) is a holomorphic function that
is invertible at z = 1. Note that for i = 1, . . . , d
(Pi + (1 − z)Qi )−1 = (Pi + (1 − z)−1 Qi ).
This shows that in a punctured neighborhood of z = 1,
Ã(z)−1 = G(z)−1 (Pd + (1 − z)−1 Qd ) · · · (P1 + (1 − z)−1 Q1 ).

(4.9)

Therefore, Ã(z)−1 is meromorphic and the lowest order of the Laurent series
is −d. Therefore Ñ (z) := (1 − z)d Ã(z)−1 can be continued holomorphically
over z = 1. By applying the linear filter induced by Ñ (z) to (4.1), we have
for all t ≥ d + 1
∆ f˜t =
d

∞
X

Ñk ε̃t−k ,

k=0

Ñk =

1 (k)
Ñ (0)
k!

(4.10)

Since (∆d f˜t , t ≥ d + 1) is stationary, it may be extended for t ≥ 1 in a
natural way and f−d+1 , . . . , f0 may be chosen accordingly. This proves the
first statement.
Additionally under Assumption (N)-(iii), Lemma 8.2 implies that d = 1.
The only remaining thing is to verify the residue formula (4.8). From (4.9)
we know that
lim (1 − z)Ã(z)−1 = G(1)−1 (idH −P1 ).

z→1

Define
A := G(1)−1 |ran(idH −P1 ) ,
B := G(1)|ker Ã(1) .
From Lemma 8.3(a,b) we know that B = −(idH −P1 )Ã(1) (1)|ker Ã(1) and
ran B ⊂ ran(idH −P1 ). Since B as a map from ker Ã(1) to ran(idH −P1 )
18

is invertible under Assumption (N)-(iii), it is clear that AB = idker Ã(1) . This
shows that
Ñ−1 := lim(1 − z)Ã(z)−1 = A(idH −P1 ) = B−1 (idH −P1 ).
z→1

Under Lemma 8.1, we may choose P1 to be any projection on ran Ã(1), so
the expression for Ñ−1 just obtained cannot depend on the particular choice
of projection.
Under Assumption (N)-(i,ii), Proposition 4.2 says that ∆d f = (∆d ft , t ≥ 1)
is a Bayes linear process under a suitable initial condition. In this sense, we
say that the AR(p) law of motion (4.1) is I(d) for some d ≥ 1. The most
common case in practice may be I(1), and such an I(1) solution is guaranteed
under an extra condition, Assumption (N)-(iii).
Remark 4.2. When H = Rn , Johansen (1991) provided a necessary and sufficient condition on the autoregressive matrix polynomial and its first derivative at one for an AR(p) law of motion to be I(1). It is commonly called the
I(1) condition. Assumption (N)-(iii) plays the same role for an AR(p) law
of motion in B 2 (λ). Beare et al. (2017) provided an I(1) condition for an
AR(p) law of motion in an arbitrary complex Hilbert space, and Assumption
(N)-(iii) is a reformulation of it.
Remark 4.3. When p = 1, an I(1) solution can be guaranteed without
requiring compactness of the autoregressive operator. See Theorem 4.1 in
Beare et al. (2017).
Under Assumption (N)-(i,ii,iii), one natural consequence of Proposition
4.2 is the following Beveridge-Nelson decomposition: for some f0 , ν0 ∈ B 2 (λ),
ft = (f0

ν0 ) ⊕ N (1)ξt ⊕ νt ,

t≥0

(4.11)

Lt
L∞
where ξt =
s=1 εs , νt =
k=0 Nk εt−k and N (1) is a finite rank operator
which can be explicitly obtained from the residue formula (4.8). If ε = (εt , t ∈
Z) has a positive definite covariance operator, one can easily identify the
cointegrating and attractor spaces of B 2 (λ) from the decomposition (4.11).
Proposition 4.3. Assume that Cε0 is positive definite and (N)-(i, ii, iii) hold.
Then C(f, B 2 (λ)) = [ker A(1)]⊥ .

19

Proof. We only provide an informal argument. For a more detailed and
rigorous proof, refer to Beare et al. (2017, Proposition 3.1). From (4.11), for
any g ∈ B 2 (λ) we have
hft , giB 2 (λ) = hf0

ν0 , giB 2 (λ) + hN (1)ξt , giB 2 (λ) + hνt , giB 2 (λ) .

By employing the initial condition f0 = ν0 , we can make the first term
vanish. Moreover, the inner product process (hνt , giB 2 (λ) , t ≥ 0) is stationary
because (νt , t ≥ 0) is stationary. The variance of hN (1)ξt , giB 2 (λ) increases in
t. Thus, hN (1)ξt , giB 2 (λ) must vanish in order for hft , giB 2 (λ) to be stationary.
Given that ε has a positive definite covariance operator, for hN (1)ξt , giB 2 (λ)
to vanish we require g ∈ ker N (1)∗ . From the closed form solution for N (1)
given in (4.8) combined with (2.1) , it can be easily deduced that ker N (1)∗ =
[ker A(1)]⊥ .
Naturally, the attractor space A(f, B 2 (λ)) = C(f, B 2 (λ))⊥ = ker A(1), which
is finite dimensional under Assumption (N).
Example 4.1 (A numerical example).
Let λ be the uniform measure on [−3, 3] and consider the Bayes Hilbert
space associated with λ. Let g be the standard normal measure truncated
on [−3, 3]. We consider the following AR(1) law of motion
g) = Ψ(ft−1 g) ⊕ εt ,
∞
M
Ψ(·) =
λj h·, ej iB 2 (λ) ej ,
(ft

j=1

where λj = 21−j . Let (uj , j ∈ N) be the Fourier basis for L2 (λ), and let e1 be
the Cauchy probability density with location zero and scale 0.25. We obtain
an orthonormal basis (ej , j ∈ N) for B 2 (λ) by applying the Gram-Schmidt
process to the basis (e1 , clr−1 u2 , clr−1 u3 , . . .). If εt has a positive definite
covariance operator, this specification makes the attractor space equal to
A(f, B 2 (λ)) = {g ⊕ (a

e1 ) : a ∈ R}.

(4.12)

Figure 4.1 shows typical elements in the attractor space. An element g ⊕
(a e1 ) shows higher concentration around 0 when a > 0, or is bimodal
when a < 0. To obtain a simulated sequence, I generated a sequence of
independent standard Brownian bridges (Bt , t ≥ 1), and set
εt = clr−1 {0.3 · Pm (λ) B t },
20

where Pm (λ) denotes orthogonal projection on theR space of mth-order poly1
nomials in L2 (λ), and B t (u) = Bt ((u + 3)/6) − 6 0 Bt (r)dr for u ∈ [−3, 3].
A consequence of using the projection Pm to construct the innovation εt is
that the attractor space may not be exactly equal to that in (4.12), since the
covariance operator of εt is not positive definite. However, this is a convenient way to force the simulated probability densities to be smooth, and by
choosing large m the attractor space should approximate that in (4.12). I
set f1 = g and generated a sequence which is shown in Figure 4.2. We can
easily see that ft floats around the attractor space in the long run.
Figure 4.1: Typical elements in the attractor
g ⊕ e1

g

e1

Figure 4.2: Realization of (ft , t ≥ 1)
t=1

t = 10

t = 30

t = 100

t = 300

t = 500

t = 800

t = 1000

21

5

Statistical inference

In this section we provide a statistical procedure to estimate the attractor
space, which is assumed here to be finite dimensional. Chang et al. (2016)
provided statistical methods based on FPCA for a cointegrated densityvalued linear process with values in L2 (λ). Even though Beare (2017) commented that their time series cannot accommodate a nontrivial attractor
space, the statistical procedure they provided is useful since we can apply it
to the time series of clr-images f˜ = (f˜t , t ≥ 1) with values in L2 (λ). Using
their procedure, the attractor space of L2 (λ) can be estimated, and then it
is easy to obtain the attractor space of B 2 (λ) using the clr−1 transformation.
We therefore do not have to develop a new statistical technique adapted to
Bayes Hilbert spaces.
Assumption (T1).
(i) For some finite rankPÑ ∈ B(L2 (λ)) and a sequence (Ñk , k ≥ 0) in
B(L2 (λ)) such that ∞
k=0 kkÑk kL2 (λ) < ∞, we have
f˜t = f˜0 − ν̃0 + Ñ ξ˜t + ν̃t ,
where ν̃t =

P∞

k=0

t ≥ 0.

Ñk ε̃t−k .

(ii) (ε̃t , t ∈ Z) is an iid sequence with E ε̃t = 0, positive definite covariance
operator Cε̃0 , and Ekε̃t kpL2 (λ) < ∞ for some p ≥ 4.
In the above assumption, it is explicitly required that the attractor space
A(f˜, L2 (λ)) is finite dimensional. An example of such an I(1) process is
an AR(p) process satisfying Assumption (N), which we considered in the
previous section. Assume that dim(A(f˜, L2 (λ))) = r > 0 and let (vi , i ∈ N)
be the orthonormal basis of L2 (λ) such that
A(f˜, L2 (λ)) = span(v1 , . . . , vr ),
C(f˜, L2 (λ)) = span(vr+1 , vr+2 , . . .).
The probability density functions (ft , t = 1, . . . , T ) are not directly observed and must be estimated from the data. We assume that there are n
cross-sectional observations that are available to estimate ft at each time
period t. One possibility would be to obtain estimated densities (fˆ1 , . . . , fˆT )
22

from a standard nonparametric kernel method. However, since we would
ˆ
Rbe applying the procedure of Chang et al. (2016) to the clr-images log ft −
ˆ
log ft dλ, this naive approach could lead to unacceptable bias if the densities are close to zero at the boundary of their support. Possibly, a better
approach is to estimate the log-densities directly. Methods for doing so include the maximum penalized likelihood method of Silverman (1982), the
spline-based method of O’Sullivan (1988) and the local likelihood method of
Loader (1996). However the clr-images are estimated, we assume that they
satisfy the following high level condition.
Assumption (T2). The estimated clr-images (f˜tE , t = 1, . . . , T ) satisfy
(i) sup1≤t≤T kf˜tE − f˜t kL2 (λ) = Op (1),
(ii) T −1

PT

t=1

kf˜tE − f˜t kL2 (λ) →p 0.

The testing procedure is based on FPCA of the empirical covariance operator
!
!
T
T
T
X
X
X
T
−1
E
−1
E
E
−1
E
Ṽ = T
f˜t − T
f˜t ⊗ f˜t − T
f˜t .
t=1

t=1

t=1

Let (ζiT , viT ), i ∈ N, be the eigenpairs of the empirical covariance operator
Ṽ T , where ζiT is assumed to be decreasing as i gets larger. Let PA be the
orthogonal projection on ran Ñ and define PC = idL2 (λ) −PA . Ṽ T can be
written as
T
T
T
T
Ṽ T = ṼAA
+ ṼAC
+ ṼCA
+ ṼCC
,
T
T
T
T
where ṼAA
= PA Ṽ T PA , ṼAC
= PA Ṽ T PC , ṼCA
= PC Ṽ T PA , and ṼCC
=
PC Ṽ T PC . Under Assumptions (T1) and (T2), it can be shown that
T
kB(L2 (λ)) = Op (T −1 ).
kT −1 Ṽ T − T −1 ṼAA
T
That is, the difference between T −1 Ṽ T and T −1 ṼAA
is a norm-vanishing opT
erator. Noting that the eigenfunctions of ṼAA are in A(f˜, L2 (λ)), it follows
that the eigenfunctions associated with the r leading eigenvalues are a consistent estimator of a spanning set of A(f˜, L2 (λ)). Formally, it can be shown
(Chang et al., 2016, Proposition 3.2 and Section 4) that

kP̂A − PA kL2 (λ) = Op (T −1 ),
23

(5.1)

where P̂A is the projection operator based on the eigenfunctions. An alternative formulation of (5.1) is
span(v1T , . . . , vrT ) →p A(f˜, L2 (λ)).

(5.2)

Needless to say, C(f˜, L2 (λ)) can be estimated as the orthogonal complement,
[span(v1T , . . . , vrT )]⊥ . These results show that if r is known, then it is easy to
estimate the attractor space, A(f˜, L2 (λ)), since the r leading eigenfunctions
of Ṽ T asymptotically span A(f˜, L2 (λ)).
In practice, the dimension of A(f˜, L2 (λ)) is typically unknown. To determine the dimension, we may consider testing the null hypothesis


2
˜
H0 : dim A(f , L (λ)) = R,
(5.3)
against the alternative


H1 : dim A(f˜, L2 (λ)) ≤ R − 1.

(5.4)

For some positive integer Rmax ≥ 1, successive tests of the null hypothesis
for R = Rmax , Rmax − 1, . . . , 1 can determine the dimension of A(f˜, L2 (λ)).
For example, if the null is rejected in favor of the alternative
for R > r, but

2
˜
not rejected for R = r , we can conclude that dim(A f , L (λ)) = r. A
feasible test statistic and its limiting distribution under the null hypothesis
are given by Chang et al. (2016). Let
*
+
*
+ !0
T
T
X
X
ẑt =
v T , f˜E − T −1
f˜E
, . . . , v T , f˜E − T −1
f˜E
1

t

t

t=1

R

H

t

t

t=1

H

for t =P
1, . . . , T , and define ẐT = (ẑ1 , . . . , ẑT )0 . Further let Q̂TR = ẐT0 ẐT and
Σ̂TR = |k|≤` w` (k)Γ̂(k), where Γ̂(k) is the sample autocovariance of ẑt and
w` (k) is a bounded weight function1 . The proposed test statistic is
τRT = T −2 ζmin (Q̂TR , Σ̂TR )
where ζmin (Q̂TR , Σ̂TR ) denotes the smallest generalized eigenvalue of Q̂TR with
respect to Σ̂TR .
1

See Chang et al. (2016) for further details. The truncation parameter ` may be chosen
as in Andrews (1991). In the empirical applications reported in Section 6 of this paper,
w` (·) is chosen to be the quadratic spectral kernel.

24

Proposition 5.1 (Theorem 4.3 in Chang et al., 2016). If Assumptions (T1)
and (T2) are satisfied then, under the null hypothesis H0 in (5.3), we have
Z

0
T
W R (r)W R (r) dr, IR ,
(5.5)
τR →d ζmin
R
where I is the R × R identity matrix, W R (r) = WR (r) − WR (s)ds, and WR
is an R × 1 vector of independent standard Brownian motions. On the other
hand, under the alternative hypothesis H1 in (5.4),
τRT →p 0.
R
That is, τRT converges in law to the smallest eigenvalue of W R (r)W R (r)0 dr
under the null hypothesis, and vanishes under the alternative. The critical
values can be easily obtained from a large number of simulated sample paths
of W R . The critical values are presented in Table 5.1.2
Table 5.1: Critical values of τRT

R
1%
5%
10%

1
2
3
4
5
6
7
0.0248 0.0163 0.0123 0.0100 0.0084 0.0073 0.0065
0.0365 0.0215 0.0156 0.0122 0.0101 0.0086 0.0075
0.0459 0.0254 0.0177 0.0136 0.0111 0.0094 0.0081

Remark 5.1. The case when f˜t is observable at each t may be easily dealt
with. Test statistics can be calculated based on true observations f˜ = (f˜t , t =
1, . . . , T ). The limit distribution is slightly different from (5.5); it is given by
Z

0
ζmin
WR (r)WR (r) dr, IR .
See Chang et al. (2016) for more details.

6

Empirical application

6.1

Example 1: Cross-sectional densities of earnings

In this section, we revisit the application of FPCA to cross-sectional densities of individual weekly earnings undertaken by Chang et al. (2016). The
2

The reported critical values for R = 1, . . . , 5 are taken from Table 1 in Chang et al.
(2016). The values for R = 6 and 7 are newly calculated for reference in the next section.

25

cross-sectional observations that are used to estimate densities are obtained
at monthly frequency from the Current Population Survey (CPS) database,
running from January 1990 to March 2017 (327 months in total). Since individual earnings in each month are reported in current dollars, I adjusted them
all to January 1990 prices.3 I excluded earnings below the 2.5th percentile
and above the 97.5th percentile since there are many near-zero4 earnings and
top-coded earnings in the raw dataset.5 Weekly earnings are censored from
above, with the threshold for censoring changing partway through the sample: the top-coded nominal earning is $1923 before January 1998, and $2885
afterward. This difference introduces significant heterogeneity; for example,
the support of the earnings distribution changes greatly after January 1998
since the number of observations contained in (1923, 2885) is always zero
before this month, but nonzero afterward. However, after dropping the top
2.5% of observations, the range of earnings is quite stable over time. Moreover, after dropping the bottom 2.5% of observations, all near-zero earnings
are excluded and the smallest earnings become reasonably sized. After applying this truncation, the number of observations for each month ranges
from 11760 to 15489. The last thing we need to notice is that the CPS individual earnings data are collected from a monthly survey of individuals, with
each individual assigned their own design weight. Therefore, any estimates
constructed from the dataset should take design weights into account. To
estimate log-densities, I used the local likelihood method of Loader (2006)
which can easily accommodate design weights. The procedure is described in
detail in the Appendix. Figure 6.2 shows the time series of density estimates
and their clr-images.
To investigate the dimension of the attractor space, I calculated the test
statistics described in the previous section. Figure 6.1 displays the test statistics for R = 1, . . . , 7 as well as a scree plot of the eigenvalues. From the scree
plot it is apparent that Rmax can be set to around 5. Referring back to Table
5.1, we see that R = 2 is rejected even at the 1% level, but R = 1 is not
rejected at the 5% level. We tentatively conclude that the dimension of the
2
attractor space is one. Viewed as a subspace of L (λ), the estimated attrac3

Monthly CPI data are obtained from the Federal Reserve Economic Data (FRED).
The dataset contains a lot of abnormal values of nominal earnings near zero. For
example, there are total 9188 observations corresponding to weekly earnings less than
$0.01 over the whole span.
5
This symmetrical trimming was considered in, for example, Autor et al. (2008) for
analysis of the CPS wage data.
4

26

tor space is the span of the leading eigenfunction v̂1 , displayed in Figure 6.1.
Alternatively, we may regard the attractor space to be the span of clr−1 (v̂1 ),
a subspace of B 2 (λ).
To describe the attractor space, I orthogonally projected the estimated
densities upon the estimated cointegrating space, obtaining the projected
series (fˆtC , t = 1, . . . , T ). I positively and negatively perturbed the sample
mean f¯TC of (fˆtC , t = 1, . . . , T ), called the stationary mean of f , in the direction associated with the leading eigenfunction clr−1 v̂1 . More specifically, the
positive perturbation is f¯TC ⊕ (ζ̂1 clr−1 (v̂1 )) and the negative perturbation
is f¯TC (ζ̂1 clr−1 (v̂1 )), where ζ̂1 is the eigenvalue associated with v̂1 . Figure
6.3 shows the results.
Figure 6.1: Test statistics — individual earnings

R
τRT

1
0.03638

2
0.01253

3
0.00908

4
0.00590

Scree plot of eigenvalues

6.2

5
0.00537

6
0.00506

7
0.00472

Leading eigenfunction

Example 2 : cross-sectional densities of wages

Now we consider cross-sectional densities of hourly wages. The cross-sectional
observations are obtained at monthly frequency from the CPS database, and
the span is the same as in the previous example. Wages are deflated using
CPI data and measured in January 1990 dollars. There are many near-zero
values and top-coded values in the raw data set, so I excluded earnings below
the 1.25th percentile and above the 98.75th percentile; unlike the earnings

27

Figure 6.2: Density estimates and clr-images - earnings

Time series of density estimates

Time series of clr-images

Figure 6.3: Description of the attractor - earnings

data set, the top-coded values are constant over time, so only a small truncation is required to get rid of abnormal or extreme values. The number
of observations for each month ranges from 7348 to 9584 after truncation.
Figure 6.5 shows the time series of density estimates and their clr-images.
Test statistics are calculated as before and reported in Figure 6.4 with
28

other graphs. From the scree plot, we may set Rmax to around 5. From
successive tests we conclude that the dimension of the attractor space is 1.6
Thus, the estimated attractor space, viewed as a subspace of L2 (λ), is the
span of the eigenfunction associated with the first leading eigenvalue. Viewed
as a subspace of B 2 (λ), it is the span of the clr−1 -image of this eigenfunction.
Figure 6.6 provides a visual impression of the attractor space, similar to
Figure 6.3 above. In this case I perturbed the stationary mean f¯TC in the
direction of clr−1 v̂1 by ±2ζ̂1 , so as to more clearly emphasize the difference
between the two perturbations.

7

Concluding remarks

In this paper we have investigated cointegrated linear processes with values
in a Bayes Hilbert space of densities. Autoregressive density-valued processes
were also studied and a version of the Granger-Johansen representation theorem is provided. We showed that the statistical methods developed by Chang
et al. (2016) can be used to estimate the attractor space associated with a
cointegrated linear process in a Bayes Hilbert space.
6

On the other hand, if we set Rmax equal to 6 (resp. 7) then our sequential testing
procedure indicates an attractor space dimension of 6 (resp. 7), which seems large.

Figure 6.4: Test statistics - wages

R
τRT

1
0.05174

2
0.01181

3
0.01081

4
0.01076

Scree plot of eigenvalues

5
0.00982

6
0.00977

7
0.00966

Leading eigenfunction

29

Figure 6.5: Density estimates and clr-images - wages

Time series of density estimates

Time series of clr-images

Figure 6.6: Description of the attractor - wages

30

References
Aitchison, J. (1982). The statistical analysis of compositional data. Journal
of the Royal Statistical Society. Series B (Methodological), 44(2):139–177.
Andrews, D. W. K. (1991). Heteroskedasticity and autocorrelation consistent
covariance matrix estimation. Econometrica, 59(3):817–858.
Autor, D. H., Katz, L. F., and Kearney, M. S. (2008). Trends in u.s. wage
inequality: Revising the revisionists. Review of Economics and Statistics,
90(2):300–323.
Beare, B. K. (2017). The Chang-Kim-Park model of cointegrated densityvalued time series cannot accommodate a stochastic trend. Econ Journal
Watch, 14(2):133 – 137.
Beare, B. K., Seo, J., and Seo, W.-K. (2017). Cointegrated linear processes
in hilbert space. Journal of Time Series Analysis, 38(6):1010–1027.
Beare, B. K. and Seo, W.-K. (2017). Representation of I(1) autoregressive
Hilbertian processes. ArXiv e-print, arXiv:1701.08149v1 [math.ST].
Bosq, D. (2000). Linear Processes in Function Spaces. Springer-Verlag New
York.
Bosq, D. (2007). General linear processes in Hilbert spaces and prediction.
Journal of Statistical Planning and Inference, 137(3):879 – 894. Special
Issue on Nonparametric Statistics and Related Topics: In honor of M.L.
Puri.
Chang, Y., Kim, C. S., and Park, J. Y. (2016). Nonstationarity in time series
of state densities. Journal of Econometrics, 192(1):152 – 167.
Conway, J. B. (1994). A Course in Functional Analysis. Springer.
Egozcue, J. J., Dı́az-Barrero, J. L., and Pawlowsky-Glahn, V. (2006). Hilbert
space of probability density functions based on Aitchison geometry. Acta
Mathematica Sinica, 22(4):1175–1182.
Granger, C. W. J. (1981). Some properties of time series data and their use
in econometric model specification. Journal of Econometrics, 16(1):121 –
130.
31

Horváth, L. and Kokoszka, P. (2012). Inference for Functional Data with
Applications. Springer-Verlag GmbH.
Howland, J. S. (1971). Simple poles of operator-valued functions. Journal of
Mathematical Analysis and Applications, 36(1):12 – 21.
Johansen, S. (1991). Estimation and hypothesis testing of cointegration vectors in Gaussian vector autoregressive models. Econometrica, 59(6):1551–
1580.
Johansen, S. (1996). Likelihood-Based Inference in Cointegrated Vector Autoregressive Models (Advanced Texts in Econometrics). Oxford University
Press.
Loader, C. R. (1996). Local likelihood density estimation. Annals of Statistics, 24(4):1602–1618.
Loader, C. R. (2006). Local Regression and Likelihood. Springer New York.
Markus, A. S. (2012). Introduction to the Spectral Theory of Polynomial
Operator Pencils (Translations of Mathematical Monographs). American
Mathematical Society.
O’Sullivan, F. (1988). Fast computation of fully automated log-density and
log-hazard estimators. SIAM Journal on Scientific and Statistical Computing, 9(2):363–379.
Petersen, A. and Müller, H.-G. (2016). Functional data analysis for density functions by transformation to a Hilbert space. Annals of Statistics,
44(1):183–218.
Silverman, B. W. (1982). On the estimation of a probability density function by the maximum penalized likelihood method. Annals of Statistics,
10(3):795–810.
van den Boogaart, K. G., Egozcue, J. J., and Pawlowsky-Glahn, V. (2014).
Bayes Hilbert spaces. Australian & New Zealand Journal of Statistics,
56(2):171–194.

32

8
8.1

Appendix
Useful lemmas

Lemma 8.1 (Theorem 1.4 in Howland (1971)). Let B(H) be the space of
bounded linear operators H → H for a separable complex Hilbert space H and
let A(z) = idH −C(z) where C(z) is an analytic family of compact operators.
If z0 ∈ σ(A) is an isolated element, there exist finite dimensional projections
Q1 , . . . , Qd such that
A(z) = (P1 + (z0 − z)Q1 ) · · · (Pd + (z0 − z)Qd )G(z),
where G(z) is analytic, G(z0 ) is invertible, and Pi = idH −Qi . We may
choose P1 to be any projection on ran A(z0 ).
Lemma 8.2. Under the conditions of Lemma 8.1,
(idH −P1 )A(1) (z0 )|ker A(z0 ) : ker A(z0 ) → ran(idH −P1 )
is invertible if and only if A(z)−1 has a simple pole at z = z0 .
Proof. From Corollary 3.5 in Howland (1971) A(z)−1 has a simple pole at
z = z0 if and only if (i) A(1) (z0 ) is injective on ker A(z0 ) and (ii) ran A(z0 ) ∩
A(1) (z0 ) ker A(z0 ) = {0}. Since A(z) = idH −C(z) and C(z) is a compact
family of operators, A(z) is an index-zero Fredholm family of operators. This
implies that dim(ker A(z0 )) = dim(ran(idH −P1 )). From this, one can easily
show that the invertibility of (idH −P1 )A(1) (z0 )|ker A(z0 ) is another equivalent
condition to (i) and (ii).
Lemma 8.3. Let d = 1 in Lemma 8.1. Then,
(a) G(z0 ) = A(z0 ) − (idH −P1 )A(1) (z0 ),
(b) ran G(z0 ) |ker A(z0 ) = ran(idH −P1 ),
Proof. We may easily deduce (b) from (a) and Lemma 8.2, so we only show
(a). We know that G(z) is holomorphic at z = z0 from Lemma 8.1, so
(1)

G(z) = G(z0 ) − G (z0 )(z0 − z) +

∞
X
G(j) (z0 )
j=2

33

j!

(z − z0 )j .

Therefore,
A(z) = [P1 + (idH −P1 )(z0 − z)][G(z0 ) − G(1) (z0 )(z0 − z) + H1 (z)],

(8.1)

where H1 (z) denotes the remainder of the Taylor series of G(z). Since A(z)
is holomorphic at z = z0 , we have
A(z) = A(z0 ) − A(1) (z0 )(z0 − z) + H2 (z),

(8.2)

where H2 (z) is the remainder of the Taylor series of A(z). Collecting terms
associated with the same powers from (8.1) and (8.2) we obtain
P1 G(z0 ) = A(z0 )

(8.3)

(idH −P1 )G(z0 ) = −A(1) (z0 ) + P1 G(1) (z0 ).

(8.4)

The left-hand side of the second equation is invariant under (idH −P1 ), implying that
(idH −P1 )G(z0 ) = −(idH −P1 )A(1) (z0 ).

(8.5)

Combining (8.3) and (8.5), we obtain (a).

8.2

Log-density estimation

In Section 6, the estimated log-density ĝt is obtained from the following
procedure. The reader is referred to Loader (1996, 2006) for more details.
Given
P survey responses X1 , . . . , Xn with design weights w1 , . . . , wn such
that i=1 wi = n, consider the weighted log-likelihood
Z

n
X
L(ft ) =
wi log(ft (Xi )) − n
ft (u)du − 1 .
i=1

Let K be the support of ft . Under some local smoothness assumptions, we
can consider a localized version of the log-likelihood and log ft (u) can be
locally approximated by a polynomial function, so follows.
Lp (ft )(x) =

n
X
i=1


wi W

Xi − x
h


Q(Xi − x; αt )


Z
−n

W

34

u−x
h


exp(Q(u − x; αt ))du (8.6)

where W is a suitable kernel function, h is a bandwidth which assumed to
be fixed, and Q(u; αt ) is polynomial in u with coefficients αt .
(1 − |u|3 )3 and Q(u; α) = α0,t + α1,t u. For fixed x ∈ K, let
I set W(u) = 70
81
(α̂0,t , α̂1,t ) be the maximizer of (8.6). Then the local likelihood log-density
estimate is given by
ĝt (x) = α̂0,t
The procedure is repeated for a fine grid of points, and then ĝt may be
obtained from an interpolation method described in (Loader, 2006, Chapter
12). Needless to say, the above estimation procedure depends on h, a fixed
−1/5
bandwidth parameter. Let bt := [qt (99) − qt (1)] nt
where qt (a) is the a-th
percentile of observations at time t and nt is the number of observations.
Within the range [0.75 bt , 1.25 bt ], bandwidth h is set to the minimizer of a
generalized BIC criterion.7

7

A generalized version of AIC as a diagnostic for the local likelihood method is given
in Loader (2006). A suitable generalization of BIC is obtained by changing the penalty
term in an obvious way.

35

