arXiv:1801.06727v1 [q-fin.ST] 20 Jan 2018

A Second Order Cumulant Spectrum Based
Test for Strict Stationarity
Douglas Patterson, 1
Department of Finance, Virginia Tech,
Melvin Hinich †,2
Denisa Roberts 3
Amazon
January 23, 2018

Abstract
This article develops a statistical test for the null hypothesis of strict stationarity
of a discrete time stochastic process. When the null hypothesis is true, the second
order cumulant spectrum is zero at all the discrete Fourier frequency pairs present
in the principal domain of the cumulant spectrum. The test uses a frame (window)
averaged sample estimate of the second order cumulant spectrum to build a test
statistic that has an asymptotic complex standard normal distribution. We derive
the test statistic, study the size and power properties of the test, and demonstrate its
implementation with intraday stock market return data. The test has conservative
size properties and good power to detect varying variance and unit root in the presence
of varying variance.

Keywords: Discrete Time, Second Order Cumulant Spectrum, Stochastic Process, Time
Varying Variance, Unit Root, High Frequency Stock Returns
1

Douglas Patterson is Professor of Finance in the Department of Finance, Pamplin College of Business,
Virginia Tech, Blacksburg, VA 24061 (email: amex@vt.edu).
2
Deceased September 6, 2010.
3
Denisa Roberts is Machine Learning Scientist at Amazon, NY, NY (email: rdenisa@amazon.com).
This work was completed before joining Amazon.

1

1.

Introduction.

A stochastic process is called strictly stationary if the joint distribution of x (t1 + τ ),
x (t2 + τ ), . . ., x (tn + τ ) is independent of τ for all time points {t1 , t2 , . . . , tn }. In other
words, the time series has the same joint distribution across all time points t. The hypothesis of stationarity plays a significant role in time series analysis and is frequently used as an
assumption in both academic research and practical applications. A unit root is one type of
violation of stationarity in the strict sense. Kwiatkowski et al. (1992) developed a test for
unit root (or level stationarity) commonly mentioned in the stationarity related literature,
hereafter referred to as KPSS. Because of its popularity, we use KPSS as a starting point
of comparison for our work.
Weak stationarity, commonly used as an assumption in financial time series forecasting
models, assumes that the first two moments of the time series are constant through time
(under suitable mixing conditions that the series is summable, and the first two moments
exist Brillinger (2001)). Hence a process is weakly stationary (or stationary in the wide
sense), if its expected value E{x (tm )} is constant and its autocovariance cx (tm ) depends
only on the difference (tn+m − tn ), cx (tm ) = cx (tn+m − tn ) (Brillinger (2001)). A weakly
stationary Gaussian process is also stationary in the strict sense Brillinger (2001). A
Gaussian process is a sequence of independent, identically normally distributed random
variables {ut }.
There are two paradigms for analyzing time series: the time domain and the frequency
domain. In the time domain, one considers the observed data {x (t)} directly and typically
makes conjectures about its moments. In the frequency domain, one decomposes the time
series into underlying frequencies and makes conjectures about spectra and cumulant spectra. The frequency domain offers its own advantages to decompose time series into trends

2

and work with higher order moments more easily. In this article we develop a statistical
test of the hypothesis that an observed stochastic process {x (t)} is stationary in the strict
sense. For convenience, we will refer to our test as PHR. PHR uses sample estimates of
the second order cumulant spectrum to test the null hypothesis of strict stationarity.
1.1.

Related Work.

The PHR second order cumulant spectrum based test of strict stationarity is different from
widely used tests for structural change of a linear time series model. For example, Farley
and Hinich (1970) proposed a control chart to detect small shifts in the mean of a stochastic
time series modeled as a white noise process plus a mean parameter that can be constant
or shifting. Their proposed chart draws no conclusions on the stationarity of the process.
Estimation and inference about the change in the mean of a stationary random process
was amply examined in literature. Many articles studies the estimation and inference about
the change in the mean of a random process. Hawkins (1977) proposed a likelihood ratio
based procedure to detect a change in the mean while assuming constant variance. In
a different approach, Dickinson et al. (2014) developed CUSUM and EWMA charts that
detect changes in probability distributions via their scale parameter. Also, tests of a unit
root with or without drift against an alternative hypothesis of trend stationarity, can be
found, for example, in Dickey and Fuller (1979) and Phillips and Perron (1988).
The space of nonstrictly stationary processes is vast and, therefore, in this manuscript
we focus the evaluation of the power of the PHR test on a few alternative hypotheses. We
do not study PHR’s power against the alternative of a process with periodicities in the
mean and random variation in the periodic structure, as in Hinich and Wild (2001). In the
space of tests for strict stationarity in the time domain, Lima and Neri (2013) proposed a
test for strict stationarity and compared its power to the KPSS test of level stationarity.
3

Similarly to Lima and Neri (2013) we compare the power of the PHR test to the power of
the KPSS test. We calculate the power for the alternative hypotheses of unit root, varying
variance, and their combination. Our PHR test excels in detecting stationarity violations
due to second moment.
We present the properties of the second order cumulant spectrum in Section 2, the test
development in Section 3, empirical size and power evaluation via Monte Carlo simulations
and comparison with the KPSS test under various alternative hypotheses in Section 4. The
test is demonstrated in 5 using high frequency (trade by trade) stock returns for the Ford
company .
2.

The Second Order Cumulant Spectrum of a Strictly Stationary Process.

In the frequency domain, the second order cumulant spectrum of a stochastic process
{x(tn )} with tn = nτ
K (f1 , f2 ) =

∞
X

