\section{ Introduction.}
\label{sec:intro}

A stochastic process is called strictly stationary if the joint distribution of $x\left(t_1 + \tau\right)$, $x\left(t_2 + \tau\right)$, $\ldots$, $x\left(t_n + \tau\right)$
is independent of $\tau$ for all time points $\{t_1,$ $t_2,$ $\ldots,$ $t_n\}$. 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.  \cite{kwiatkowski1992testing} 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  \cite{brillinger2001time}).
Hence a process is weakly stationary (or stationary in the wide sense), if its
expected value $E\{x\left(t_m\right)\} $ is constant and its autocovariance $c_x\left(t_m\right)$ depends only on the difference $\left(t_{n+m} - t_{n}\right)$, $c_x\left(t_m\right) = c_x\left(t_{n+m} - t_{n}\right)$ (\cite{brillinger2001time}). A weakly stationary Gaussian process is also stationary in the strict sense \cite{brillinger2001time}. A Gaussian process is a sequence of independent, identically normally distributed random variables $\{u_t\}$.

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\left(t\right)\}$ 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 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\left(t\right)\}$ 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.

\subsection{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, \cite{farley1970detecting} 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. \cite{hawkins1977testing} proposed a likelihood ratio based procedure to detect a change in the mean while assuming constant variance. In a different approach, \cite{dickinson2014cusum} 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 
\cite{dickey1979distribution} and \cite{phillips1988testing}.

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 \cite{hinich2001testing}. In the space of tests for strict stationarity in the time domain, \cite{lima2013test} proposed a test for strict stationarity and compared its power to the KPSS test of level stationarity. Similarly to \cite{lima2013test} 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 \ref{sec:cum2}, the test development in Section \ref{sec:derivation}, empirical size and power evaluation via Monte Carlo simulations and comparison with the KPSS test under various alternative hypotheses in Section \ref{sec:evaluation}. The test is demonstrated in \ref{sec:sec5} using high frequency (trade by trade) stock returns for the Ford company . 


\section{The Second Order Cumulant Spectrum of a Strictly Stationary Process.}
\label{sec:cum2}

In the frequency domain, the second order cumulant spectrum of a stochastic process $\{x(t_n)\}$  with $t_n = n\tau$
\begin{equation}
K\left(f_1, f_2\right) = \sum_{n_1 = -\infty}^{\infty}\sum_{n_1 = -\infty}^{\infty}E\left[x\left(t_{n_1}\right)x\left(t_{n_2}\right)\right]exp\left[-2\pi\left(f_1t_{n_1} + f_2t_{n_2}\right)\right],
\end{equation}

\noindent where $f_i$ are discrete Fourier frequencies. 

The cumulant spectrum $K$ is defined for frequency values in the square $\{- f_0 \le f_1 \le f_0, - f_0 \le f_2 \le f_0 \}$,  where $f_0 = 1 / \left(2 \tau\right)$ is called the Nyquist frequency. Fourier frequencies display the following properties and symmetries:
\begin{itemize}
\item Frequencies are between zero and one, normalized by Nyquist frequency. The Nyquist frequency $f_0$ 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 $\tau$ is one. When $\tau = 1$ the Nyquist frequency  is $f_0 = 0.5$. 
\item The discrete Fourier transform of a sample of length T at frequency $f$, $X\left(f\right) = \sum_{t=0}^{T - 1}x\left(t \right)\exp\{-i2\pi f t \}$ is equal to the conjugate of the Fourier transform at $-f$, $X(f) = X^{*}(-f)$.
\item The Fourier transform is periodic at intervals of the sampling frequency for discrete time processes. Higher frequencies fold down into the $\left(0, f_0\right)$ interval.
\end{itemize}

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 (\cite{bloomfield2004fourier}), also called fundamental frequencies. The principal domain forms the triangle $\{0 < f_1 \le f_0, -f_1 < f_2 \le f_1\}$.

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 \cite{hinich1994higher}, if the time series $\{x\left(t\right)\}$ is strictly stationary then the second order cumulant spectrum for a pair of frequencies in the principal domain $\left(f_1, f_2\right)$, $K\left(f_1, f_2 \right) $is equal to

\begin{equation} \label{eq:4}
K\left(f_1, f_2 \right)= E[X(f_1)X(f_2)] = 
\left\{
	\begin{array}{ll}
		2S\left( f _1\right) + O\left( 1\right) & \mbox{if } f_1+f_2 = \mbox{ 0 mod  1 }\\			
		O\left(1\right) & \mbox{ otherwise, }
	\end{array}
\right.
\end{equation}

\noindent where $X(f_i)$ denote the discrete Fourier transform of the time series at the frequencies included in the principal domain. The line $f_2 = - f_1$ is
not in the principal domain (\cite{hinich1994higher} ). From Equation \ref{eq:4}, $K\left(f_1, -f_1\right) = \delta\left(f_1\right)S\left(f_1\right)$,
where $\delta\left(f_1\right)$ is the Dirac delta function, and $S\left(f_1\right) = \sum_{m = -\infty}^{\infty}c_x\left(t_m\right)\exp\left(-i2\pi f_1t_m\right)$ is the spectrum at frequency $f_1$ and  $c_x\left(t_m\right)$ is the autocovariance of the process. 

In summary, if a time series is strictly stationary, then its second order cumulant spectrum

\begin{equation} \label{eq:2}
K\left(f_1, f_2\right) = 0
\end{equation}

\noindent for all frequency pairs $\left(f_1, f_2 \right)$ in the principal domain \cite{hinich1994higher}. 



\section{ Derivation of The Test Statistic for the Second Order Cumulant Spectrum Test of Strict Stationarity.}
\label{sec:derivation}

We now use the second order cumulant spectrum to construct the test for the null hypothesis that a stochastic process $\{x(t_n)\}$ 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(f_1, f_2)$ and for the spectrum $S(f)$, where $(f_1, f_2)$ 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

\begin{equation} \label{eq:9}
\hat{K}\left(f_1 , f_2 \right) = \frac{1}{P} \sum_{p=1}^{P} K_p\left(f_1 , f_2 \right),
\end{equation}

We can also express the principal domain for $K$ in terms of integers $\left(k_1, k_2\right)$, with $f = k/L$, and $0 < k_1 \le L/2$, $-k_1 \le k_2 \le k_1$.  The steps to get the frame averaged estimates of the cumulant spectrum, $\hat{K}$, and of the spectrum, $\hat{S}$, are calculated following the steps: 

\begin{itemize}
\item Use a reasonable length time series, ideally of no less than $T = 150$ observations.
\item 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 $\{x_p(1), x_p(2), ..., x_p(L)\}$.
\item Calculate the discrete Fourier transform (DFT) for the $p$ frame of length $L$ as 

\begin{equation} \label{eq:6}
X_p\left(k\right) = \sum_{t=0}^{L - 1}x\left[t + \left(p L \right)\right]exp\{-i2\pi f_k \left[t + \left(p  L\right)\right]\},
\end{equation}

where $f_k = \frac{k}{L}$.

\item Compute the spectrum at each frequency in the principal domain of $K$ for each frame $p$

\begin{equation} \label{eq:7}
S_p\left(k\right) = \frac{1}{L} |X_p \left(k\right)|^2 = \frac{1}{L} X_p\left(k\right)X_p\left(-k\right)
\end{equation}

\item Compute the second order cumulant spectrum for each frame $p$

\begin{equation} \label{eq:8}
K_p\left(k_1 , k_2 \right) = \frac{1}{L} X_p \left(k_1\right)X_p \left(k_2\right)
\end{equation}

for pairs of frequencies defined by $0 < k_1 \le L/2$ and $-k_1 < k_2 \le k_1$.

\item Get the frame averaged $\hat S$ and $\hat K$ over the P frames for each $k$ and $(k_1, k_2)$ pair respectively.

\end{itemize}

From the theory in \cite{hinich1994higher} and in \cite{brillinger2001time}, 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})$ ,