∞
X

E [x (tn1 ) x (tn2 )] exp [−2π (f1 tn1 + f2 tn2 )] ,

(1)

n1 =−∞ n1 =−∞

where fi are discrete Fourier frequencies.
The cumulant spectrum K is defined for frequency values in the square {−f0 ≤ f1 ≤
f0 , −f0 ≤ f2 ≤ f0 }, where f0 = 1/ (2τ ) is called the Nyquist frequency. Fourier frequencies
display the following properties and symmetries:
• Frequencies are between zero and one, normalized by Nyquist frequency. The Nyquist
frequency f0 is equal to the inverse of twice the sampling rate of a discrete time
process. The time series must be sampled at specific time intervals. We assume
for simplicity that the sampling interval length denoted τ is one. When τ = 1 the
Nyquist frequency is f0 = 0.5.
4

• The discrete Fourier transform of a sample of length T at frequency f , X (f ) =
PT −1
t=0 x (t) exp{−i2πf t} is equal to the conjugate of the Fourier transform at −f ,
X(f ) = X ∗ (−f ).
• The Fourier transform is periodic at intervals of the sampling frequency for discrete
time processes. Higher frequencies fold down into the (0, f0 ) interval.
Because of these symmetries, we focus the study of the second order cumulant spectrum
in the principal domain only. The principal domain contains a subset of frequencies that
are not aliased or indistinguishable from each other (Bloomfield (2004)), also called fundamental frequencies. The principal domain forms the triangle {0 < f1 ≤ f0 , −f1 < f2 ≤ f1 }.
In general, to cross from the time domain into the frequency domain, we use the discrete
Fourier transform of the time series realization at a series of frequencies in the principal
domain. From Hinich (1994), if the time series {x (t)} is strictly stationary then the
second order cumulant spectrum for a pair of frequencies in the principal domain (f1 , f2 ),
K (f1 , f2 )is equal to

 2S (f ) + O (1) if f + f = 0 mod 1
1
1
2
K (f1 , f2 ) = E[X(f1 )X(f2 )] =
 O (1)
otherwise,

(2)

where X(fi ) denote the discrete Fourier transform of the time series at the frequencies
included in the principal domain. The line f2 = −f1 is not in the principal domain (Hinich
(1994) ). From Equation 2, K (f1 , −f1 ) = δ (f1 ) S (f1 ), where δ (f1 ) is the Dirac delta
P
function, and S (f1 ) = ∞
m=−∞ cx (tm ) exp (−i2πf1 tm ) is the spectrum at frequency f1 and
cx (tm ) is the autocovariance of the process.
In summary, if a time series is strictly stationary, then its second order cumulant spectrum

5

K (f1 , f2 ) = 0

(3)

for all frequency pairs (f1 , f2 ) in the principal domain Hinich (1994).
3.

Derivation of The Test Statistic for the Second Order Cumulant Spectrum Test
of Strict Stationarity.

We now use the second order cumulant spectrum to construct the test for the null hypothesis that a stochastic process {x(tn )} is strictly stationary. The test uses frame averaged
cumulant spectra at fundamental frequency pairs and has an asymptotic complex normal
distribution. To get frame averaged estimates for the cumulant spectrum K(f1 , f2 ) and for
the spectrum S(f ), where (f1 , f2 ) are frequencies in the principal domain, we partition the
time series into frames of equal length, calculate the cumulant spectrum and the spectrum
for each frame and then average over all the frames. A time series of length T is partitioned
into P = [T /L] complete nonoverlapping frames of length L. We omit the last frame if it
has less than L observations. Then the frame averaged cumulant spectrum estimate is
P
1 X
Kp (f1 , f2 ) ,
K̂ (f1 , f2 ) =
P p=1

(4)

We can also express the principal domain for K in terms of integers (k1 , k2 ), with
f = k/L, and 0 < k1 ≤ L/2, −k1 ≤ k2 ≤ k1 . The steps to get the frame averaged estimates
of the cumulant spectrum, K̂, and of the spectrum, Ŝ, are calculated following the steps:
• Use a reasonable length time series, ideally of no less than T = 150 observations.
• Split the time series into P = [T /L] complete consecutive nonoverlapping frames or
windows of L elements, where L is sufficiently large. Let the p frame of length L be
6

{xp (1), xp (2), ..., xp (L)}.
• Calculate the discrete Fourier transform (DFT) for the p frame of length L as

Xp (k) =

L−1
X

x [t + (pL)] exp{−i2πfk [t + (pL)]},

(5)

t=0

where fk = Lk .
• Compute the spectrum at each frequency in the principal domain of K for each frame
p

Sp (k) =

1
1
|Xp (k) |2 = Xp (k) Xp (−k)
L
L

(6)

• Compute the second order cumulant spectrum for each frame p

Kp (k1 , k2 ) =

1
Xp (k1 ) Xp (k2 )
L

(7)

for pairs of frequencies defined by 0 < k1 ≤ L/2 and −k1 < k2 ≤ k1 .
• Get the frame averaged Ŝ and K̂ over the P frames for each k and (k1 , k2 ) pair
respectively.
From the theory in Hinich (1994) and in Brillinger (2001), if L and P are sufficiently
large, the expected value of the frame averaged spectrum is equal to its theoretical value
up to O(L−1 ) ,
h

i


E Ŝ (fk ) = S (fk ) + O L−1 .

7

(8)

Similarly to Equation 8, under strict stationarity, from Hinich (1994)
h
i

E K̂ (k1 , k2 ) = K (k1 , k2 ) + O L−1 .

(9)

We stated in Equation 3 that the second order cumulant spectrum K (k1 , k2 ) for strictly
stationary time series is zero for pairs of frequencies in the frequency domain. Furthermore,
the variance of K̂, considering P frames, is equal to P −1 S (k1 ) S (k2 ) as (L, P → ∞) from
Hinich (1994).
Next, we define the normalized second order cumulant spectrum as
Γ̂ (k1 , k2 ) = p

K̂ (k1 , k2 )
S (fk1 ) S (fk2 )

.

(10)

We make a further simplifying assumption used in Hinich and Rothman (1998) that the
time series has been prewhitened, hence the theoretical spectrum S is assumed constant
across frequencies and, without loss of generality, equal to one. With this assumption, the
√
√
standard error of K̂ becomes 1/ P . Then the estimate for Γ, Γ̂ simplifies to P K̂.
Next we construct the complex valued pivotal quantity Y as
√
Ŷ (k1 , k2 ) =

h

i

2P Γ̂ (k1 , k2 ) − Γ (k1 , k2 ) .

(11)

From the central limit theorem, the pivotal quantity Y has an asymptotic complex
normal distribution. Hence, the real value <(Y ) and the imaginary value =(Y ) are bivariate standard normal random variables asymptotically. Under the alternative hypothesis
of nonstationarity, the theoretical distribution depends on the type of strict stationarity
violation.
Notice that we calculate a Y value for each pair of (k1 , k2 ) in the principal domain. The
number of calculated cumulants spectra equals the number of pairs (k1 , k2 ) in the principal
8

domain, which depends on the frame length L. We assume a sufficiently large frame length
L so that the asymptotic theory holds.
We can test the normality of the stacked vector of real and imaginary values using the
Kolmogorov Smirnoff (KS) one sample test. Let Femp (x) denote the empirical cumulative
distribution of the vector of stacked (<(Y ), =(Y )). We calculate Femp (x) as the fraction
of values that are less than a given value x in the unit interval. In our case the vector
x represent the real and imaginary parts for the PHR test statistic and the cumulative
distribution of interest is that of a standard normal variable. The KS statistic is D̂ =
sup0<x≤1 |Femp (x) − x| (Press et al. (1992)).
4.

Empirical Evaluation of Size and Power of the PHR Test. Comparison to the
KPSS Test.

Next we proceed to evaluate the empirical size and power of PHR to detect departures
from the hypothesis of strict stationarity. We compare PHR to the KPSS test for level
stationarity implemented in the Trapletti and Hornik (2015) tseries R package. Developed
in Kwiatkowski et al. (1992), the KPSS test tests the null hypothesis of first order (level
or absence of unit root) stationarity using the test statistic
!2
j
T
1 X X
(x(t) − x̄(t)) ,
KP SS = 2 2
σ̂ T j=1 t=1

(12)

where xt is the time series, with t = 1 . . . T , where T is the length of the times series, x̄(t)
is the mean of the time series, and σ̂ is a consistent estimator of the long run variance.
Lima and Neri (2013) studied the size and power of KPSS to detect various departures
from the strict stationarity hypothesis. Cavaliere and Taylor (2005) evaluated KPSS performance to detect unit root in the presence of varying second moment. In the following
9

sections we replicate several of the data generating processes in Lima and Neri (2013) and
Cavaliere and Taylor (2005) and compare PHR with KPSS in terms of power to detect
departures from strict stationarity.
We first generated strictly stationary data and evaluated how closely the empirical test
size tracked the nominal test size of five percent and one percent respectively. PHR has
conservative size properties. We then evaluated the power of the test. The null hypothesis
of strict stationarity can be violated in various ways, for example unit root, seasonal root,
time varying second moment, and long memory. We expect that under large sample size
conditions PHR detects a wide variety of violations. We limit ourselves to evaluating the
power of PHR to detect unit root, varying second moment, and the combination of the
two. The power of the test to detect time varying second moment in almost all considered
cases as well as the power to detect unit root in certain cases surpasses the KPSS test.
We built an R R Core Team (2015) package to implement PHR (available upon request).
In the R package we use the base R implementation of the one sample KS test to test that
the vector of real and imaginary parts of the Ŷ test statistic in Equation 11 likely comes
from the standard normal distribution. The KS test assumes a continuous distribution
of the tested sample and a large sample size. The approximate p value of the KS test is
calculated using the technique described in Wang et al. (2003). The p value is approximate
for two reasons. First, the test assume that the sample is large enough for the asymptotic
theory to hold, and second, the test assumes that sample values have no ties (but ties are
possible due to the discrete nature of the values).
Each data generation process in the size and power study was replicated10,000 times
and PHR and KPSS size and power were calculated as the percentage of times each test
rejected the null of strict stationarity under each data generation scenario. For reproducible
results, the random number generator seed was set and reused for each evaluated scenario.
10

We performed all the simulations in the R open source software and we created all the
tables using the stargazer Hlavac (2015) R package. Simulation code is available upon
request.
Our PHR is built on the assumptions that time series are free from outliers, have
finite moments, do not trend, are prewhitened and have zero means. We evaluated the
performance of both PHR and KPSS test to detect a unit root when data are generated from
a normal distribution or a t distribution with three, five and ten degrees of freedom. We also
evaluated the performance of the tests after first trimming (drop top and bottom 0.5, one
or two percent of the data respectively), demeaning (subtract the mean), and detrending
the data. The Borchers (2015) pracma package detrended the data by estimating and
removing a least squares fit. Demeaning, detrending, and trimming are commonly used in
fast Fourier transform applications.
Lastly, we generated data with time varying second moment, as well as a combination
of varying second moment and unit root, and evaluated the power of PHR and KPSS. For
the sake of brevity, we only included partial and essential tables in this manuscript. All
simulation results are available upon request.
4.1.