\begin{equation} \label{eq:10}
E\left[\hat{S}\left(f_k\right)\right] = S\left(f_k\right) + O\left( L^{-1}\right).
\end{equation}

Similarly to Equation \ref{eq:10}, under strict stationarity, from \cite{hinich1994higher}

\begin{equation} \label{eq:11}
E\left[\hat{K}\left(k_1 , k_2 \right)\right] = K\left(k_1 , k_2 \right) + O\left(L^{-1}\right).
\end{equation}

We stated in Equation ~\ref{eq:2} that the second order cumulant spectrum $K\left(k_1 , k_2 \right)$ for strictly stationary time series is zero for pairs of frequencies in the frequency domain. Furthermore, the variance of $\hat{K}$, considering $P$ frames, is equal to $P^{-1}S\left({k_1}\right)S\left({k_2}\right)$ as $\left(L, P \to \infty\right)$ from \cite{hinich1994higher}.

Next, we define the normalized second order cumulant spectrum as 
\begin{equation} \label{eq:12}
\hat{\Gamma} \left(k_1, k_2 \right) = \frac{\hat{K}\left(k_1 , k_2 \right)}{\sqrt{S\left(f_{k_1}\right)S\left(f_{k_2}\right)}}.
\end{equation}

We make a further simplifying assumption used in \cite{hinich1998frequency} 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 $\hat{K}$ becomes $1/\sqrt P$. Then the estimate for $\Gamma$, $\hat{\Gamma}$ simplifies to $\sqrt{P}\hat{K}$. 
 