Empirical Test Size.

First, we generated a white noise process and evaluated PHR empirical size sensitivity to
sample size and frame length L. Simulation results in Table 1 show that it is best to have
sample sizes above 500, at least 50 frames and a frame length of at least 10. A rule of
thumb is to set L and P equal to the square root of sample size T .
Next, we generated strictly stationary data in order to evaluate the empirical sizes of
PHR and KPSS tests. We generated samples from a t distribution with ∞ (normal), three,
five, ten and 15 degrees of freedom. Furthermore, we included an autocorrelation structure
11

Table 1: Empirical Test Size for Various Sample Sizes T and Frame Lengths L (corresponding number of frames P=[T/L]). Nominal Test Size of 0.05.
L
10

T.250

T.500

0.062

0.059

T.1000

T.5000

0.062

0.058

20

0.062

0.058

0.060

0.058

50

0.195

0.084

0.063

0.057

100

1

0.716

0.183

0.059

200

1

1

1

0.122

of +0.3 or −0.3 to evaluate the robustness of the test to the assumption of prewhitening.
We also evaluated PHR and KPSS empirical test size sensitivity to sample size, frame
length, demeaning, detrending and trimming. We computed and compared the empirical
test size of PHR and KPSS corresponding to nominal test sizes of five and one percent.
4.1.1

Data Generation Processes for Size of Test Simulations.

The stationary data generation process is

x(t) = u(t).

(13)

u(t) = ρu(t − 1) + v(t),

(14)

The innovations u(t) follow the process

where v(t) are independently and identically distributed with one of the t distributions.
When ρ is zero, u(t) in Equation 14 process reduces to just the v(t) component. If ρ is
different from zero the time series has short term memory.

12

4.1.2

Simulation Results for Empirical Test Size.

Table 2 shows size of test results for a white noise process with zero autocorrelation in the
error term. We expect the empirical size of the test to be close to the nominal size of the
test. Each entry in the table represents the percentage of times PHR and KPSS rejected
the null hypothesis of strict stationarity when the null was in fact true, in 10,000 runs. We
also compared the size of PHR and KPSS after demeaning, trimming and detrending the
data.
When we do not preprocess the data, PHR showed slightly higher empirical size than
KPSS, too frequently rejecting the null hypothesis. For a larger sample size of T = 5000
the empirical size of PHR is closer to the nominal size of the test. When we trim the data
at different rates, we noticed that PHR’s empirical test size decreased while KPSS size
remained unchanged. PHR is sensitive to outliers. On the other hand KPSS is extremely
sensitive to detrending, achieving zero empirical test size in all cases where we included
detrending. Neither PHR nor KPSS reacted excessively to demeaning. PHR’s empirical
size varies with the sample size T , while KPSS’s empirical size is less dependent on sample
size.
The empirical sizes of PHR and KPSS are consistent across t distributions with three,
five, ten and 15 degrees of freedom and zero autocorrelation. We only included here the
table for five degrees of freedom but results are similar across degrees of freedom. PHR is
sensitive to fat tails and fares better with appropriate trimming. Table 3 illustrates that a
one percent trimming in the tails tames the empirical Type I error rate for PHR. Results
(not included) for two percent trimming were similar. The KPSS is less sensitive to fat
tails but reacts strongly to detrending, as observed in previous tables.
Next, we evaluated PHR and KPSS empirical size sensitivity to a ρ = 0.3 autocorre-

13

Table 2: Empirical Size of Test Comparisons for Normal Data with Zero Autocorrelation
in Error Term(ρ = 0) and Nominal Size α = 0.05.
itrim

idetrend

idemean

0

0

0

KPSS.250
0.049

PHR.250
0.062

KPSS.500

PHR.500

0.052

0.084

KPSS.1000

PHR.1000

0.050

0.063

KPSS.5000
0.052

PHR.5000
0.059

1

0

0

0.050

0.051

0.051

0.070

0.051

0.044

0.051

0.043

0

1

0

0

0.062

0

0.085

0

0.064

0

0.059

1

1

0

0

0.049

0

0.069

0

0.046

0

0.042

0

0

1

0.049

0.062

0.052

0.084

0.050

0.063

0.052

0.059

1

0

1

0.050

0.051

0.051

0.070

0.051

0.044

0.051

0.043

0

1

1

0

0.062

0

0.085

0

0.064

0

0.059

1

1

1

0

0.049

0

0.069

0

0.046

0

0.042

Table 3: Empirical Size of Test Comparisons for t5 Data with Zero Autocorrelation in
Error Term(ρ = 0) and Nominal Size α = 0.05.
itrim

idetrend

idemean

0

0

0

KPSS.250
0.046

PHR.250
0.152

KPSS.500

PHR.500

0.048

0.266

KPSS.1000

PHR.1000

0.048

0.247

KPSS.5000
0.048

PHR.5000
0.280

1

0

0

0.048

0.066

0.048

0.106

0.047

0.068

0.048

0.064

0

1

0

0

0.152

0

0.268

0

0.246

0

0.280

1

1

0

0.001

0.067

0.001

0.105

0.001

0.066

0.001

0.064

0

0

1

0.046

0.152

0.048

0.266

0.048

0.247

0.048

0.280

1

0

1

0.048

0.066

0.048

0.106

0.047

0.068

0.048

0.064

0

1

1

0

0.152

0

0.268

0

0.246

0

0.280

1

1

1

0.001

0.067

0.001

0.105

0.001

0.066

0.001

0.064

14

lation in the error term. Prewhitening of data is recommended before applying PHR. For
normal data the empirical size for both PHR and KPSS fared similarly to the case with
zero autocorrelation. When data were generated from a t distribution with three, five, ten
and 15 degrees of freedom we notice the necessity for trimming. We also notice the strong
KPSS reaction to detrending. Notice in Table 4 how when we have large sample size of
T = 5000, as required by the asymptotic properties of PHR, PHR’s empirical size is very
close to nominal size and stable over demeaning, detrending and trimming scenarios while
KPSS’s empirical size varies.
Table 4: Empirical Size of Test Comparisons for Normal Data with Autocorrelation in
Error Term ρ = 0.3 and Nominal Size of Test α = 0.05.
itrim

idetrend

idemean

KPSS.250

PHR.250

KPSS.500

PHR.500

KPSS.1000

PHR.1000

KPSS.5000

PHR.5000

0

0

0

0.048

0.062

0.050

0.084

0.047

0.061

0.047

0.058

1

0

0

0.045

0.055

0.051

0.067

0.048

0.044

0.049

0.040

0

1

0

0

0.061

0

0.084

0

0.061

0

0.058

1

1

0

0

0.055

0

0.068

0

0.045

0

0.040

0

0

1

0.048

0.062

0.050

0.084

0.047

0.061

0.047

0.058

1

0

1

0.045

0.055

0.051

0.067

0.048

0.044

0.049

0.040

0

1

1

0

0.061

0

0.084

0

0.061

0

0.058

1

1

1

0

0.055

0

0.068

0

0.045

0

0.040

When the data generation process included negative −0.3 autocorrelation in the error
term, the empirical size simulation results were similar to Table 4 results.
4.2.

Power of Test to Detect Unit Root.

Next we evaluate the power of PHR to detect unit root and compare its power to the power
of KPSS test.

15

4.2.1

Data Generation Processes With Unit Root.

We replicated the unit root data generation process from Lima and Neri (2013). The
generated data included a unit root component y(t) and a stationary component u(t). The
tuning parameter λ moderates the strength of the unit root process. Then the time series
x(t) is

x(t) = λy(t) + u(t),
with y(t) =

Pt

j=1

(15)

v(j).

The innovations u(t) follow the process in Equation 14. The y(t) process is a random
walk with v(t) i.i.d from the same distribution as u(t) but independently from u(t). We
compared the power of the PHR to the power of KPSS test to detect a unit root in processes
with relatively low λ of 0.1 and 0.3. A larger λ led to 100% power for both tests, as observed
in additional simulation studies not included here.
4.2.2

Simulation Results Unit Root.

For a process with high lambda (strong presence of unit root), the estimated cumulant
spectrum values depart strongly from zero, the spectrum estimate explodes at frequencies
near zero and the distribution of the resulting test statistic has much heavier tails than a
standard normal. Under the null hypothesis of strict stationarity, the test statistic follows
an asymptotic normal distribution.
The power of PHR is greater for stronger unit root presence (larger λ value), having
difficulty with smaller λ values. We also evaluated the impact of trimming, demeaning and
detrending on the power of the test.
Table 5 illustrates the good performance of PHR in large sample sizes for normal data.

16

Table 5: Power of Test Comparisons for Normal Data with a Unit Root (λ = 0.3), Zero
Autocorrelation in Error Term(ρ = 0), and Nominal Size of Test α = 0.05.
itrim

idetrend

idemean

KPSS.5000

0

0

0

1

PHR.5000
0.997

KPSS.250

PHR.250

0.963

0.080

KPSS.1000
0.997

PHR.1000
0.420

1

0

0

1

0.996

0.956

0.072

0.996

0.406

0

1

0

0.998

0.997

0.620

0.078

0.936

0.393

1

1

0

0.998

0.997

0.597

0.072

0.928

0.376

0

0

1

1

0.997

0.963

0.080

0.997

0.420

1

0

1

1

0.996

0.956

0.072

0.996

0.406

0

1

1

0.998

0.997

0.620

0.078

0.936

0.393

1

1

1

0.998

0.997

0.597

0.072

0.928

0.376

The power is very close to one for both PHR and KPSS in sample sizes of 5000. The power
of PHR deteriorates at lower sample sizes while KPSS holds its ground. Power deteriorates
for both KPSS and PHR with detrending and trimming.
Table 6: Power of Test Comparisons for t5 Data with Unit Root (λ = 0.3), Zero Autocorrelation in Error Term(ρ = 0), and Nominal Size of Test α = 0.05.
itrim

idetrend

idemean

0

0

0

KPSS.250
0.968

PHR.250
0.148

KPSS.1000
0.998

PHR.1000
0.573

1

0

0

0.965

0.106

0.997

0.534

0

1

0

0.628

0.148

0.939

0.550

1

1

0

0.623

0.098

0.935

0.493

0

0

1

0.968

0.148

0.998

0.573

1

0

1

0.965

0.106

0.997

0.534

0

1

1

0.628

0.148

0.939

0.550

1

1

1

0.623

0.098

0.935

0.493

In Table 6, we generated data from a t distribution with five degrees of freedom, no
17