Next we construct the complex valued pivotal quantity $Y$ as

\begin{equation} \label{eq:13}
\hat Y\left(k_1, k_2\right) = \sqrt {2P} \left[ \hat{\Gamma}\left(k_1, k_2 \right) - \Gamma \left(k_1, k_2 \right) \right].
\end{equation}

From the central limit theorem, the pivotal quantity $Y$ has an asymptotic complex normal distribution. Hence, the real value $\Re{(Y)}$ and the imaginary value $\Im{(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 $\left(k_1, k_2\right)$ in the principal domain. The number of calculated cumulants spectra equals the number of pairs $\left(k_1, k_2\right)$ in the principal 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 $F_{emp}\left(x\right)$ denote the empirical cumulative distribution of the vector of stacked $\left(\Re{(Y)}, \Im{(Y)} \right)$. We calculate  $F_{emp}\left(x\right)$  as the fraction of values that are less than a given
value $\emph{x}$ in the unit interval. In our case the vector $\mathbf 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 $\hat{D} = sup_{0 < x \le 1}|{F_{emp}\left(\mathbf x\right) - \mathbf x}|$  (\cite{press1992kolmogorov}).


\section{Empirical Evaluation of Size and Power of the PHR Test. Comparison to the KPSS Test.}
\label{sec:evaluation}

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 \cite{tseries15} tseries R package. Developed in \cite{kwiatkowski1992testing}, the KPSS test tests the null hypothesis of first order (level or absence of unit root) stationarity using the test statistic

\begin{equation} \label{eq:15}
KPSS = \frac{1}{\hat{\sigma}^2 T^2}\sum_{j=1}^{T}\left(\sum_{t=1}^{j}\left(x(t) - \bar {x}(t)\right)\right)^2,
\end{equation}

\noindent where $x_t$ is the time series, with $t = 1 \ldots T$,  where $T$ is the length of the times series, $\bar{x}(t)$ is the mean of the time series, and $\hat{\sigma}$ is a consistent estimator of the long run variance. 

 \cite{lima2013test} studied the size and power of KPSS to detect various departures from the strict stationarity hypothesis. \cite{cavaliere2005stationarity} evaluated KPSS performance to detect unit root in the presence of varying second moment. In the following sections we replicate several of the data generating processes in \cite{lima2013test} and \cite{cavaliere2005stationarity} 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 \cite{R15} 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 $\hat{Y}$ test statistic in Equation \ref{eq:13} 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 \cite{wang2003evaluating}. 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. We performed all the simulations in the R open source software and we created all the tables using the stargazer \cite{stargazer15} 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 \cite{pracma15} 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. 
 
\subsection{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 \ref{tab:tabthree} 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$. 
 
\begin{table}[!htbp] \centering 
  \caption{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.} 
  \label{tab:tabthree} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{L} & \multicolumn{1}{c}{T.250} & \multicolumn{1}{c}{T.500} & \multicolumn{1}{c}{T.1000} & \multicolumn{1}{c}{T.5000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{10} & 0.062 & 0.059 & 0.062 & 0.058 \\ 
\multicolumn{1}{c}{20} & 0.062 & 0.058 & 0.060 & 0.058 \\ 
\multicolumn{1}{c}{50} & 0.195 & 0.084 & 0.063 & 0.057 \\ 
\multicolumn{1}{c}{100} & 1 & 0.716 & 0.183 & 0.059 \\ 
\multicolumn{1}{c}{200} & 1 & 1 & 1 & 0.122 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

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 $\infty$ (normal), three, five, ten and 15 degrees of freedom. Furthermore, we included an autocorrelation structure 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.

\subsubsection{Data Generation Processes for Size of Test Simulations.}

The stationary data generation process is

\begin{equation} \label{eq:16}
x(t) = u(t).
\end{equation}

\noindent The innovations $u(t)$ follow the process

\begin{equation} \label{eq:17}
u(t)= \rho u({t-1}) + v(t),
\end{equation}

\noindent where $v(t)$ are independently and identically distributed with one of the  $t$ distributions. When $\rho$ is zero, $u(t)$ in Equation \ref{eq:17} process reduces to just the $v(t)$ component. If $\rho$ is different from zero the time series has short term memory. 

\subsubsection{Simulation Results for Empirical Test Size.} 
	
Table \ref{tab:tabfour} 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. 
	
\begin{table}[!htbp] \centering 
  \caption{Empirical Size of Test Comparisons for Normal Data with Zero Autocorrelation in Error Term($\rho = 0$) and Nominal Size $\alpha = 0.05$.} 
  \label{tab:tabfour} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{itrim} & \multicolumn{1}{c}{idetrend} & \multicolumn{1}{c}{idemean} & \multicolumn{1}{c}{KPSS.250} & \multicolumn{1}{c}{PHR.250} & \multicolumn{1}{c}{KPSS.500} & \multicolumn{1}{c}{PHR.500} & \multicolumn{1}{c}{KPSS.1000} & \multicolumn{1}{c}{PHR.1000} & \multicolumn{1}{c}{KPSS.5000} & \multicolumn{1}{c}{PHR.5000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.049 & 0.062 & 0.052 & 0.084 & 0.050 & 0.063 & 0.052 & 0.059 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.050 & 0.051 & 0.051 & 0.070 & 0.051 & 0.044 & 0.051 & 0.043 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0 & 0.062 & 0 & 0.085 & 0 & 0.064 & 0 & 0.059 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0 & 0.049 & 0 & 0.069 & 0 & 0.046 & 0 & 0.042 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.049 & 0.062 & 0.052 & 0.084 & 0.050 & 0.063 & 0.052 & 0.059 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.050 & 0.051 & 0.051 & 0.070 & 0.051 & 0.044 & 0.051 & 0.043 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0 & 0.062 & 0 & 0.085 & 0 & 0.064 & 0 & 0.059 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0 & 0.049 & 0 & 0.069 & 0 & 0.046 & 0 & 0.042 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

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 \ref{tab:tabfive} 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. 

\begin{table}[!htbp] \centering 
  \caption{Empirical Size of Test Comparisons for t5 Data with Zero Autocorrelation in Error Term($\rho = 0$) and Nominal Size $\alpha = 0.05$.} 
  \label{tab:tabfive} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{itrim} & \multicolumn{1}{c}{idetrend} & \multicolumn{1}{c}{idemean} & \multicolumn{1}{c}{KPSS.250} & \multicolumn{1}{c}{PHR.250} & \multicolumn{1}{c}{KPSS.500} & \multicolumn{1}{c}{PHR.500} & \multicolumn{1}{c}{KPSS.1000} & \multicolumn{1}{c}{PHR.1000} & \multicolumn{1}{c}{KPSS.5000} & \multicolumn{1}{c}{PHR.5000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.046 & 0.152 & 0.048 & 0.266 & 0.048 & 0.247 & 0.048 & 0.280 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.048 & 0.066 & 0.048 & 0.106 & 0.047 & 0.068 & 0.048 & 0.064 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0 & 0.152 & 0 & 0.268 & 0 & 0.246 & 0 & 0.280 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0.001 & 0.067 & 0.001 & 0.105 & 0.001 & 0.066 & 0.001 & 0.064 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.046 & 0.152 & 0.048 & 0.266 & 0.048 & 0.247 & 0.048 & 0.280 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.048 & 0.066 & 0.048 & 0.106 & 0.047 & 0.068 & 0.048 & 0.064 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0 & 0.152 & 0 & 0.268 & 0 & 0.246 & 0 & 0.280 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0.001 & 0.067 & 0.001 & 0.105 & 0.001 & 0.066 & 0.001 & 0.064 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

Next, we evaluated PHR and KPSS empirical size sensitivity to a $\rho = 0.3$ autocorrelation 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 \ref{tab:tabsix}  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.

\begin{table}[!htbp] \centering 
  \caption{Empirical Size of Test Comparisons for Normal Data with Autocorrelation in Error Term $\rho = 0.3$ and Nominal Size of Test $\alpha = 0.05$.} 
  \label{tab:tabsix} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{itrim} & \multicolumn{1}{c}{idetrend} & \multicolumn{1}{c}{idemean} & \multicolumn{1}{c}{KPSS.250} & \multicolumn{1}{c}{PHR.250} & \multicolumn{1}{c}{KPSS.500} & \multicolumn{1}{c}{PHR.500} & \multicolumn{1}{c}{KPSS.1000} & \multicolumn{1}{c}{PHR.1000} & \multicolumn{1}{c}{KPSS.5000} & \multicolumn{1}{c}{PHR.5000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.048 & 0.062 & 0.050 & 0.084 & 0.047 & 0.061 & 0.047 & 0.058 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.045 & 0.055 & 0.051 & 0.067 & 0.048 & 0.044 & 0.049 & 0.040 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0 & 0.061 & 0 & 0.084 & 0 & 0.061 & 0 & 0.058 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0 & 0.055 & 0 & 0.068 & 0 & 0.045 & 0 & 0.040 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.048 & 0.062 & 0.050 & 0.084 & 0.047 & 0.061 & 0.047 & 0.058 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.045 & 0.055 & 0.051 & 0.067 & 0.048 & 0.044 & 0.049 & 0.040 \\  
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0 & 0.061 & 0 & 0.084 & 0 & 0.061 & 0 & 0.058 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0 & 0.055 & 0 & 0.068 & 0 & 0.045 & 0 & 0.040 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

When the data generation process included negative $-0.3$ autocorrelation in the error term, the empirical size simulation results were similar to Table \ref{tab:tabsix} results.


\subsection{Power of Test to Detect Unit Root.}
\label{sec:ur}

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

\subsubsection{Data Generation Processes With Unit Root.}

 We replicated the unit root data generation process from \cite{lima2013test}. The generated data included a unit root component $y(t)$ and a stationary component $u(t)$. The tuning parameter $\lambda$ moderates the strength of the unit root process. Then the time series $x(t)$ is 

\begin{equation} \label{eq:23}
x(t) = \lambda  y(t) + u(t),
\end{equation}

 \noindent with $ y(t) = \sum_{j=1}^{t}v(j)$. 

The innovations $u(t)$ follow the process in Equation \ref{eq:17}. 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 $\lambda$ of $0.1$ and $0.3$. A larger $\lambda$ led to $100\%$ power for both tests, as observed in additional simulation studies not included here.


\subsubsection{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 $\lambda$ value), having difficulty with smaller $\lambda$ values. We also evaluated the impact of trimming, demeaning and detrending on the power of the test.

\begin{table}[!htbp] \centering 
  \caption{Power of Test Comparisons for Normal Data with a Unit Root ($\lambda = 0.3$), Zero Autocorrelation in Error Term($\rho = 0$), and Nominal Size of Test $\alpha = 0.05$.} 
  \label{tab:tabnine} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{itrim} & \multicolumn{1}{c}{idetrend} & \multicolumn{1}{c}{idemean} & \multicolumn{1}{c}{KPSS.5000} & \multicolumn{1}{c}{PHR.5000} & \multicolumn{1}{c}{KPSS.250} & \multicolumn{1}{c}{PHR.250} & \multicolumn{1}{c}{KPSS.1000} & \multicolumn{1}{c}{PHR.1000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 1 & 0.997 & 0.963 & 0.080 & 0.997 & 0.420 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 1 & 0.996 & 0.956 & 0.072 & 0.996 & 0.406 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0.998 & 0.997 & 0.620 & 0.078 & 0.936 & 0.393 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0.998 & 0.997 & 0.597 & 0.072 & 0.928 & 0.376 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 1 & 0.997 & 0.963 & 0.080 & 0.997 & 0.420 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 1 & 0.996 & 0.956 & 0.072 & 0.996 & 0.406 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0.998 & 0.997 & 0.620 & 0.078 & 0.936 & 0.393 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0.998 & 0.997 & 0.597 & 0.072 & 0.928 & 0.376 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

Table ~\ref{tab:tabnine} illustrates the good performance of PHR in large sample sizes for normal data. 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. 

\begin{table}[!htbp] \centering 
  \caption{Power of Test Comparisons for t5 Data with Unit Root ($\lambda = 0.3$), Zero Autocorrelation in Error Term($\rho = 0$), and Nominal Size of Test $\alpha = 0.05$.} 
  \label{tab:tabten} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{itrim} & \multicolumn{1}{c}{idetrend} & \multicolumn{1}{c}{idemean} & \multicolumn{1}{c}{KPSS.250} & \multicolumn{1}{c}{PHR.250} & \multicolumn{1}{c}{KPSS.1000} & \multicolumn{1}{c}{PHR.1000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.968 & 0.148 & 0.998 & 0.573 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & 0.965 & 0.106 & 0.997 & 0.534 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0.628 & 0.148 & 0.939 & 0.550 \\  
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & 0.623 & 0.098 & 0.935 & 0.493 \\  
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.968 & 0.148 & 0.998 & 0.573 \\  
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & 0.965 & 0.106 & 0.997 & 0.534 \\ 
\multicolumn{1}{c}{0} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0.628 & 0.148 & 0.939 & 0.550 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & \multicolumn{1}{c}{1} & 0.623 & 0.098 & 0.935 & 0.493 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

In Table ~\ref{tab:tabten}, we generated data from a t distribution with five degrees of freedom, no correlation in error term and a unit root with $\lambda = 0.3$. PHR performed similarly to the normal case, displaying the already established sensitivity to sample size and fat tails. In Table \ref{tab:tabten} 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 $\lambda$ decreased to $0.1$ (weaker presence of the unit root process). Both PHR and KPSS performed well for larger sample sizes. 

\subsection{Power of the Test to Detect Time Varying Second Moment.}
\label{sec:vol}

Another alternative to the null hypothesis of stationarity is a varying second moment.  The authors of \cite{cavaliere2005stationarity} 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 \cite{cavaliere2005stationarity}, and similar to \cite{lima2013test}. 
 
 \cite{lima2013test} also evaluated several tests for strict stationarity and their power against the alternative of unconditional heteroskedasticity.  \cite{lima2013test} 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 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. \cite{lima2013test} only discussed one pattern of heteroskedasticity or varying second moment while we considered the three patterns suggested in \cite{cavaliere2005stationarity}. 


\subsubsection{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 \ref{eq:21} is for a level stationary (no unit root) process with variance $\sigma_t$. Under the null hypothesis of time constant variance, $\sigma_t ^ 2= 1$. Under the alternative hypothesis of varying variance, $\sigma_t^2$ follows one of the three processes in Equation \ref{eq:19}.  DGP1 is

 \begin{equation} \label{eq:21}
 x(t) =  \sigma_t u(t), 
 \end{equation}
 
\noindent with $t = 1 \ldots T$ and $u(t)$ a white noise.
  
 The second data generating process (DGP2) in Equation \ref{eq:22}  generates a process with unit root and variance $\sigma_t ^2$. Similarly to DGP1, $\sigma_t ^ 2= 1$ under the null hypothesis of constant variance. Under the alternative hypothesis of varying variance, $\sigma_t ^ 2$ follows one of the three processes in Equation \ref{eq:19}.  DGP2 is
  
 \begin{equation} \label{eq:22}
 x(t) = x({t-1}) + \sigma_t u(t),
 \end{equation}
 
\noindent 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 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 $\sigma_t ^ 2$ to vary according to the three patterns described below, while the initial paper by \cite{cavaliere2005stationarity} 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.
 


 
\begin{itemize}
\item Single Break 
\begin{equation*}  \label{eq:18}
 \sigma^2_t = \begin{cases} {\sigma_0}^2 + \left( {\sigma_1}^2 - {\sigma_0}^2\right),  & t \ge m \\
 {\sigma_0}^2, & t < m
\end{cases}
\end{equation*}

where $\sigma_0 = 1$ and $\sigma_1 = c\sigma_0$.
 
 \item Smooth Transition
 \begin{equation} \label{eq:19}
  \sigma^2_t = {\sigma_0}^2 + \left( {\sigma_1}^2 - {\sigma_0}^2\right)W\left(t\right),
 \end{equation}
 
 where $W\left(t\right) = \left(1+\exp\left(-\gamma\left(t-m\right)\right)\right)^{-1}$,  $\sigma_0 = 1$ and $\sigma_1 = c\sigma_0$.
 \item
  Piecewise Linear Trend
  \begin{equation*}  \label{eq:20}
 \sigma^2_t = \begin{cases} {\sigma_0}^2 + \left( {\sigma_1}^2 - {\sigma_0}^2\right)\left(t-m\right)  \left(1-m\right)^{-1},  & t \ge m \\
 {\sigma_0}^2, & t < m
\end{cases}
\end{equation*}

where $\sigma_0 = 1$ and $\sigma_1 = c\sigma_0$.

\end{itemize}
 
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 $\gamma = 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 $\sigma$ 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$. 

 \begin{figure} [!htbp] \label{fig:fig1}
\begin{center}
\includegraphics[height=3in, width=3in, keepaspectratio]{Rplot02.pdf}
\end{center}
\caption{The Three Functions for  a Variance Increase and T = 1000. Left: Single Break and Smooth Transition. Right: Piecewise Linear Trend.}
\end{figure}


Figure 2 illustrates both DGP1 and DGP2 processes with time varying variance according to the three patterns. 
 
 \begin{figure} [!htbp] \label{fig:fig2}
\begin{center}
\includegraphics[height=5in, width=5in, keepaspectratio]{VaryingVariance.pdf}
\end{center}
\caption{Time Series with Time Varying Variance}
\end{figure}


\subsubsection{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 \ref{tab:tabseven} 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 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.

\begin{table}[!htbp] \centering 
  \caption{Power of Tests to Detect Piecewise Linear Trend at Nominal Test Size $\alpha = 0.05$.} 
  \label{tab:tabseven} 
\scriptsize
\begin{tabular}{@{\extracolsep{0pt}} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} } 
\\[-1.8ex]\hline \\[-1.8ex] 
\multicolumn{1}{c}{dgp1\_2} & \multicolumn{1}{c}{m} & \multicolumn{1}{c}{c} & \multicolumn{1}{c}{PHR.250} & \multicolumn{1}{c}{KPSS.250} & \multicolumn{1}{c}{PHR.1000} & \multicolumn{1}{c}{KPSS.1000} & \multicolumn{1}{c}{PHR.5000} & \multicolumn{1}{c}{KPSS.5000} \\ 
\hline \\[-1.8ex] 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{4} & 0.090 & 0.054 & 0.742 & 0.056 & 1 & 0.055 \\ 
\multicolumn{1}{c}{2} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{4} & 0.995 & 0.982 & 1 & 0.998 & 1 & 1 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0.25} & 0.094 & 0.050 & 0.762 & 0.051 & 1 & 0.054 \\ 
\multicolumn{1}{c}{2} & \multicolumn{1}{c}{0} & \multicolumn{1}{c}{0.25} & 0.996 & 0.982 & 1 & 0.997 & 1 & 1 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0.5} & \multicolumn{1}{c}{4} & 0.273 & 0.098 & 1 & 0.099 & 1 & 0.100 \\ 
\multicolumn{1}{c}{2} & \multicolumn{1}{c}{0.5} & \multicolumn{1}{c}{4} & 0.987 & 0.966 & 1 & 0.995 & 1 & 1 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0.5} & \multicolumn{1}{c}{0.25} & 0.078 & 0.036 & 0.378 & 0.036 & 0.956 & 0.039 \\ 
\multicolumn{1}{c}{2} & \multicolumn{1}{c}{0.5} & \multicolumn{1}{c}{0.25} & 0.997 & 0.984 & 1 & 0.998 & 1 & 1 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0.9} & \multicolumn{1}{c}{4} & 0.447 & 0.114 & 1 & 0.120 & 1 & 0.125 \\ 
\multicolumn{1}{c}{2} & \multicolumn{1}{c}{0.9} & \multicolumn{1}{c}{4} & 0.990 & 0.974 & 1 & 0.996 & 1 & 1 \\ 
\multicolumn{1}{c}{1} & \multicolumn{1}{c}{0.9} & \multicolumn{1}{c}{0.25} & 0.066 & 0.042 & 0.084 & 0.042 & 0.105 & 0.047 \\ 
\multicolumn{1}{c}{2} & \multicolumn{1}{c}{0.9} & \multicolumn{1}{c}{0.25} & 0.998 & 0.984 & 1 & 0.998 & 1 & 1 \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 

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 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 \ref{tab:tabseven} and we skipped the tables for the sake of brevity but have them available upon request.

\section{Application. High Frequency Stock Market Returns.}
\label{sec:sec5}

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 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 \ref{tab:tabtwo} 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.

 
 \begin{figure} [!htbp]
 \label{fig:figone}
\begin{center}
\includegraphics[height=3in, width=3in, scale = 0.5]{Ford02.pdf}
\end{center}
\caption{High Frequency Ford Returns 2002 \label{fig:figone}}
\end{figure}


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 one percent level.  

When we adjusted the data using a rolling window estimate of variance like in \cite{lima2013test}, 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.

\begin{table}[!htbp] \centering 
  \caption{Descriptive Statistics for Ford Stock 5-Minute Returns} 
  \label{tab:tabtwo} 
  \scriptsize
\begin{tabular}{@{\extracolsep{5pt}}lcccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
Year & \multicolumn{1}{c}{N} & \multicolumn{1}{c}{Mean} & \multicolumn{1}{c}{St. Dev.} & \multicolumn{1}{c}{Min} & \multicolumn{1}{c}{Max} & \multicolumn{1}{c}{Skew.} & 
\multicolumn{1}{c}{Kurt.}\\ 
\hline \\[-1.8ex] 
2001 & 19,344 & 0.00000 & 0.010 & $-$0.088 & 0.081 & $-$0.136 & 6.47\\ 
2002 & 19,656 & $-$0.00001 & 0.031 & $-$0.586 & 0.268 & $-$0.664 & 13.46\\ 
2003 & 19,655 & 0.00000 & 0.037 & $-$0.369 & 0.239 & $-$0.397 & 4.86\\
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 





\section{Conclusion.}
\label{sec:conc}

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 \cite{ashley2010apparent}. PHR can be included in a battery of tests used to verify the strict stationarity assumptions of a time series.



\bibliographystyle{apalike}