correlation in error term and a unit root with λ = 0.3. PHR performed similarly to the
normal case, displaying the already established sensitivity to sample size and fat tails. In
Table 6 we see slightly higher power to detect the unit root process in smaller sample sizes
as a reaction to fat tails. It is also interesting to notice how the power of KPSS to detect
unit root decreases significantly with small sample sizes such as T = 250 when we apply
higher trimming, and markedly decreases with detrending.
PHR and KPSS performed similarly when we included autocorrelation in the error term,
as well as when λ decreased to 0.1 (weaker presence of the unit root process). Both PHR
and KPSS performed well for larger sample sizes.
4.3.

Power of the Test to Detect Time Varying Second Moment.

Another alternative to the null hypothesis of stationarity is a varying second moment.
The authors of Cavaliere and Taylor (2005) evaluated the power of the KPSS test and
two other tests to detect a unit root in the presence of varying variance. We employed
their data generating processes in a similar but not identical fashion and evaluated the
power of PHR to detect the alternative hypothesis of time varying second moment for
three different patterns of nonconstant variance. The null hypothesis is the hypothesis of
strict stationarity of the time series, when the second order cumulant spectrum is equal to
zero. The alternative hypothesis is nonstationarity due to varying second moment. Hence,
our set up is different from Cavaliere and Taylor (2005), and similar to Lima and Neri
(2013).
Lima and Neri (2013) also evaluated several tests for strict stationarity and their power
against the alternative of unconditional heteroskedasticity. Lima and Neri (2013) stated
that tests for level stationarity (or unit root) such as the KPSS test (and its robust version
iKPSS) have low power against the alternative of time varying variance. Furthermore they
18

added that KPSS is designed based on the fluctuation of the data around the sample mean
and it has low power to detect how the scale of the distribution varies over time. A process
with time varying scale has varying second moment and hence is not strictly stationary.
Therefore we are interested in the power of PHR to detect a varying second moment. Lima
and Neri (2013) only discussed one pattern of heteroskedasticity or varying second moment
while we considered the three patterns suggested in Cavaliere and Taylor (2005).
4.3.1

Data Generation Processes With Time Varying Second Moment.

Next, we describe the data generating processes with time varying variance patterns and
the simulation results for the power of PHR and KPSS test. The first data generating
process (DGP1) in Equation 16 is for a level stationary (no unit root) process with variance
σt . Under the null hypothesis of time constant variance, σt2 = 1. Under the alternative
hypothesis of varying variance, σt2 follows one of the three processes in Equation 18. DGP1
is

x(t) = σt u(t),

(16)

with t = 1 . . . T and u(t) a white noise.
The second data generating process (DGP2) in Equation 17 generates a process with
unit root and variance σt2 . Similarly to DGP1, σt2 = 1 under the null hypothesis of constant
variance. Under the alternative hypothesis of varying variance, σt2 follows one of the three
processes in Equation 18. DGP2 is

x(t) = x(t − 1) + σt u(t),

(17)

where u(t) are white noise and x(0) = 0. When we evaluated the power of PHR and
KPSS tests to detect violations in DGP2, we effectively study two confounding sources
19

of nonstationarity, both unit root and varying second moment. In this case, the tests do
not indicate which one of the sources of nonstationarity caused the rejection of the null
hypothesis.
For clarity, we only allow σt2 to vary according to the three patterns described below,
while the initial paper by Cavaliere and Taylor (2005) allowed the same patterns in the
irregular component u(t) as well. Hence, our simulation results will naturally differ from
theirs. Equation 18 defines the three variance patterns as: single and multiple breaks in
variance, smooth transition variance breaks and linear trending variance. Notice how the
variance is constant up to a certain point in time, m, after which it evolves according to
one of the three functions.
• Single Break
σt2 =



σ0 2 + (σ1 2 − σ0 2 ) , t ≥ m

σ0 2 ,

t<m

where σ0 = 1 and σ1 = cσ0 .
• Smooth Transition

σt2 = σ0 2 + σ1 2 − σ0 2 W (t) ,
where W (t) = (1 + exp (−γ (t − m)))−1 , σ0 = 1 and σ1 = cσ0 .
• Piecewise Linear Trend


σ0 2 + (σ1 2 − σ0 2 ) (t − m) (1 − m)−1 , t ≥ m
2
σt =

σ0 2 ,
t<m
where σ0 = 1 and σ1 = cσ0 .
20

(18)

The tuning parameters c and m dictate the scale and time of the change in variance
under the alternative hypothesis of varying variance. We considered c values of {0.25, 4},
corresponding to a decreased or increased variance, respectively. For the single break
and smooth transition cases, m takes the values {0.1, 0.5, 0.9}. For a sample size T of
1000 observations, m = 0.1 indicates that after 100 observations we should see a variance
change. The γ = 10 parameter in Equation 18 gives the speed or smoothness of transition
from constant variance for the smooth transition variance pattern.
For the piecewise linear trend in Equation 19 m takes the values {0, 0.5, 0.9}. In this
case, when m > 0 the variance is fixed and equal to one until time mT and then it changes
to an increasing (when c = 4) or decreasing (when c = 0.25) linear trend. Higher m dictates
a later linear trend with a higher slope. When m = 0 the variance process follows a linear
trend starting at t = 1.
Figure 1 displays the patterns in σ for each of single break, smooth transition and
piecewise linear trend. In each of the four plots, variance increases, corresponding to the
shift parameter c = 4 and the time series sample is 1000 observations long. For the single
break and smooth transition cases, m = 0.5 indicates a variance shift at the time series
midpoint (after 500 observations). For the piecewise linear trend, the plot displays both
the trend for m = 0.5 and for m = 0.
Figure 2 illustrates both DGP1 and DGP2 processes with time varying variance according to the three patterns.
4.3.2

Simulation Results Power to Detect Time Varying Second Moment.

We present next Monte Carlo simulation results to compare the power of PHR and KPSS
test to detect time varying variance. Table 7 contains the power of the two tests to detect
piecewise linear trend in second moment (increase if c = 4 and decrease if c = 0.25) that
21

Piecewise Linear Trend

4.0

4.0

m = 0.5

Single Break

m = 0.5

2.5

sigma

3.0

3.5

m=0

1.0

1.5

2.0

2.5
1.0

1.5

2.0

sigma

3.0

3.5

Smooth Transition

0

200

400

600

800

1000

time index

0

200

400

600

800

1000

time index

Figure 1: The Three Functions for a Variance Increase and T = 1000. Left: Single Break
and Smooth Transition. Right: Piecewise Linear Trend.
starts at different points in the sample (early if m = 0, half way if m = 0.5 and late if
m = 0.9), for different sample sizes, T = {250, 1000, 5000}, and at five percent nominal
size of the test. PHR outperformed KPSS. Both PHR and KPSS disappointed for a few
combinations of parameters m and c. These cases are for DGP1 process (no unit root) and
c = 4 (increase in variance). Even smaller sample sizes fared well for PHR power. For large
sample sizes, the two tests agreed.
For T = 1000, PHR achieved 100% power in almost all cases, except when detecting a
variance decrease that started at a later stage. The KPSS test had close to 99% power for
DGP2, which is a process that includes both a unit root and time varying variance. For
the DGP1 process, when we have time varying variance and no unit root, KPSS had power
close to 99%. The power patterns for both PHR and KPSS are almost identical for the
T = 5000 sample size. PHR power was slightly less but still above 95% when the sample
size was T = 250. However, when the sample size was T = 250 with varying variance
22

Unit Root and Single Break

5

dgp2

0

−5

−3

−2

0

−1

dgp1

1

2

10

Level Stationary and Single Break in Variance

0

50

100

150

200

250

0

50

100

150

200

250

time points

Level Stationary and Smooth Transition Variance

Unit Root and Smooth Transition Variance

5

dgp2

0
−3

−5

−2

0

−1

dgp1

1

10

2

15

time points

0

50

100

150

200

250

0

50

100

150

200

250

time points

Level Stationary and Piecewise Linear Variance

Unit Root and Piecewise Linear Variance

0

dgp2

−2

0

−4

−1

−6

−2

dgp1

2

1

4

2

6

time points

0

50

100

150

200

250

0

time points

50

100

150

200

250

time points

Figure 2: Time Series with Time Varying Variance

23

Table 7: Power of Tests to Detect Piecewise Linear Trend at Nominal Test Size α = 0.05.
dgp1 2

m

c

PHR.250

KPSS.250

PHR.1000

KPSS.1000

PHR.5000

KPSS.5000

1

0

4

0.090

0.054

0.742

0.056

1

0.055

2

0

4

0.995

0.982

1

0.998

1

1

1

0

0.25

0.094

0.050

0.762

0.051

1

0.054

2

0

0.25

0.996

0.982

1

0.997

1

1

1

0.5

4

0.273

0.098

1

0.099

1

0.100

2

0.5

4

0.987

0.966

1

0.995

1

1

1

0.5

0.25

0.078

0.036

0.378

0.036

0.956

0.039

2

0.5

0.25

0.997

0.984

1

0.998

1

1

1

0.9

4

0.447

0.114

1

0.120

1

0.125

2

0.9

4

0.990

0.974

1

0.996

1

1

1

0.9

0.25

0.066

0.042

0.084

0.042

0.105

0.047

2

0.9

0.25

0.998

0.984

1

0.998

1

1

and no unit root, PHR had lower power. We know that the construction of PHR relies on
asymptotical properties and large enough sample sizes and frame lengths so a diminished
power is not surprising for small sample sizes. Since KPSS was designed to detect violations
due to unit root, perhaps its underperformance detecting varying variance is not surprising.
However, PHR had greater power to detect both the combination of unit root and varying
variance violations as well as only varying variance violation.
Simulation results for the cases with a single break in variance and with a smooth
transition in variance are consistent with the power results in Table 7 and we skipped the
tables for the sake of brevity but have them available upon request.
5.

Application. High Frequency Stock Market Returns.

When conducting statistical analyses of stock market rates of return, finance academics
typically assume that the return series are normally distributed and weakly stationary. If
24

a time series is stationary, the k time step ahead forecasts of the time series converge to
the time series mean, and the variance of the forecasting error converges to the time series
variance, as the forecasting horizon lengthens. A test that detects evidence against weak
stationarity also detects evidence against strict stationarity. On the other hand, if we fail to
reject weak stationarity, we can still reject strict stationarity, as strict stationarity requires
that higher order moments are also time invariant. Next we illustrate the test power to
detect nonstationarity in high frequency returns for one New York Stock Exchange (NYSE)
traded stock, the Ford Motor Company. We obtained our data from Wharton Research
Data Services.
We considered every Ford stock trade for the years 2001, 2002, and 2003. With trades
occurring at random times during market trading hours (9:30AM to 4:00PM), we sampled
prices at equal length intervals. We then transformed the prices by taking the first difference
of the natural logarithm in order to produce the stock rates of return. When conducting a
statistical analysis of stock market rates of return, finance academics typically assume that
the data are normally distributed and weakly stationary. In order to guard against outliers
unduly affecting the test we trimmed the lowest and largest one percent of the observations.
We had over 19,000 observations per year, after trimming. We also detrended the data in
order to remove any potential subtle trend in the trimmed returns. Table 8 shows the
Ford stock descriptive statistics, before any data processing. The means and standard
deviations are typical for high frequency return data. Positive kurtosis is also a common
characteristic. The series has small negative skewness in all years. Figure 3 shows one year
of return data with apparent time varying variance.
For each evaluated year of data, PHR rejected the null hypothesis of stationarity at five
percent significance level before and after demeaning, detrending and trimming. The KPSS
test of stationarity did not detect this pattern before trimming but did after trimming at
25

−0.2
−0.6

−0.4

5−minute returns

0.0

0.2

High Frequency Return Data for Ford Stock For 2002

0

5000

10000

15000

20000

time points

Figure 3: High Frequency Ford Returns 2002
one percent level.
When we adjusted the data using a rolling window estimate of variance like in Lima
and Neri (2013), the time varying variance violation was smoothed and PHR concurred
with KPSS in failing to reject stationarity. We estimated the smoothing correction factor
using a rolling window of 120 days.
Table 8: Descriptive Statistics for Ford Stock 5-Minute Returns
Year

N

Mean

St. Dev.

Min

Max

Skew.

2001
2002
2003

19,344

0.00000

0.010

−0.088

0.081

−0.136

6.47

19,656

−0.00001

0.031

−0.586

0.268

−0.664

13.46

19,655

0.00000

0.037

−0.369

0.239

−0.397

4.86

26

Kurt.

6.

Conclusion.

In this article we used the frame averaged estimates of the second order cumulant spectrum
calculated for pairs of fundamental Fourier frequencies in the principal domain to build a
statistical test of the null hypothesis that a stochastic process is strictly stationary. We
assume that the process has zero mean, is free of outliers, has no trend and has been
prewhitened. One example of a stationary process is white noise. We based the test
development on the property that the second order cumulant spectrum of a time series in
the frequency domain is zero under strict stationarity. An application of the test to five
minute sampled high frequency Ford stock returns showed nonstationarity.
We evaluated the size and power properties of PHR and KPSS via Monte Carlo simulations. We also evaluated the sensitivity of PHR to sample sizes, fat tails in the time
series distribution, and autocorrelation in the error term. Some limitations of PHR are
sensitivity to outliers, lower power to detect subtle unit root processes in smaller sample
sizes, and higher Type I error rate in certain data scenarios. PHR had good power to
detect second moment variations, and outperformed KPSS test in all considered scenarios. PHR also showed good power to detect a unit root in the presence of a time varying
second moment. In future research it would be interesting to study the power of PHR to
detect other types of stationarity violations, such as long memory violations as described
in Ashley and Patterson (2010). PHR can be included in a battery of tests used to verify
the strict stationarity assumptions of a time series.
References
Ashley, R. A. and Patterson, D. M. (2010). Apparent long memory in time series as an
artifact of a time-varying mean: Considering alternatives to the fractionally integrated
27

model. Macroeconomic Dynamics, 14(S1):59–87.
Bloomfield, P. (2004). Fourier Analysis of Time Series: an Introduction. John Wiley &
Sons.
Borchers, H. W. (2015). pracma: Practical Numerical Math Functions. R package version
1.8.8.
Brillinger, D. R. (2001). Time Series: Data Analysis and Theory, volume 36. Siam.
Cavaliere, G. and Taylor, A. R. (2005). Stationarity tests under time-varying second moments. Econometric Theory, 21(06):1112–1129.
Dickey, D. A. and Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root.

Journal of the American Statistical Association,

74(366a):427–431.
Dickinson, R. M., Roberts, D. A. O., Driscoll, A. R., Woodall, W. H., and Vining, G. G.
(2014). Cusum charts for monitoring the characteristic life of censored Weibull lifetimes.
Journal of Quality Technology, 46(4):340.
Farley, J. U. and Hinich, M. J. (1970). Detecting small mean shifts in time series. Management Science, 17(3):189–199.
Hawkins, D. M. (1977). Testing a sequence of observations for a shift in location. Journal
of the American Statistical Association, 72(357):180–186.
Hinich, M. A. and Wild, P. (2001). Testing time-series stationarity against an alternative
whose mean is periodic. Macroeconomic Dynamics, 5(03):380–412.

28

Hinich, M. J. (1994). Higher order cumulants and cumulant spectra. Circuits, Systems and
Signal Processing, 13(4):391–402.
Hinich, M. J. and Rothman, P. (1998). Frequency-domain test of time reversibility. Macroeconomic Dynamics, 2(01):72–88.
Hlavac, M. (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
Harvard University, Cambridge, USA. R package version 5.2.
Kwiatkowski, D., Phillips, P. C., Schmidt, P., and Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic
time series have a unit root? Journal of Econometrics, 54(1):159–178.
Lima, L. R. and Neri, B. (2013). A test for strict stationarity. Uncertainty Analysis in
Econometrics with Applications, pages 17–30.
Phillips, P. C. and Perron, P. (1988). Testing for a unit root in time series regression.
Biometrika, pages 335–346.
Press, W., Flannery, B., Teukolsky, S., and Vetterling, W. (1992). Kolmogorov-Smirnov
test. Numerical Recipes in FORTRAN: The Art of Scientific Computing, pages 617–620.
R Core Team (2015). R: A Language and Environment for Statistical Computing. R
Foundation for Statistical Computing, Vienna, Austria.
Trapletti, A. and Hornik, K. (2015). tseries: Time Series Analysis and Computational
Finance. R package version 0.10-34.
Wang, J., Tsang, W. W., and Marsaglia, G. (2003). Evaluating Kolmogorov’s distribution.
Journal of Statistical Software, 8(18).
29

