# Non-ergodic statistics and spectral density estimation for stationary real
harmonizable symmetric $\alpha$-stable processes

Ly Viet Hoang<EMAIL_ADDRESS>Evgeny Spodarev<EMAIL_ADDRESS>Ulm University

###### Abstract

We consider stationary real harmonizable symmetric $\alpha$-stable processes
$X=\left\\{X(t):t\in\mathbb{R}\right\\}$ with a finite control measure.
Assuming the control measure is symmetric and absolutely continuous with
respect to the Lebesgue measure on the real line, we refer to its density
function as the spectral density of $X$. Standard methods for statistical
inference on stable processes cannot be applied as harmonizable stable
processes are non-ergodic.

A stationary real harmonizable symmetric $\alpha$-stable process $X$ admits a
LePage series representation and is conditionally Gaussian which allows us to
derive the non-ergodic limit of sample functions on $X$. In particular, we
give an explicit expression for the non-ergodic limits of the empirical
characteristic function of $X$ and the lag process
$\left\\{X(t+h)-X(t):t\in\mathbb{R}\right\\}$ with $h>0$, respectively.

The process admits an equivalent representation as a series of sinusoidal
waves with random frequencies whose probability density function is in fact
the (normalized) spectral density of $X$. Based on strongly consistent
frequency estimation using the periodogram we present a strongly consistent
estimator of the spectral density. The periodogram’s computation is fast and
efficient, and our method is not affected by the non-ergodicity of $X$. Most
notably no prior knowledge on parameters of the process such as its index of
stability $\alpha$ is needed.

###### keywords:

non-ergodic statistics , stationary process , stable process , harmonizable
process , non-ergodic process , Fourier analysis , frequency estimation ,
spectral density estimation , periodogram

††journal: Bernoulli

##  1 Introduction

Stationarity is often a key property in the analysis of dependence structures
of time series or more generally stochastic processes. Especially stationary
Gaussian processes have been extensively studied [4, 5, 10]. It is well known
that $\alpha$-stable processes, where $\alpha\in(0,2]$ is the so-called index
of stability, generalize Gaussian processes [17]. In particular, symmetric
$\alpha$-stable ($S\alpha S$) processes are of interest. It is the infinite
second moment of $\alpha$-stable distributions (in the case $\alpha<2$) which
sets them apart from Gaussian processes ($\alpha=2$), and allows them to be
used for example in models of heavy-tailed random phenomena.

It can be shown that stationary $S\alpha S$ processes fall into one of the
following three subclasses – $S\alpha S$ moving average processes,
harmonizable $S\alpha S$ processes and a third class, which consists of
processes characterized by a non-singular conservative flow and the
corresponding cocycle [14]. We will omit further details on this third class
here, as our focus lies on the class of harmonizable $S\alpha S$ processes,
and refer to [14, 16] for the in-depth structural analysis of stationary
$S\alpha S$ processes as well as the definitions of ergodic flows and
cocycles, respectively.

The classes of moving averages and harmonizable processes are disjoint when
the index of stability $\alpha$ is less then $2$. Only in the case $\alpha=2$,
i.e. in the Gaussian case, one may find both a moving average representation
and a harmonizable representation for the same process [17, Theorem 6.7.2].
Between moving average and the harmonizable processes, the latter have
attracted little attention from a statistical point of view so far. This is
mainly due to their non-ergodicity [3].

In this paper, we consider the class of stationary real harmonizable symmetric
$\alpha$-stable (SRH $S\alpha S$) processes $X=\\{X(t):t\in\mathbb{R}\\}$ with
index of stability $0<\alpha<2$, which are defined by

$\displaystyle X(t)=Re\left(\int_{\mathbb{R}}e^{itx}M(dx)\right),$

where $M$ is an isotropic complex random measure with circular control measure
$k$. This measure is a product measure on the space $(\mathbb{R}\times
S^{1},\mathcal{B}(\mathbb{R})\times\mathcal{B}(S^{1}))$ and admits the form
$k=m\cdot\gamma$. The measure $m$ is called the _control measure_ of $M$ and
$\gamma$ is the uniform probability measure on $S^{1}$. The finite dimensional
distributions of the process $X$ are determined by $m$, hence the process $X$
is completely characterized by $m$ as well. We assume that the control measure
$m$ is an absolutely continuous symmetric probability measure on $\mathbb{R}$
with symmetric probability density function $f$, which we refer to as the
_spectral density_ of the process $X$.

In this paper, an approach for the estimation of the spectral density based on
classical methods from spectral analysis and signal processing [2, 13] is
presented. Furthermore, we examine the asymptotic behaviour of time-averages
of observables on SRH $S\alpha S$ processes and derive an ergodic theorem.
This allows us to give an explicit expression for the non-ergodic limit of the
finite-dimensional empirical characteristic function of the process $X$.

First, Section 2 establishes the basics on SRH $S\alpha S$ processes. An in-
depth definition of the process as well as its properties are given. In
particular, the non-ergodicity of SRH $S\alpha S$ processes is discussed. We
also give a short introductory example involving $\alpha$-sine transforms on
$\mathbb{R}$, and how their inversion can be used to estimate the spectral
density.

In Section 3, the non-ergodicity of SRH $S\alpha S$ processes is examined in
more detail. From a LePage-type series representation it follows that SRH
$S\alpha S$ processes are conditionally Gaussian, albeit still non-ergodic.
This underlying Gaussian structure can be used to study the non-ergodic
behaviour of time-averages of observables of $X$. In particular, explicit
expressions for the non-ergodic limit of the finite-dimensional empirical
characteristic functions of the process $X$ and the lag process
$\\{X(t+h)-X(t):t\in\mathbb{R}\\}$ are derived.

The fourth section dives deeper in the underlying Gaussian structure of SRH
$S\alpha S$ processes. As a consequence of the LePage type series
representation, SRH $S\alpha S$ processes can be generated by a series of
sinusoidal waves with random amplitudes, phases and frequencies. The
frequencies are independently and identically distributed with the spectral
density $f$ as their probability density function. The periodogram is a
standard tool for the estimation of frequencies in signal theory. The
locations of the peaks in the periodogram are strongly consistent estimators
of the aforementioned frequencies, and in conjunction with kernel density
estimators they can be used to estimate the spectral density $f$ of the SRH
$S\alpha S$ process $X$. Under minimal assumptions on the kernel function the
strong consistency of the frequency estimators translates to strong or weak
consistency of the kernel density estimator for $f$.

In the final section of this paper we present a thorough numerical analysis
and a collection of numerical examples. Furthermore, some minor numerical
challenges of the spectral density estimation are discussed.

##  2 Preliminaries

Consider the probability space $(\Omega,\mathcal{F},P)$, and denote by
$L^{0}(\Omega)$ the space of real-valued random variables on this probability
space. Furthermore, define the space of complex-valued random variables on
$(\Omega,\mathcal{F},P)$ by
$L^{0}_{c}(\Omega)=\left\\{X=X_{1}+iX_{2}:X_{1},X_{2}\in
L^{0}(\Omega)\right\\}.$ A real-valued random variable $X\in L^{0}(\Omega)$ is
said to be _symmetric $\alpha$-stable_ if its characteristic function is of
the form

$\displaystyle\mathbb{E}\left[e^{isX}\right]=\exp\left\\{-\sigma^{\alpha}|s|^{\alpha}\right\\},$

where $\sigma>0$ is called the _scale parameter_ of $X$ and $\alpha\in(0,2]$
its _index of stability_. We write $X\sim S\alpha S(\sigma)$. In the
multivariate case, a real-valued symmetric $\alpha$-stable random vector
$\bm{X}=(X_{1},\dots,X_{n})$ is defined by its joint characteristic function

$\displaystyle\mathbb{E}\left[\exp\left\\{i(\bm{s},\bm{X})\right\\}\right]=\exp\left\\{-\int_{S^{n-1}}\left|(\bm{\theta},\bm{s})\right|^{\alpha}\Gamma(d\bm{\theta})\right\\},$

where $(\bm{x},\bm{y})$ denotes the scalar product of two vectors
$\bm{x},\bm{y}\in\mathbb{R}^{n}$, and $S^{n-1}$ is the unit sphere in
$\mathbb{R}^{n}$. The measure $\Gamma$ is called the _spectral measure_ of
$\bm{X}$. It is unique, finite and symmetric for $0<\alpha<2$ [17, Theorem
2.4.3]. A random variable $X\in L^{0}_{c}(\Omega)$ has _complex symmetric
$\alpha$-stable_ distribution if its real and imaginary parts form a $S\alpha
S$ random vector, i.e. if the vector $(Re(X),Im(X))$ is jointly $S\alpha S$.

To give a rigorous definition of harmonizable $S\alpha S$ processes, the
notion of complex random measures needs to be introduced. Let
$(E,\mathcal{E})$ be a measurable space, and let $(S^{1},\mathcal{B}(S^{1}))$
be the measurable space on the unit circle $S^{1}$ equipped with the Borel
$\sigma$-algebra $\mathcal{B}(S^{1})$. Let $k$ be a measure on the product
space $(E\times S^{1},\mathcal{E}\times\mathcal{B}(S^{1}))$, and let

$\mathcal{E}_{0}=\left\\{A\in\mathcal{E}:k(A\times S^{1})<\infty\right\\}.$

A _complex-valued $S\alpha S$ random measure_ on $(E,\mathcal{E})$ is an
independently scattered, $\sigma$-additive, complex-valued set function

$M:\mathcal{E}_{0}\rightarrow L^{0}_{c}(\Omega)$

such that the real and imaginary part of $M(A)$, i.e. the vector
$\left(Re(M(A)),Im(M(A))\right)$, is jointly $S\alpha S$ with spectral measure
$k(A\times\cdot)$ for every $A\in\mathcal{E}_{0}$ [17, Definition 6.1.2]. We
refer to $k$ as the _circular control measure_ of $M$, and denote by
$m(A)=k(A\times S^{1})$ the _control measure_ of $M$. Furthermore, $M$ is
isotropic if and only if its circular control measure is of the form

$k=m\cdot\gamma,$

where $\gamma$ is the uniform probability measure on $S^{1}$ [17, Example
6.1.6].

Define the space of $\alpha$-integrable functions on $E$ with respect to the
measure $m$ by

$L^{\alpha}(E,m)=\left\\{u:E\rightarrow\mathbb{C}:\int_{E}|u(x)|^{\alpha}m(dx)<\infty\right\\}.$

A stochastic integral with respect to a complex $S\alpha S$ random measures
$M$ is defined by

$\displaystyle I(u)=\int_{E}u(x)M(dx),$

for all $u\in L^{\alpha}(E,m)$. This stochastic integration is well-defined on
the space $L^{\alpha}(E,m)$ [17, Chapter 6.2]. In fact, for simple functions
$u(x)=\sum_{j=1}^{n}c_{j}\mathbbm{1}_{A_{j}}(x)$, where $c_{j}\in\mathbb{C}$
and $A_{j}\in\mathcal{E}_{0}$, it is easily seen that the integral $I(u)$ is
well-defined. Moreover, for any function $u\in L^{\alpha}(E,m)$ one can find a
sequence of simple functions $\\{u_{n}\\}_{n=1}^{\infty}$ which converges to
the function $u$ almost everywhere on $E$ with $|u_{n}(x)|\leq v(x)$ for some
function $v\in L^{\alpha}(E,m)$ for all $n\in\mathbb{N}$ and $x\in E$. The
sequence $\\{I(u_{n})\\}_{n=1}^{\infty}$ converges in probability, and $I(u)$
is then defined as its limit.

Setting $(E,\mathcal{E})=(\mathbb{R},\mathcal{B})$ and $u(t,x)=e^{itx}$, the
definition of harmonizable $S\alpha S$ processes is as follows.

######  Definition 1.

The stochastic process $X=\left\\{X(t):t\in\mathbb{R}\right\\}$ defined by

$\displaystyle X(t)=Re\left(\int_{\mathbb{R}}e^{itx}M(dx)\right),$

where $M$ is a complex $S\alpha S$ random measure on
$(\mathbb{R},\mathcal{B})$ with finite circular control measure $k$
(equivalently, with finite control measure $m$), is called a _real
harmonizable $S\alpha S$ process_.

A real harmonizable $S\alpha S$ process is stationary if and only if $M$ is
isotropic, i.e. its spectral measure is of the form $k=m\cdot\gamma$ [17,
Theorem 6.5.1]. In this case $X$ is called a _stationary real harmonizable
$S\alpha S$ process_. Furthermore, by [17, Proposition 6.6.3] the finite-
dimensional characteristic function of a SRH $S\alpha S$ process $X$ is given
by

$\displaystyle\mathbb{E}\left[\exp\left\\{i\sum_{i=1}^{n}s_{i}X(t_{i})\right\\}\right]=\exp\left\\{-\lambda_{\alpha}\int_{\mathbb{R}}\left|\sum_{j,k=1}^{n}s_{j}s_{k}\cos\left(\left(t_{k}-t_{j}\right)x\right)\right|^{\alpha/2}m(dx)\right\\}$
(1)

with constant
$\lambda_{\alpha}=\frac{1}{2\pi}\int_{0}^{2\pi}|\cos\left(x\right)|^{\alpha}dx$
$s_{1},\dots,s_{n}\in\mathbb{R}$ and $t_{1},\dots,t_{n}\in\mathbb{R}$ for all
$n\in\mathbb{N}$. Clearly, the process $X$ is uniquely characterized by its
control measure $m$ since all its finite-dimensional distributions are
determined by $m$.

We assume that the control measure $m$ is absolutely continuous with respect
to the Lebesgue measure on $\mathbb{R}$ with symmetric density function $f$,
which we refer to as the _spectral density_ of the process $X$. The goal is to
estimate $f$ from one single sample of observations
$\left(X(t_{1}),\dots,X(t_{n})\right)$.

###  2.1 Ergodicity and series representation of stationary real harmonizable
$S\alpha S$ processes

In statistical physics the study of ergodic properties of random processes is
motivated by the fundamental question whether long-term empirical observations
of a random process evolving in time, e.g. the motion of gas molecules,
suffice to estimate the mean value of the observable on the state space of the
process. In other words, are time-averages of observables equal to their so-
called phase or ensemble averages? First groundbreaking results, i.e. the
pointwise or strong ergodic theorem and mean ergodic theorem were proven by
Birkhoff in 1931 and von Neumann in 1932, respectively. These results also
sparked the general study of ergodic theory in mathematics, in particular the
field of dynamical systems.

In practical applications, a finite sample of the SRH $S\alpha S$ process $X$
at discrete time points will be the basis for all statistical evaluations.
Hence, in the following we will state the ergodic theorems for sequences of
random variables only. We first give a brief introduction to ergodic theory,
in particular we mention the notions of measure preserving transforms and
their invariant sets.

######  Definition 2.

  1. (a)

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space. The map
$T:\Omega\rightarrow\Omega$ is called _measure preserving transformation_ if
it is measurable with $\mathbb{P}(T^{-1}A)=\mathbb{P}(A)$ for all
$A\in\mathcal{F}$.

  2. (b)

A set $A\in\mathcal{F}$ is called _invariant_ with respect to $T$ if
$T^{-1}A=A$.

######  Lemma 1.

Let $\left\\{X_{k}\right\\}_{k=0}^{\infty}\subset L^{0}(\Omega)$ be a
stationary sequence of random variables. Then, there exists a measure
preserving transform $T$ such that $X_{k}(\omega)\stackrel{{\scriptstyle
d}}{{=}}X_{0}\left(T^{k}\omega\right)$ for all $k\geq 0$ and
$\omega\in\Omega$. We say that $\left\\{X_{k}\right\\}_{k=0}^{\infty}$ is
generated by $X_{0}$ and $T$.

Denote by $\mathcal{I}$ the $\sigma$-algebra of invariant sets with respect to
$T$. We can now state the ergodic theorem of Birkhoff-Khinchin, which is the
probabilistic version of the pointwise ergodic theorem.

######  Theorem 1 (Ergodic Theorem, Birkhoff-Khinchin).

Let $\left\\{X_{k}\right\\}_{k=0}^{\infty}$ be a stationary sequence of random
variables generated by a random variable $X_{0}$ and a measure preserving
transform $T$. Then, the limit

$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{k=0}^{n-1}h(X_{k})=\mathbb{E}\left[h(X_{0})|\mathcal{I}\right]$

exists almost surely, if $h:\mathbb{R}\rightarrow\mathbb{R}$ is a measurable
function with $\mathbb{E}\left[|h(X_{0})|\right]<\infty$.

######  Definition 3.

A process $X$ is ergodic if and only if the $\sigma$-algebra of invariant sets
is trivial.

As a consequence, the limit in the above ergodic theorem is deterministic with
$\mathbb{E}\left[h(X_{0})|\mathcal{I}\right]=\mathbb{E}\left[h(X_{0})\right]$.
Note that in the time-continuous case the arithmetic mean in Theorem 1 is
replaced by $\frac{1}{T}\int_{0}^{T}h(X(\tau))d\tau$.

######  Proposition 1.

Harmonizable $S\alpha S$ process with $0<\alpha<2$ are non-ergodic.

This result was proven by Cambanis, Hardin Jr., Weron [3, Theorem 4]. The
authors also show that $S\alpha S$ moving averages are ergodic and Sub-
Gaussian processes are non-ergodic, see [17, Chapter 3.6, 3.7] for the
respective definitions of these processes.

###  2.2 Independent path realizations, ensemble averages and $\alpha$-sine
transform

In the case that samples of independent path realizations of a SRH $S\alpha S$
process $X$ are available, consistent estimation of a symmetric spectral
density $f$ of $X$ can be performed without any problems caused by the non-
ergodicity of the process. Consider the codifference function, which is the
stable law analog of the covariance function when second moments do not exist.
It is defined by
$\tau(t)=2\|X(0)\|_{\alpha}^{\alpha}-\|X(t)-X(0)\|_{\alpha}^{\alpha}$, where
$\|X(0)\|_{\alpha}$ and $\|X(t)-X(0)\|_{\alpha}$ are the scale parameters of
$X(0)$ and $X(t)-X(0)$, respectively, see [17, Chapter 2.10]. Using the
finite-dimensional characteristic function in Equation (1) these scale
parameters can be computed explicitly which yields

$\displaystyle\tau(t)=2\sigma^{\alpha}-2^{\alpha}\lambda_{\alpha}\int_{\mathbb{R}}\left|\sin\left(\frac{tx}{2}\right)\right|^{\alpha}m(dx),$
(2)

see [17, Chapter, 6.7] and [8, Section 5.2]. Following Equation (2) and
assuming the control measure $m$ to have a symmetric density function $f$ we
define the $\alpha$-sine transform of $f$ as

$\displaystyle
T_{\alpha}f(t):=\int_{0}^{\infty}\left|\sin\left(tx\right)\right|^{\alpha}f(x)dx=\frac{2\sigma^{\alpha}-\tau(2t)}{2^{\alpha+1}\lambda_{\alpha}}=\frac{\|X(2t)-X(0)\|_{\alpha}^{\alpha}}{2^{\alpha+1}\lambda_{\alpha}}.$
(3)

Consistent estimation of the right hand side of (3) (based on independent
realizations of $X$) and the inversion of the above integral transform, which
was studied in [8], yield an estimate of the spectral density function $f$.

In more detail, let $X^{(1)},\dots,X^{(L)}$ denote a sample of independent
realizations of the path of $X$ finitely observed at points
$0=t_{0}<t_{1}<\dots<t_{n}<\infty$, i.e.
$X^{(l)}=(X^{(l)}(t_{0}),\dots,X^{(l)}(t_{n}))$ for $l=1,\dots,L$. Then, for
any fixed time instant $t_{i}$ the sample
$(X^{(1)}(t_{i}),\dots,X^{(L)}(t_{i}))$ consists of i.i.d. $S\alpha S(\sigma)$
random variables. In particular, the vector of differences
$(X^{(1)}(t_{i})-X^{(1)}(t_{0}),\dots,X^{(L)}(t_{i})-X^{(L)}(t_{0}))$ consists
of i.i.d. random samples of $X(t_{i})-X(t_{0})\sim S\alpha S(\sigma(t_{i}))$
with scale parameter $\sigma(t_{i})=\|X(t_{i})-X(t_{0})\|_{\alpha}$ which
depends on the lag $t_{i}-t_{0}$.

There are a handful of consistent parameter estimation techniques available in
literature, e.g. McCulloch’s quantile based method [11] or the regression-type
estimators by Koutrouvelis [9]. These methods can be used to estimate the
index of stability $\alpha$ and the scale parameters $\sigma(t_{i})$. Let
$(\hat{\alpha},\hat{\sigma}(t_{i}))$ be a consistent estimator of
$(\alpha,\sigma(t_{i}))$. Then,

$\displaystyle\widehat{T}_{\alpha}f(t_{i}/2)=\frac{\hat{\sigma}(t_{i})}{2^{\hat{\alpha}+1}\lambda_{\hat{\alpha}}}$
(4)

is consistent for $T_{\alpha}f(t_{i}/2)$ for all $i=1,\dots,n$. Using the
results of [8] we can estimate the Fourier transform of the spectral density
$f$ at equidistant points, which then allows us to reconstruct $f$ itself
using interpolation methods from sampling theory and Fourier inversion. An
illustration of the described estimation method for independent paths is given
in Section 5.1, Figure 1. We will not go into further details at this point as
our main interest lies within single path statistics.

##  3 Non-ergodic limit of sample functions

Trying to estimate the spectral density $f$ from a single path of a SRH
$S\alpha S$ process, determining the non-ergodic almost sure limit of
empirical functions is essential. We make use of the following LePage type
series representation of $X$ which stems from the series representation of
complex $S\alpha S$ random measures [17, Section 6.4]. As a consequence, SRH
$S\alpha S$ processes are in fact conditionally Gaussian which allows us to
use their underlying Gaussian structure for further analysis.

######  Proposition 2 (LePage Series representation).

Let $X=\\{X(t):t\in\mathbb{R}\\}$ be a SRH $S\alpha S$ process with finite
control measure $m$. Then, X is conditionally stationary centered Gaussian
with

$\displaystyle X(t)\stackrel{{\scriptstyle
d}}{{=}}\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}G(t),\quad
t\in\mathbb{R},$ (5)

where

$\displaystyle
G(t)=\sum_{k=1}^{\infty}\Gamma_{k}^{-1/\alpha}\left(G_{k}^{(1)}\cos(tZ_{k})+G_{k}^{(2)}\sin(tZ_{k})\right),\quad
t\in\mathbb{R},$

and the constants $b_{\alpha}$ and $C_{\alpha}$ are given by

$\displaystyle b_{\alpha}$
$\displaystyle=2^{\alpha/2}\Gamma\left(1+\alpha/2\right),$ $\displaystyle
C_{\alpha}$
$\displaystyle=\left(\int_{0}^{\infty}x^{-\alpha}\sin(x)dx\right)^{-1}=\begin{cases}\frac{1-\alpha}{\Gamma(2-\alpha)\cos(\pi\alpha/2)},&\alpha\neq
1,\\\ 2/\pi,&\alpha=1.\end{cases}$

Furthermore, the sequence $\left\\{\Gamma_{k}\right\\}_{k=1}^{\infty}$ denotes
the arrival times of a unit rate Poisson process,
$\\{G_{k}^{(i)}\\}_{k=1}^{\infty}$, $i=1,2$, are sequences of i.i.d. standard
normally distributed random variables, and
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$ is a sequence of i.i.d. random
variables with law $m(\cdot)/m(\mathbb{R})$.

######  Assumption.

Additionally to the assumption that the control measure $m$ is absolutely
continuous with respect to the Lebesgue measure on $\mathbb{R}$, let its
density $f$ be symmetric with $m(\mathbb{R})=1$, i.e. the spectral density $f$
is in fact a probability density function.

The process $G$ in Proposition 2, conditionally on the sequences
$\left\\{\Gamma_{k}\right\\}_{k=1}^{\infty}$ and
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$, is a stationary centered Gaussian
process with autocovariance function

$\displaystyle
r(t)=\mathbb{E}\left[G(0)G(t)\right]=\sum_{k=1}^{\infty}\Gamma_{k}^{-2/\alpha}\cos\left(tZ_{k}\right),\quad
t\in\mathbb{R}$

[17, Proposition 6.6.4]. Note that by the Wiener-Khinchin theorem the
autocovariance function $r$ and the spectral measure $\zeta$ of $G$ are
directly related by the Fourier transform, i.e.
$r(t)=\int_{\mathbb{R}}e^{it\omega}\zeta(d\omega)$, and one can compute

$\displaystyle\zeta(d\omega)=\sum_{k=1}^{\infty}\frac{\Gamma_{j}^{-2/\alpha}}{2}\left(\delta\left(d\omega-
Z_{j}\right)+\delta\left(d\omega+Z_{j}\right)\right)$

via Fourier inversion. The spectral measure $\zeta$ is purely discrete, hence
$G$ is non-ergodic, as Gaussian processes are ergodic if and only if their
corresponding spectral measure is absolutely continuous [10, Chapter 6.5].

The asymptotic behavior of time averages for non-ergodic Gaussian processes is
studied in [21]. There, the so-called Harmonic Gaussian processes of the form

$\displaystyle X(t)=\sum_{k}A_{k}e^{it\omega_{k}},$ (6)

are considered, where $A_{k}$ are independent complex Gaussian random
variables, and $\omega_{k}$ are the frequencies of the process. The
autocovariance function $r$ and spectral measure $\sigma$ of the above
harmonic Gaussian process are given by

$\displaystyle r(t)$
$\displaystyle=\sum_{k}\mathbb{E}|A_{k}|^{2}\cos(\omega_{k}t),$
$\displaystyle\sigma(d\omega)$
$\displaystyle=\sum_{k}\frac{\mathbb{E}|A_{k}|^{2}}{2}\left(\delta\left(d\omega-\omega_{k}\right)+\delta\left(d\omega+\omega_{k}\right)\right)$

[21, Equation (22), (23)]. Since SRH $S\alpha S$ processes are conditionally
Gaussian, we can embed them into the above non-ergodic harmonic Gaussian case
and apply the results in [21].

######  Theorem 2.

Let $X$ be a SRH $S\alpha S$ process.

  1. (a)

The process $X$ admits the series representation

$\displaystyle X(t)\stackrel{{\scriptstyle
d}}{{=}}\sum_{k=1}^{\infty}R_{k}\cos\left(\Theta_{k}+tZ_{k}\right),$ (7)

where $\left\\{\Theta_{k}\right\\}$ are i.i.d. uniformly distributed on
$(0,2\pi)$,

$\displaystyle
R_{k}=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\Gamma_{k}^{-1/\alpha}\sqrt{\left(G_{k}^{(1)}\right)^{2}+\left(G_{k}^{(2)}\right)^{2}}$

and
$\\{\Gamma_{k}\\}_{k=1}^{\infty},\\{G^{(1)}_{k}\\}_{k=1}^{\infty},\\{G^{(2)}_{k}\\}_{k=1}^{\infty},\\{Z_{k}\\}_{k=1}^{\infty}$
are the sequences in the series representation of $X$ from Proposition 2.

  2. (b)

The (time-average) empirical characteristic function of the summands
$\left(R_{k}\cos\left(\Theta_{k}+tZ_{k}\right)\right)_{k\in\mathbb{N}}$
converges almost surely, i.e. for all
$\bm{\lambda}=(\lambda_{k})_{k\in\mathbb{N}}$ it holds

$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\exp\left\\{i\sum_{k=1}^{\infty}\lambda_{k}R_{k}\cos\left(\Theta_{k}+\tau
Z_{k}\right)\right\\}d\tau=\prod_{k=1}^{\infty}J_{0}(\lambda_{k}R_{k})\quad
a.s.,$

where $J_{0}(s)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{is\cos(u)}du$ is the Bessel
function of the first kind of order $0$.

###### Proof.

  1. (a)

The series representation of $X$ in Equation (5) clearly shows a strong
resemblance to a Fourier series. Similarly to the equivalence of the sine-
cosine and exponential form of Fourier series, the process X admits the form

$\displaystyle X(t)\stackrel{{\scriptstyle
d}}{{=}}\sum_{k=-\infty}^{\infty}A_{k}e^{itZ_{k}},$ (8)

where $\\{A_{k}\\}_{k=1}^{\infty}$ conditionally on
$\left\\{\Gamma_{k}\right\\}_{k=1}^{\infty}$ is a sequence of complex Gaussian
random variables with

$\displaystyle
A_{k}=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\frac{\Gamma_{k}^{-1/\alpha}}{2}G_{k}^{(1)}-i\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\frac{\Gamma_{k}^{-1/\alpha}}{2}G_{k}^{(2)}$
(9)

for $k\geq 1$, and $A_{0}=0$. Furthermore, we set $A_{-k}=A_{k}^{\ast}$, the
complex conjugate of $A_{k}$, and $Z_{-k}=-Z_{k}$. One can easily verify that
the exponential form in Equation (8) is indeed equal in distribution to the
series representation (5) of $X$ from Proposition 2. To do so, denote
$\tilde{C}_{k}=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\frac{\Gamma_{k}^{-1/\alpha}}{2}$
for ease of notation, and compute

$\displaystyle A_{k}e^{itZ_{k}}$
$\displaystyle=\tilde{C}_{k}\left(G_{k}^{(1)}+iG_{k}^{(2)}\right)\left(\cos(tZ_{k})+i\sin(tZ_{k})\right)$
$\displaystyle=\tilde{C}_{k}\left(G_{k}^{(1)}\cos(tZ_{k})-G_{k}^{(2)}\sin(tZ_{k})+iG_{k}^{(1)}\sin(tZ_{k})+iG_{k}^{(2)}\cos(tZ_{k})\right),$
$\displaystyle A_{-k}e^{-Z_{-k}}$
$\displaystyle=\tilde{C}_{k}\left(G_{k}^{(1)}-iG_{k}^{(2)}\right)\left(\cos(tZ_{k})-i\sin(tZ_{k})\right)$
$\displaystyle=\tilde{C}_{k}\left(G_{k}^{(1)}\cos(tZ_{k})-G_{k}^{(2)}\sin(tZ_{k})-iG_{k}^{(1)}\sin(tZ_{k})-iG_{k}^{(2)}\cos(tZ_{k})\right).$

Consequently,

$\displaystyle A_{k}e^{itZ_{k}}+A_{-k}e^{-Z_{-k}}$
$\displaystyle=2\tilde{C}_{k}\left(G_{k}^{(1)}\cos(tZ_{k})-G_{k}^{(2)}\sin(tZ_{k})\right)$
$\displaystyle=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\Gamma_{k}^{-1/\alpha}\left(G_{k}^{(1)}\cos(tZ_{k})-G_{k}^{(2)}\sin(tZ_{k})\right)$

is equal in distribution to the $k$-th summand of the series representation in
Proposition 2, taking into account that $-G_{k}^{(2)}\stackrel{{\scriptstyle
d}}{{=}}G_{k}^{(2)}$ by symmetry of the standard normal distribution.

Conditionally on $\left\\{\Gamma_{k}\right\\}_{k=1}^{\infty}$ and
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$, the SRH $S\alpha S$ process X is in
fact a non-ergodic harmonic Gaussian process as in Equation (6). Analogously
to [21, Equation (27)], setting

$\displaystyle R_{k}$
$\displaystyle=2|A_{k}|=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\Gamma_{k}^{-1/\alpha}\sqrt{\left(G_{k}^{(1)}\right)^{2}+\left(G_{k}^{(1)}\right)^{2}}$

yields an alternative series representation for the process $X$ with

$\displaystyle X(t)\stackrel{{\scriptstyle
d}}{{=}}\sum_{k=1}^{\infty}R_{k}\cos\left(\Theta_{k}+tZ_{k}\right),$

where $R_{k}$ are the amplitudes of the spectral points at the frequencies
$Z_{k}$, and $\Theta_{k}$ are i.i.d. uniformly distributed phases of $A_{k}$
on $(0,2\pi)$.

  2. (b)

We first compute the characteristic function of
$\left(R_{k}\cos\left(\Theta_{k}+tZ_{k}\right)\right)_{k\in\mathbb{N}}$
conditional on $\\{R_{k}\\}_{k=1}^{\infty}$ and $\\{Z_{k}\\}_{k=1}^{\infty}$:

$\displaystyle\mathbb{E}\left[\exp\left(i\sum_{k=1}^{\infty}\lambda_{k}R_{k}\cos\left(\Theta_{k}+tZ_{k}\right)\right)|\left\\{R_{k}\right\\},\left\\{Z_{k}\right\\}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod_{k=1}^{\infty}\exp\left(i\lambda_{k}R_{k}\cos\left(\Theta_{k}+tZ_{k}\right)\right)|\left\\{R_{k}\right\\},\left\\{Z_{k}\right\\}\right]$
$\displaystyle=$
$\displaystyle\prod_{k=1}^{\infty}\underbrace{\mathbb{E}\left[\exp\left(i\lambda_{k}R_{k}\cos\left(\Theta_{k}+tZ_{k}\right)\right)|\left\\{R_{k}\right\\},\left\\{Z_{k}\right\\}\right]}_{\frac{1}{2\pi}\int_{-\pi}^{\pi}\exp\left(i\lambda_{k}R_{k}\cos(u)\right)du=J_{0}\left(\lambda_{k}R_{k}\right)}=\prod_{k=1}^{\infty}J_{0}(\lambda_{k}R_{k}).$

Note that the the infinite product and expectation in the above can be
interchanged by Lebesgue’s dominated convergence theorem and the fact that
$\Theta_{k}$ are i.i.d. Furthermore, for the computation of the Bessel
function we use the fact that $\Theta_{k}+c$ modulo $2\pi$ stays uniformly
distributed on $(0,2\pi)$ for any constant $c\in\mathbb{R}$.

Next, we consider the time-average characteristic function and show that it
converges almost surely to the above conditional characteristic function. In
the following, we will use the Jacobi-Anger identity
$e^{is\cos(u)}=\sum_{m\in\mathbb{Z}}i^{m}e^{imu}J_{m}(s)$, where $J_{m}$ is
the Bessel function of the first kind of order $m$. Furthermore, we denote by
$M_{N}$ the family of subsets $S\subset\mathbb{Z}$ with cardinality $N$. The
computation is similar to the one in [21, Proof of Theorem 5]:

$\displaystyle I_{T}$
$\displaystyle=\int_{0}^{T}\exp\left\\{i\sum_{k=1}^{\infty}\lambda_{k}R_{k}\cos\left(\Theta_{k}+\tau
Z_{k}\right)\right\\}d\tau$
$\displaystyle=\int_{0}^{T}\lim_{N\rightarrow\infty}\prod_{k=1}^{N}\exp\left\\{i\lambda_{k}R_{k}\cos\left(\Theta_{k}+\tau
Z_{k}\right)\right\\}d\tau$
$\displaystyle=\lim_{N\rightarrow\infty}\int_{0}^{T}\prod_{k=1}^{N}\sum_{m\in\mathbb{Z}}i^{m}e^{im\left(\Theta_{k}+\tau
Z_{k}\right)}J_{m}(\lambda_{k}R_{k})d\tau$
$\displaystyle=\lim_{N\rightarrow\infty}\int_{0}^{T}\sum_{S\in
M_{N}}\prod_{m_{k}\in S}i^{m_{k}}e^{im_{k}\left(\Theta_{k}+\tau
Z_{k}\right)}J_{m_{k}}(\lambda_{k}R_{k})d\tau$
$\displaystyle=\lim_{N\rightarrow\infty}\sum_{S\in
M_{N}}\exp\left\\{i\sum_{m_{k}\in
S}m_{k}\Theta_{k}\right\\}\left(\prod_{m_{k}\in
S}i^{m_{k}}J_{m_{k}}(\lambda_{k}R_{k})\right)\underbrace{\int_{0}^{T}\exp\Bigg{\\{}i\tau\underbrace{\sum_{m_{k}\in
S}m_{k}Z_{k}}_{=\Omega_{S}}\Bigg{\\}}d\tau}_{=\widetilde{I_{T}}(\Omega_{S})}.$
(10)

It holds that $\frac{1}{T}\widetilde{I_{T}}(\Omega_{S})$ converges almost
surely to $0$ as $T\rightarrow\infty$ for all $\Omega_{S}\neq 0$ [21, Lemma
1]. Since the $Z_{k}$ are i.i.d. random variables drawn from the absolutely
continuous distribution with probability density function $f$ and are
therefore almost surely not multiples of each other, the summands in Equation
((b)) above vanish for all $S\neq(0,0,\dots)$. In the case $S=(0,0,\dots)$ it
holds that $\frac{1}{T}\widetilde{I_{T}}(0)$ converges almost surely to 1.
Thus,

$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}I_{T}=\prod_{k=1}^{\infty}J_{0}(\lambda_{k}R_{k})$

almost surely. Taking limits and integrating can be interchanged by Lebesgue’s
dominated convergence theorem since all functions are bounded by $1$.

∎

As a direct consequence of the series representation of a SRH $S\alpha S$
process $X$ and the ergodic Theorem 2 (b) above we can compute the almost sure
limit of the characteristic functions of $X$ and the lag process
$\left\\{X(t+h)-X(t):t\in\mathbb{R}\right\\}$, $h>0$.

######  Corollary 1.

  1. (a)

The empirical characteristic function $\frac{1}{T}\int_{0}^{T}e^{i\lambda
X(\tau)}d\tau$ converges almost surely to $\prod_{k=1}^{\infty}J_{0}(\lambda
R_{k})$ as $T\rightarrow\infty$.

  2. (b)

The empirical characteristic function of the lag process
$\left\\{X(t+h)-X(t):t\in\mathbb{R}\right\\}$, $h>0$, i.e.
$\frac{1}{T}\int_{0}^{T}e^{i\lambda X(\tau+h)-X(\tau)}d\tau$, converges almost
surely to $\prod_{k=1}^{\infty}J_{0}\left(2\lambda
R_{k}\sin\left(hZ_{K}/2\right)\right)$ as $T\rightarrow\infty$.

###### Proof.

  1. (a)

The almost sure limit of the empirical characteristic function follows
immediately from Theorem 2 (b).

  2. (b)

The computation of the almost sure limit
$\mathbb{E}\left[\exp\left\\{i\lambda(X(t+h)-X(t))\right\\}|\left\\{R_{k}\right\\},\left\\{Z_{k}\right\\}\right]$
is straightforward using
$\cos(\Theta_{k}+(t+h)Z_{k})-\cos(\Theta_{k}+tZ_{k})=-2\sin\left(hZ_{k}/2\right)\sin\left(\Theta_{k}+tZ_{k}+hZ_{k}/2\right)$,
where the latter $-\sin\left(\Theta_{k}+tZ_{k}+hZ_{k}/2\right)$ is equal in
distribution to $\cos\left(\Theta_{k}\right)$ with $\Theta_{k}\sim U(0,2\pi)$
i.i.d.

∎

######  Corollary 2.

  1. (a)

Let $h:\mathbb{R}\rightarrow\mathbb{R}$ be a $p$-periodic function with
absolutely convergent Fourier series
$h(x)=\sum_{n\in\mathbb{Z}}c_{n}e^{i\frac{p}{2\pi}nx}$. Then,

$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}h(X(\tau))d\tau=\sum_{n\in\mathbb{Z}}c_{n}\prod_{k=1}^{\infty}J_{0}\left(n\frac{p}{2\pi}R_{k}\right)\quad
a.s.$

  2. (b)

Let $h:\mathbb{R}\rightarrow\mathbb{R}$ be a function in $L^{1}(\mathbb{R})$
with integrable Fourier transform $\mathcal{F}h\in L^{1}(\mathbb{R})$ such
that $h(x)=\int_{\mathbb{R}}\mathcal{F}h(y)e^{ixy}dy$. Then,

$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}h(X(\tau))d\tau=\int_{\mathbb{R}}\mathcal{F}h(y)\left(\prod_{k=1}^{\infty}J_{0}(yR_{k})\right)dy\quad
a.s.$

###### Proof.

We prove part (b) only since all arguments are applicable for part (a) as
well. By Fubini’s theorem and Lebesgue’s dominated convergence theorem we can
interchange the order of integration and limit, which yields

$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}h(X(\tau))d\tau$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\int_{\mathbb{R}}\mathcal{F}h(y)\exp\left\\{iy\sum_{k=1}^{\infty}R_{k}\cos\left(\Theta_{k}+\tau
Z_{k}\right)\right\\}dyd\tau$
$\displaystyle=\int_{\mathbb{R}}\mathcal{F}h(y)\left(\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\exp\left\\{iy\sum_{k=1}^{\infty}R_{k}\cos\left(\Theta_{k}+\tau
Z_{k}\right)\right\\}d\tau\right)dy$
$\displaystyle=\int_{\mathbb{R}}\mathcal{F}h(y)\left(\prod_{k=1}^{\infty}J_{0}(yR_{k})\right)dy\quad
a.s.$

∎

Summarizing, it becomes clear that, conditionally on
$\left\\{R_{k}\right\\}_{k=1}^{\infty}$ and
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$, SRH $S\alpha S$ processes are non-
ergodic harmonic Gaussian processes from Equation (6). For SRH $S\alpha S$
processes additional randomness is introduced by the random frequencies
$\\{Z_{k}\\}_{k=1}^{\infty}$ and the arrival times
$\\{\Gamma_{k}\\}_{k=1}^{\infty}$, which are in a sense random variances to
the Gaussian random variables $\\{G_{k}^{(i)}\\}_{k=1}^{\infty}$, $i=1,2$.
Indeed, the sequence $\left\\{A_{k}\right\\}_{k=1}^{\infty}$ in Equation (9)
consists of variance mixtures of Gaussian random variables.

##  4 Spectral density estimation

The underlying Gaussian structure of the SRH $S\alpha S$ process $X$ plays a
central role in the estimation of the spectral density. The LePage series
representation and Theorem 2 demonstrate that the randomness of the paths of
$X$ is generated by the Gaussian variance mixture amplitudes
$\left\\{R_{k}\right\\}_{k=1}^{\infty}$ as well as the random frequencies
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$. Although these quantities are
inherently random and not directly observable, they are fixed for a given
path. It is therefore possible to consider a path of the harmonizable process
$X$ as a signal generated by frequencies
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$ with amplitudes
$\left\\{R_{k}\right\\}_{k=1}^{\infty}$ and phases
$\left\\{U_{k}\right\\}_{k=1}^{\infty}$ (see Equation (7)), and use standard
frequency estimation techniques to estimate
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$.

Recall that $Z_{1},Z_{2},\dots$ are i.i.d. with probability density function
$f$, which is the spectral density function of the process $X$. A density
kernel estimate will then yield the desired estimate of $f$. As the goal is to
estimate the spectral density $f$ alone, we can neglect the estimation of
amplitudes $\left\\{R_{k}\right\\}_{k=1}^{\infty}$ and phases
$\left\\{U_{k}\right\\}_{k=1}^{\infty}$.

###  4.1 Periodogram method

Spectral density estimation is a widely studied subject in the field of signal
processing, and we will rely on classical techniques provided there [2, 6,
13]. In particular, the periodogram is a standard estimate of the density
function $f_{\zeta}$ of an absolutely continuous spectral measure $\zeta$,
often also referred to as power spectral density. When the spectral measure of
a signal is purely discrete, i.e. the signal is produced by sinusoidal waves
alone, the periodogram still proves to be a powerful tool for the estimation
of frequencies [13]. Such a sinusoidal signal
$X=\left\\{X(t):t\in\mathbb{R}\right\\}$ is modeled by

$\displaystyle
X(t)=\mu+\sum\limits_{k=1}^{\infty}r_{k}\cos\left(\phi_{k}+\lambda_{k}t\right)+Y(t),$
(11)

where $\mu$ is called the overall mean of the signal and
$r_{k},\phi_{k},\lambda_{k}$ are the amplitude, phase and frequency of the
$k$-th sinusoidal. The process $Y=\left\\{Y(t):t\in\mathbb{R}\right\\}$ is a
noise process usually assumed to be a stationary zero mean process with
spectral density $f_{y}$ [13, p.5].

Note that in the above model the amplitudes, phases and frequencies are
deterministic and randomness is introduced by the noise process $Y$. For
harmonic Gaussian processes it is assumed that the phases are i.i.d. uniformly
distributed on $(0,2\pi)$ and the amplitudes are given by $R_{k}=|A_{k}|$,
where $A_{k}$ are complex Gaussian random variables as seen in [21]. For SRH
$S\alpha S$ process these Gaussian random variables are replaced by variance
mixtures of Gaussian variables, and, additionally, the frequencies are i.i.d.
random variables with density $f$. The amplitudes
$\left\\{R_{k}\right\\}_{k=1}^{\infty}$ and frequencies
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$ in (7) are fixed for a given path of
the SRH $S\alpha S$ process $X$. Furthermore, we have $\mu=0$. Also, the
phases $\left\\{U_{k}\right\\}_{k=1}^{\infty}$ do not play any role in the
following computation of the periodogram.

We define the discrete Fourier transform and the periodogram as follows.

######  Definition 4.

Let $x=\left(x(1),\dots,x(n)\right)\in\mathbb{C}^{n}$ be a (complex) vector.
Then, the transform

$\displaystyle
F_{n}\left(\theta_{j}\right)=n^{-1}\sum\limits_{k=1}^{n}x(k)e^{-ik\theta_{j}}$

defines the _discrete Fourier transform_ of $x$ at the Fourier frequencies
$\\{\theta_{j}\\}_{j=1}^{n}=\left\\{2\pi j/n\right\\}_{j=1}^{n}$. Furthermore,
the _periodogram_ of $x$ at the Fourier frequencies
$\\{\theta_{j}\\}_{j=1}^{n}$ is defined by

$\displaystyle
I_{n}\left(\theta_{j}\right)=\left|F_{n}\left(\theta_{j}\right)\right|^{2}=\left|n^{-1}\sum\limits_{k=1}^{n}x(k)e^{-ik\theta_{j}}\right|^{2}.$

Since the discrete Fourier transform can be viewed as a discretization of the
continuous time Fourier transform, the definition of the periodogram can
simply be extended to arbitrary frequencies by replacing the Fourier
frequencies $\theta_{j}$ with $\theta>0$ in the above. The periodogram is a
standard tool for the estimation of the spectral density if the spectral
measure is absolutely continuous with respect to the Lebesgue measure on
$\mathbb{R}$. In the case that the spectral measure is purely discrete it can
be used to estimate the underlying frequencies of a signal in the following
way.

Let $\left(x(1),\dots x(n)\right)=\left(X(t_{1}),\dots,X(t_{n})\right)$ be a
sample generated from the SRH $S\alpha S$ process
$X=\left\\{X(t):t\in\mathbb{R}\right\\}$ at equidistant points
$t_{j}=j\delta$, $\delta>0$, for $j=1,\dots,n$. Then, according to Equation
(7) the sample $x(j)$ is of the form

$\displaystyle x(j)=\sum\limits_{l=1}^{\infty}R_{l}\cos\left(U_{l}+j\delta
Z_{l}\right)$

for all $j=1,\dots,n$. For $\theta>0$, the discrete Fourier transform is given
by

$\displaystyle F_{n}(\theta\delta)$
$\displaystyle=n^{-1}\sum\limits_{j=1}^{n}x\left(j\right)e^{-ij\theta\delta}=n^{-1}\sum\limits_{j=1}^{n}\sum\limits_{l=1}^{\infty}R_{l}\cos\left(U_{l}+j\delta
Z_{l}\right)e^{-ij\theta\delta}$ (12)
$\displaystyle=n^{-1}\sum\limits_{l=1}^{\infty}R_{l}\sum\limits_{j=1}^{n}e^{-ij\theta\delta}\underbrace{\cos\left(U_{l}+j\delta
Z_{l}\right)}_{=\frac{e^{i\left(U_{l}+j\delta
Z_{l}\right)}+e^{-i\left(U_{l}+j\delta Z_{l}\right)}}{2}}$
$\displaystyle=n^{-1}\sum\limits_{l=1}^{\infty}\Bigg{(}\frac{e^{iU_{l}}R_{l}}{2}\sum\limits_{j=1}^{n}e^{ij\delta\left(Z_{l}-\theta\right)}+\frac{e^{-iU_{l}}R_{l}}{2}\sum\limits_{j=1}^{n}e^{-ij\delta\left(Z_{l}+\theta\right)}\Bigg{)}.$
(13)

Taking the squared modulus of the right-hand side expression in Equation (12)
yields the periodogram. For the full computation see Appendix A, and for more
details we refer to [13].

For any fixed $k\in\mathbb{N}$ and $\theta$ close to $|Z_{k}|$ the
periodogram’s behavior for increasing $n$ can be described by the following
expression.

$\displaystyle I_{n}(\theta\delta)$
$\displaystyle=\frac{R_{k}^{2}}{4}\Biggl{(}\underbrace{\frac{\sin^{2}\left(\frac{n}{2}\delta(Z_{k}-\theta)\right)}{n^{2}\sin^{2}\left(\frac{1}{2}\delta(Z_{k}-\theta)\right)}}_{(\ast)}+\underbrace{\frac{\sin^{2}\left(\frac{n}{2}\delta(Z_{k}+\theta)\right)}{n^{2}\sin^{2}\left(\frac{1}{2}\delta(Z_{k}+\theta)\right)}}_{(\ast\ast)}\Biggr{)}+O\left(n^{-2}\right).$
(14)

Analogously to [13, pp. 35-36], for $\theta>0$ and $Z_{k}>0$ the second term
$(\ast\ast)$ vanishes with rate $O(n^{-2})$ as $n\rightarrow\infty$.
Furthermore, the first term $(\ast)$ converges to $1$ as $\theta$ tends to
$Z_{k}$ or to $0$ uniformly with rate $O(n^{-2})$ for any $\theta>0$ with
$\left|\theta-Z_{k}\right|>\varepsilon$ and $\varepsilon>0$ [7, p. 515]. On
the other hand, for $\theta>0$ and $Z_{k}<0$ the first term $(\ast)$ in
Equation (14) behaves like $O(n^{-2})$ as $n\rightarrow\infty$ with uniform
convergence to $0$ as before, and $(\ast\ast)$ converges to $1$ as
$\theta\rightarrow|Z_{k}|$.

Therefore, the periodogram $I_{n}$ shows pronounced peaks close to the
absolute values of the true frequencies $\left\\{Z_{k}\right\\}$, and the
locations of those peaks will be used as estimates. Note that the amplitudes
are not of interest, hence we do not need to estimate
$\left\\{R_{k}\right\\}$.

######  Theorem 3.

Let $X=\left\\{X(t):t\in\mathbb{R}\right\\}$ be a SRH $S\alpha S$ process with
symmetric spectral density $f$, and
$\left(x(1),\dots,x(n)\right)=\left(X(t_{1}),\dots,X(t_{n})\right)$ be a
sample of the process $X$ at equidistant points $t_{1}<\dots<t_{n}$,
$n\in\mathbb{N}$ with $t_{j}=j\delta$, $\delta>0$ for $j=1,\dots,n$. The
frequency estimators $\left\\{\hat{Z}_{k,n}\right\\}_{k=1}^{N}$,
$N\in\mathbb{N}$, acquired from the locations of the $N$ largest peaks of the
periodogram of $(x(1),\dots,x(n))$, are strongly consistent for the absolute
frequencies $\left\\{|Z_{k}|\right\\}_{k=1}^{N}$, i.e.
$\hat{Z}_{k,n}\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}|Z_{k}|$ as
$n\rightarrow\infty$ for $k=1,\dots,N$.

###### Proof.

Recall that $\left(\Omega,\mathcal{F},P\right)$ denotes the probability space
on which the process $X$ lives. For any $\omega\in\Omega$ we denote the
corresponding path of $X$ by
$X(\omega)=\left\\{X(t;\omega):t\in\mathbb{R}\right\\}$. The periodogram of
the sample $\left(x(1),\dots,x(n)\right)$ at $\theta>0$ is the squared modulus
of the discrete Fourier transform at $\theta$, and from Equation (14) we know
that for $k=1,\dots,N$

$\displaystyle
I_{n}(\delta\theta;\omega)=\frac{R_{k}(\omega)^{2}}{4}\frac{\sin^{2}\left(\frac{n}{2}\delta(|Z_{k}(\omega)|-\theta)\right)}{n^{2}\sin^{2}\left(\frac{1}{2}\delta(|Z_{k}(\omega)|-\theta)\right)}+O\left(n^{-2}\right).$

The above converges to $R_{k}(\omega)^{2}/4$ as $\theta$ tends to
$\left|Z_{k}(\omega)\right|$ and to $0$ for
$\theta\neq\left|Z_{k}(\omega)\right|$. The convergence to $0$ is uniform for
$\left|\theta-|Z_{k}(\omega)|\right|>\varepsilon$ for any $\varepsilon>0$.
Similar to [7] it follows that the estimator $\hat{Z}_{k,n}(\omega)$ converges
to $\left|Z_{k}(\omega)\right|$ as $n\rightarrow\infty$ for any
$\omega\in\Omega$. To see this, assume that $\hat{Z}_{k,n}(\omega)$ does not
converge to $\left|Z_{k}(\omega)\right|$. Then, there exists a subsequence
$\hat{Z}_{k,n_{l}}$ which converges to some $Z^{\prime}\neq|Z_{k}(\omega)|$.
Along this subsequence $I_{n_{l}}(\hat{Z}_{k,n_{l}})$ converges to $0$ but
since the estimator $\hat{Z}_{k,n}$ is a local maximum of the periodogram
$I_{n}$ it holds that $I_{n_{l}}(\hat{Z}_{k,n_{l}})\geq
I_{n_{l}}(|Z_{k}(\omega)|)\rightarrow R_{k}(\omega)^{2}/4>0$, which leads to a
contradiction. ∎

###  4.2 Kernel density estimator and weak consistency

We denote by $\\{\hat{Z}_{k,n}\\}_{k=1}^{N}$ the estimators of the absolute
frequencies $\\{\left|Z_{k}\right|\\}_{k=1}^{N}$, where $N\in\mathbb{N}$ is
the number of the largest peaks located from the periodogram. Recall that
frequencies $Z_{k}$ are drawn independently from a distribution with the
sought-after symmetric spectral density $f$ as probability density function.

######  Remark 1.

Let $Z$ be a random variable with symmetric probability density function $f$.
Note that the probability density function of $|Z|$ satisfies the relation
$f_{|Z|}(x)=2f(x)$ for any $x\geq 0$.

Hence, applying a kernel density estimator on the estimates
$\\{\hat{Z}_{k,n}\\}_{k=1}^{N}$ yields an estimate of $2f$. Denote by
$\hat{f}_{N}$ the kernel density estimator with kernel function $\kappa$ and
bandwidth $h_{N}>0$:

$\displaystyle\hat{f}_{N}(x)=\frac{1}{Nh_{N}}\sum\limits_{k=1}^{N}\kappa\left(\frac{x-|Z_{k}|}{h_{N}}\right),$
(15)

and define the kernel density estimator $\hat{f}_{N,n}$ based on
$\left\\{\hat{Z}_{k,n}\right\\}_{k=1}^{N}$ by

$\displaystyle\hat{f}_{N,n}(x)=\frac{1}{Nh_{N}}\sum\limits_{k=1}^{N}\kappa\left(\frac{x-\hat{Z}_{k,n}}{h_{N}}\right).$
(16)

Kernel density estimators have been studied extensively, and the following
consistency results for $\hat{f}_{N}$ are well known, see [20, 22]. Let
$\|u\|_{\infty}=\sup_{x\in\mathbb{R}}|u(x)|$ denote the uniform norm of a
function $u$ on $\mathbb{R}$. Furthermore, let $\ \stackrel{{\scriptstyle
P}}{{\longrightarrow}}$ denote convergence in probability and $\
\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}$ almost sure convergence.

######  Lemma 2.

  1. (a)

If $\lim_{N\rightarrow\infty}Nh_{N}=\infty$, then $\hat{f}_{N}$ is weakly
pointwise consistent for $2f$ at every point of continuity $x$ of $f$, i.e.
$\hat{f}_{N}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}2f(x)$ as
$N\rightarrow\infty$.

  2. (b)

If $\lim_{N\rightarrow\infty}Nh_{N}^{2}=\infty$, then $\hat{f}_{N}$ is weakly
uniformly consistent for $2f$, i.e.
$\left\|\hat{f}_{N}-2f\right\|_{\infty}\stackrel{{\scriptstyle
P}}{{\longrightarrow}}0$.

  3. (c)

Let $\kappa$ be a right-continuous kernel function with bounded variation and
$\kappa(x)\rightarrow 0$ as $|x|\rightarrow\infty$, and let $f$ be uniformly
continuous. Choose the bandwidth $h_{N}>0$ such that
$\sum_{N=1}^{\infty}\exp(-\gamma Nh_{N}^{2})<\infty$ for any $\gamma>0$. Then,
the kernel density estimator $\hat{f}_{N}$ is strongly uniformly consistent
for $2f$, i.e. $\|\hat{f}_{N}-2f\|_{\infty}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}0$.

######  Theorem 4.

Let $X=\left\\{X(t):t\in\mathbb{R}\right\\}$ be a SRH $S\alpha S$ process with
symmetric spectral density $f$, and
$\left(x(1),\dots,x(n)\right)=\left(X(t_{1}),\dots,X(t_{n})\right)$ be a
sample of the process $X$ at equidistant points $t_{1}<\dots<t_{n}$,
$n\in\mathbb{N}$. Furthermore, let $\left\\{\hat{Z}_{k,n}\right\\}_{k=1}^{N}$,
$N\in\mathbb{N}$, be the strongly consistent frequency estimators of
$\left\\{|Z_{k}|\right\\}_{k=1}^{N}$ from Theorem 3. Consider the kernel
density estimator $\hat{f}_{N,n}$ defined in Equation (16).

  1. (a)

Let $\kappa$ be a continuous kernel function and $h=h_{N}>0$ such that
$\lim_{N\rightarrow\infty}Nh_{N}=\infty$. Then, at any continuity point $x\geq
0$ of $f$ the estimator $\hat{f}_{N,n}$ is weakly consistent for $2f$, i.e.
$\hat{f}_{N,n}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}2f(x)$ as
$n\rightarrow\infty$ and $N\rightarrow\infty$.

  2. (b)

Let $\kappa$ be a Lipschitz continuous kernel function and $h=h_{N}>0$ such
that $\lim_{N\rightarrow\infty}Nh_{N}^{2}=\infty$. Then, the estimator
$\hat{f}_{N,n}$ is weakly uniformly consistent for $2f$, i.e.
$\left\|\hat{f}_{N,n}-2f\right\|_{\infty}\stackrel{{\scriptstyle
P}}{{\longrightarrow}}0$ as $n\rightarrow\infty$ and $N\rightarrow\infty$.

  3. (c)

Under the conditions of Lemma 2 (c) the kernel density estimator
$\hat{f}_{N,n}(x)$ is strongly consistent for $2f(x)$ for any $x\in\mathbb{R}$
as $n\rightarrow\infty$ and $N\rightarrow\infty$. If the kernel function
$\kappa$ is additionally Lipschitz continuous, then $\hat{f}_{N,n}$ is
strongly uniformly consistent for $2f$, i.e.
$\left\|\hat{f}_{N,n}-2f\right\|_{\infty}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}0$ as $n\rightarrow\infty$ and $N\rightarrow\infty$.

###### Proof.

  1. (a)

Let $x\in\mathbb{R}$ be a continuity point of the spectral density $f$. We
need to show

$\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left|\hat{f}_{N,n}(x)-2f(x)\right|>\varepsilon\right)=0$

for all $\varepsilon>0$. Note that the continuity of the kernel function
$\kappa$ ensures the continuity of the kernel density estimator. Hence, by the
continuous mapping theorem
$\left(\hat{Z}_{1,n},\dots,\hat{Z}_{N,n}\right)\stackrel{{\scriptstyle
P}}{{\longrightarrow}}(|Z_{1}|,\dots,|Z_{N}|)$ as $n\rightarrow\infty$ implies
that $\hat{f}_{N,n}(x)\stackrel{{\scriptstyle
P}}{{\longrightarrow}}\hat{f}_{N}(x)$ as $n\rightarrow\infty$. Additionally,
the weak consistency of $\hat{f}_{N}$ for $2f$ is given at any point of
continuity of $f$.

Applying the triangle inequality yields

$\displaystyle\hskip
8.61108pt\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left|\hat{f}_{N,n}(x)-2f(x)\right|>\varepsilon\right)$
$\displaystyle\leq\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left|\hat{f}_{N,n}(x)-\hat{f}_{N}(x)\right|+\left|\hat{f}_{N}(x)-2f(x)\right|>\varepsilon\right)$
$\displaystyle\leq\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left\\{\left|\hat{f}_{N,n}(x)-\hat{f}_{N}(x)\right|>\frac{\varepsilon}{2}\right\\}\cup\left\\{\left|\hat{f}_{N}(x)-2f(x)\right|>\frac{\varepsilon}{2}\right\\}\right)$
$\displaystyle\leq\lim\limits_{N\rightarrow\infty}\underbrace{\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left|\hat{f}_{N,n}(x)-\hat{f}_{N}(x)\right|>\frac{\varepsilon}{2}\right)}_{=0}+\underbrace{\lim\limits_{N\rightarrow\infty}\mathbb{P}\left(\left|\hat{f}_{N}(x)-2f(x)\right|>\frac{\varepsilon}{2}\right)}_{=0}=0$

for all $\varepsilon>0$.

  2. (b)

For the weak uniform consistency of $\hat{f}_{N,n}$ we need to show

$\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left\|\hat{f}_{N,n}-2f\right\|_{\infty}>\varepsilon\right)=0$

for any $\varepsilon$. Similar to part (a) applying the triangle inequality
yields

$\displaystyle\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left\|\hat{f}_{N,n}-2f\right\|_{\infty}>\varepsilon\right)$
$\displaystyle\leq$
$\displaystyle\underbrace{\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left\|\hat{f}_{N,n}-\hat{f}_{N}\right\|_{\infty}>\frac{\varepsilon}{2}\right)}_{(\ast)}+\underbrace{\lim\limits_{N\rightarrow\infty}\mathbb{P}\left(\left\|f_{N}-2f\right\|_{\infty}>\frac{\varepsilon}{2}\right)}_{=0},$

where the convergence of the second summand in the above to $0$ follows from
the weak uniform consistency of $\hat{f}_{N}$. As for the first summand
$(\ast)$, following the Lipschitz continuity of the kernel $\kappa$, we can
bound

$\displaystyle\left\|\hat{f}_{N,n}-\hat{f}_{N}\right\|_{\infty}$
$\displaystyle=\sup\limits_{x\in\mathbb{R}}\left|\frac{1}{Nh}\sum\limits_{k=1}^{N}\kappa\left(\frac{x-\hat{Z}_{k,n}}{h}\right)-\frac{1}{Nh}\sum\limits_{k=1}^{N}\kappa\left(\frac{x-|Z_{k}|}{h}\right)\right|$
(17)
$\displaystyle\leq\sup\limits_{x\in\mathbb{R}}\frac{1}{Nh}\sum\limits_{k=1}^{N}\underbrace{\left|\kappa\left(\frac{x-\hat{Z}_{k,n}}{h}\right)-\kappa\left(\frac{x-|Z_{k}|}{h}\right)\right|}_{\leq
L\left|\frac{x-\hat{Z}_{k,n}}{h}-\frac{x-|Z_{k}|}{h}\right|=\frac{L}{h}\left||Z_{k}|-\hat{Z}_{k,n}\right|}$
$\displaystyle\leq\frac{L}{Nh^{2}}\sum\limits_{k=1}^{N}\left||Z_{k}|-\hat{Z}_{k,n}\right|,$
(18)

and consequently

$\displaystyle(\ast)$
$\displaystyle\leq\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\frac{L}{Nh^{2}}\sum\limits_{k=1}^{N}\left|\hat{Z}_{k,n}-|Z_{k}|\right|>\frac{\varepsilon}{2}\right)\leq\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\bigcup\limits_{k=1}^{N}\left\\{\frac{L}{Nh^{2}}\left|\hat{Z}_{k,n}-|Z_{k}|\right|>\frac{\varepsilon}{2N}\right\\}\right)$
$\displaystyle\leq\lim\limits_{N\rightarrow\infty}\sum\limits_{k=1}^{N}\underbrace{\lim\limits_{n\rightarrow\infty}\mathbb{P}\left(\left|\hat{Z}_{k,n}-|Z_{k}|\right|>\frac{\varepsilon
h^{2}}{L}\right)}_{=0}=0.$

  3. (c)

Again, the triangle inequality yields

$\displaystyle\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\left|\hat{f}_{N,n}(x)-2f(x)\right|>0\right)$
$\displaystyle\leq$
$\displaystyle\underbrace{\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\left|\hat{f}_{N,n}(x)-\hat{f}_{N}(x)\right|>0\right)}_{=0}+\underbrace{\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}\left|\hat{f}_{N}(x)-2f(x)\right|>0\right)}_{=0}=0.$

The first term on the right-hand side of the above inequality is equal to $0$
by the continuous mapping theorem and the strong consistency of
$\hat{Z}_{k,n}$ for $|Z_{k}|$, and the second term is equal to $0$ by the
strong uniform consistency of $\hat{f}_{N}$ for $2f$ imposed by the conditions
of Lemma 2 (c).

Under Lipschitz continuity of $\kappa$ we can make use of inequality (17) and
compute

$\displaystyle\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\left\|\hat{f}_{N,n}-2f\right\|_{\infty}>0\right)$
$\displaystyle\leq$
$\displaystyle\underbrace{\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}\lim\limits_{n\rightarrow\infty}\left\|\hat{f}_{N,n}-\hat{f}_{N}\right\|_{\infty}>0\right)}_{(\ast\ast)}+\underbrace{\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}\left\|\hat{f}_{N}-2f\right\|_{\infty}>0\right)}_{=0}=0,$

where the second summand is $0$ again by the strong uniform consistency of
$\hat{f}_{N}$, and for $(\ast\ast)$ note that the events
$A_{N}:=\left\\{\lim\limits_{n\rightarrow\infty}\frac{L}{Nh^{2}}\sum\limits_{k=1}^{N}\left||Z_{k}|-\hat{Z}_{k,n}\right|>0\right\\}$
are nested with $A_{N}\subseteq A_{N+1}$, hence

$\displaystyle(\ast\ast)$
$\displaystyle=\mathbb{P}\left(\lim\limits_{N\rightarrow\infty}A_{N}\right)=\lim\limits_{N\rightarrow\infty}\mathbb{P}\left(A_{N}\right)=\lim_{N\rightarrow\infty}\mathbb{P}\left(\lim\limits_{n\rightarrow\infty}\frac{L}{Nh^{2}}\sum\limits_{k=1}^{N}\left||Z_{k}|-\hat{Z}_{k,n}\right|>0\right)$
$\displaystyle\leq\lim\limits_{N\rightarrow\infty}\sum\limits_{k=1}^{N}\underbrace{\mathbb{P}\left(\lim\limits_{n\rightarrow\infty}\left||Z_{k}|-\hat{Z}_{k,n}\right|>0\right)}_{=0}=0.$

∎

##  5 Numerical results

In this section, we aim to verify our theoretical results with numerical
inference. We will consider the following four examples for the spectral
density of a SRH $S\alpha S$ process $X$.

######  Example 1.

We consider the following symmetric probability density functions on
$\mathbb{R}$ as examples.  (a)
$f_{1}(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$ , (b)
$f_{2}(x)=x^{2}e^{-\frac{|x|}{4}}$ ,  (c)
$f_{3}(x)=\frac{1}{x^{2}}\mathbbm{1}_{[1,\infty)}(|x|)$ , (d)
$f_{4}(x)=\frac{1}{2}\mathbbm{1}_{[-1,1]}(x)$ .

For the simulation of the SRH $S\alpha S$ process $X$ we make use of the
series representation in Proposition 2, i.e.

$\displaystyle
X_{N}(t)=\left(C_{\alpha}b_{\alpha}^{-1}\right)^{1/\alpha}\sum\limits_{k=1}^{N}\Gamma_{k}^{-1/\alpha}\left(G_{k}^{(1)}\cos(tZ_{k})+G_{k}^{(2)}\sin(tZ_{k})\right),$
(19)

where $\Gamma_{k}$ are the arrival times of a unit rate Poisson point process,
$G_{k}^{(i)}$, $i=1,2$, are i.i.d. $N(0,1)$ and $Z_{k}$ are i.i.d. with
probability density function $f_{i}$, $i=1,\dots,4$ from the example above.
The constants $C_{\alpha}$ and $b_{\alpha}$ are given in Proposition 2. Note
that in the above series representation of $X$ the summation is finite up to
$N\in\mathbb{N}$, where we choose $N$ large enough such that
$\Gamma_{k}^{-1/\alpha}$ are negligibly small for $k\geq N$.

###  5.1 Independent paths

In the introductory example in Section 2.2 we considered $L\in\mathbb{N}$
independent path realizations of a SRH $S\alpha S$ process $X$ with symmetric
spectral density $f$. These paths build the basis for the estimation of the
$\alpha$-sine transform of $f$. Applying the inversion method of the
$\alpha$-sine transform described in [8] allows us to reconstruct the spectral
density.

For the results in Figure 1 we simulated $L=100$ paths of the process $X$ with
index of stability $\alpha=1.5$ and spectral density $f_{1}$ of Example 1. The
paths were sampled at $n=101$ equidistant points $0=t_{1}<...<t_{n}=T$ on the
interval $[0,T]=[0,10]$. Across all paths samples of lags
$X^{(l)}(t_{i})-X^{(l)}(0)\sim S1.5S(\sigma(t_{i}))$ for $i=1,\dots,n$,
$l=1,\dots,L$ were generated. The regression-type estimators of Koutrouvelis
[9] were used to estimate $\alpha$ and the scale parameters
$\sigma(t_{i})=\left\|X(t_{i})-X(0)\right\|_{\alpha}$ of the lags. With
Equation (4) we get estimates of the $\alpha$-sine transform of $f$ at the
points $\\{t_{i}/2\\}_{i=1}^{n}$, and the inversion method in [8] yields an
estimate of the spectral density $f$.

(a) Estimate of the $\alpha$-sine transform $T_{\alpha}f_{1}$

(b) Estimate of the spectral density $f_{1}$.

Figure 1: Multiple paths inference. The solid blue line shows the estimates
and the dashed red line shows the theoretical function. The inference is based
on $L=100$ independent paths of a SRH $S1.5S$ process with spectral density
$f_{1}$ from Example 1 (a). The paths are sampled on $[0,10]$ at $n=101$
equidistant points.

One might incorporate smoothing for better estimation results, see [8, Section
6.2.1], but further analysis was omitted here, since we focused on statistics
based on a single path of SRH $S\alpha S$ processes.

###  5.2 Periodogram frequency estimation

The periodogram estimate can be jeopardized by errors caused by _aliasing_ and
_spectral leakage_ [12, Chapters 4 and 10], when both the range $T$ on which
the signal is sampled as well as number of sample points $n$ are small.
Determining all peaks in the periodogram at once proves to be a difficult
task. Distinguishing between peaks in the periodogram which actually stem from
frequencies and not spectral leakage requires a meticulous setup of tuning
parameters such as the minimum distance between peaks or their minimum height
when employing algorithms such as findpeaks in Matlab. Instead we utilize an
iterative approach to estimate the peak locations from the periodogram [13,
Chapter 3.2, p.53].

Set
$(\alpha_{k},\beta_{k})=(C_{\alpha}b_{\alpha}^{-1})^{1/\alpha}\Gamma_{k}^{-1/\alpha}(G_{k}^{(1)},G_{k}^{(2)})$
for $k=1,\dots,N$, hence Equation (19) can be written as

$\displaystyle
X_{N}(t)=\sum\limits_{k=1}^{N}\left(\alpha_{k}\cos(tZ_{k})+\beta_{k}\sin(tZ_{k})\right),$

and consider the sample $(x(1),\dots,x(n))=(X_{N}(t_{1}),\dots,X_{N}(t_{n}))$.
For the ease of notation, let $t_{k}=k$, $k=1,\dots,n$, and denote by
$I_{n}^{(0)}=I_{n}$ the periodogram computed from the sample
$x(1),\dots,x(n)$. Determine the estimate $\hat{Z}_{1,n}$ from the location of
the largest peak of $I_{n}^{(0)}$. Estimate
$(\hat{\alpha}_{1},\hat{\beta}_{1})$ by regressing
$\alpha_{1}\cos(t\hat{Z}_{1,n})+\beta_{1}\sin(t\hat{Z}_{1})$, $t=1,\dots,n$,
on $x(1),\dots,x(n)$, and compute the periodogram $I_{n}^{(1)}$ of the
residuals
$x^{(1)}(t)=x(t)-\hat{\alpha}_{1}\cos(t\hat{Z}_{1,n})-\hat{\beta}_{1}\sin(t\hat{Z}_{1})$,
$t=1,\dots,n$. Determining the location of the largest peak of $I_{n}^{(1)}$
yields the estimate $\hat{Z}_{2,n}$. Compute the regression estimates
$(\hat{\alpha}_{2},\hat{\beta}_{2})$ and the periodogram $I_{n}^{(2)}$ from
$x^{(2)}(t)=x^{(1)}(t)-\hat{\alpha}_{2}\cos(t\hat{Z}_{2,n})-\hat{\beta}_{2}\sin(t\hat{Z_{2}})$,
$t=1,\dots,n$. We repeat this process until the periodogram
$I_{n}^{(N^{\prime})}$, $N^{\prime}<N$, does not show any pronounced peaks
anymore. The resulting estimate $\hat{Z}_{1,n},\dots,\hat{Z}_{N^{\prime},n}$
will be used for the kernel density estimate of the spectral density $f$.

We use the MinPeakProminence parameter in Matlab’s findpeaks algorithm for the
iteration’s stopping criterion. The specific value of MinPeakProminence is
chosen manually depending on the desired sensitivity of findpeaks on the
iteration’s periodograms.

We consider one path realization of a SRH $S1.5S$ process for all examples
$f_{1},\dots,f_{4}$ in Figures 2 and 3. Each path is sampled on the interval
$[0,T]$, $T=500$, at $n=10^{3}$ equidistant points. The number of frequencies
in the series representation of $X$ is set to $N=10^{4}$. The histograms of
the frequency estimates $\hat{Z}_{k,n}$ from the peak locations in the
periodogram for example $f_{1},\dots,f_{4}$ are depicted in Figure 2.

The computation of the periodogram is performed with _zero-padding_ [12,
Chapter 8], i.e. zeros are added at the end of the sample resulting in a new
sample $(x(1),\dots,x(n),0,\dots,0)$ of length $\tilde{n}>n$. Hence, the
discrete Fourier transform is computed on a finer grid $\left\\{2\pi
j/\tilde{n}\right\\}_{j=1}^{\tilde{n}}$ resulting in an interpolation of the
periodogram between the actual Fourier frequencies $\left\\{2\pi
j/n\right\\}_{j=1}^{n}$, which allows us to better distinguish between peaks.
Zero-padding has no noticable effect on the computation time of the
periodogram estimate. It is practical to choose $\tilde{n}$ to be a power of 2
to make use of the (Cooley-Tukey or radix-2) fast Fourier transform and its
reduced complexity of $O(\tilde{n}\log\tilde{n})$ compared to the direct
computation of the discrete Fourier transform [12, Chapter 9]. For our
examples the samples were zero-padded with $\tilde{n}=2^{13}=8192$.

We compute kernel density estimator $\hat{f}_{N,n}$ given in Equation (16)
with the Gaussian kernel $\kappa(x)=e^{-x^{2}/2}/\sqrt{2\pi}$. As for the
bandwidth, any fixed bandwidth $h=h_{N}$ for all sample lengths $N$
immediately fulfills all the conditions for the consistency results in Lemma
2. But kernel density estimation is highly sensitive to bandwidth selection,
and a poor choice of $h$ naturally leads leads to poor performance of the
estimator. A global optimal bandwidth, which minimizes the mean integrated
squared error of the kernel density estimator, i.e. the
$\text{MISE}(\hat{f}_{N})=\int_{\mathbb{R}}\text{MSE}(\hat{f}_{N}(x))dx$,
where $\text{MSE}(\hat{f}_{N}(x))=\mathbb{E}[(\hat{f}_{N}(x)-2f(x))^{2}]$ is
given by

$\displaystyle
h^{\ast}_{N}=N^{-1/5}\left(\int_{\mathbb{R}}\kappa(t)dt\right)^{1/5}\left(\int_{\mathbb{R}}t^{2}\kappa(t)dt\right)^{-2/5}\left(\int_{\mathbb{R}}2f^{\prime\prime}(x)dx\right)^{-1/5}=O\left(N^{-1/5}\right).$
(20)

For the above optimal bandwidth the MISE and the MSE are of order
$O(N^{-5/4})$, see [20, Chapter 3.3].

The optimal bandwidth $h^{\ast}_{N}$ in Equation (20) is nice for theoretical
purposes but is not applicable in practice due to the simple fact that it
depends on the unknown spectral density $f$. Instead, many methods to
approximate the optimal bandwidth can be found in literature. _Scott’s_ or
_Silverman’s rules of thumb_ are widespread in practices for their ease of use
but assume the sample to be drawn from a Gaussian distribution [20, Section
3.4]. These methods are efficient and fast but their accuracy can only be
guaranteed in the Gaussian case or for the estimation of unimodal and close to
Gaussian densities. When the form of the sought-after density is unknown,
which is usually the case, methods like the _unbiased cross-validation_ [1,
15] or the _Sheather and Jones plug-in method_ [19] are far better suited.
There are many more cross-validation and plug-in methods at hand, see e.g.
[18] for an overview, but we applied the Sheather and Jones method as it is
already implemented in R. Similar results were achieved with other methods.

The results are given in Figure 3. Similar results were achieved using other
well-known kernel functions such as the Epanechnikov kernel or triangle
kernel.

(a) $f_{1}$

(b) $f_{2}$

(c) $f_{3}$

(d) $f_{4}$

Figure 2: Periodogram frequency estimators. Histogram of the periodogram
frequency estimators $\hat{Z}_{k}$(blue) against true absolute frequencies
$|Z_{k}|$ (orange). The estimation is based on one path of a SRH $S1.5S$
process $X$ with spectral density $f_{1},\dots,f_{4}$ sampled on $[0,500]$ at
$n=10^{3}$ equidistant points. For the simulation, the number of frequencies
in the series representation of $X$ is set to $N=10^{4}$.

(a) $f_{1}$

(b) $f_{2}$

(c) $f_{3}$

(d) $f_{4}$

Figure 3: Spectral density estimation (Gaussian kernel, Sheather-Jones
bandwidth). Kernel density estimator (blue) against spectral density (red) for
examples $f_{1},\dots,f_{4}$ with index of stability $\alpha=1.5$, sample
range $T=500$ and sample size $n=10^{3}$.

Figure 4 illustrates our results for the estimation of example function
$f_{2}$ in the cases $\alpha\in\\{1.75,1.25,0.75,0.25\\}$. As expected, the
estimation becomes more difficult as the index of stability $\alpha$
decreases. Recall that the amplitudes in the series representation of the SRH
$S\alpha S$ process $X$ in Equation (7) are given by

$R_{k}=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\Gamma_{k}^{-1/\alpha}\sqrt{\left(G_{k}^{(1)}\right)^{2}+\left(G_{k}^{(2)}\right)^{2}}.$

The smaller $\alpha$, the faster the variance mixture term
$\nu_{k}=\left(C_{\alpha}b_{\alpha}^{-1}m(\mathbb{R})\right)^{1/\alpha}\Gamma_{k}^{-1/\alpha}$
tends to $0$. This directly translates to the amplitudes $R_{k}$, which in a
sense are the effect of the frequency $Z_{k}$ on the whole process. Since the
$R_{k}$ decay much faster it becomes increasingly difficult to estimate the
corresponding frequencies $Z_{k}$ with smaller values of $\alpha$. In the
cases $\alpha=0.75$ and $\alpha=0.25$ in Figure 4 (c) and (d) the
MinPeakProminence parameter was decreased for a higher sensitivity of Matlab’s
findpeaks function.

(a) $\alpha=1.25$

(b) $\alpha=1.75$

(c) $\alpha=0.75$

(d) $\alpha=0.25$

Figure 4: Spectral density estimation (Gaussian kernel, Sheather-Jones
bandwidth). Kernel density estimator (blue) against spectral density (red) for
example $f_{2}$ with index of stability $\alpha\in\\{1.75,1.25,0.75,0.25\\}$,
sample range $T=500$ and sample size $n=10^{3}$.

Increasing the samples’ window length $T$ and therefore the sample size $n$
naturally leads to better estimation results. We increased the sampling window
to $[0,5000]$ and the sampling size to $n=10^{4}$. The results of the
estimation for all examples $f_{1},\dots,f_{4}$ and $\alpha=1.5$ are depicted
in Figure 5. Furthermore, Figure 6 shows the estimation for $\alpha=0.75$.

(a) $f_{1}$

(b) $f_{2}$

(c) $f_{3}$

(d) $f_{4}$

Figure 5: Spectral density estimation (Gaussian kernel, Sheather-Jones
bandwidth). Kernel density estimator (blue) against spectral density (red) for
examples $f_{1},\dots,f_{4}$ with index of stability $\alpha=1.5$, sample
range $T=5000$ and sample size $n=10^{4}$.

(a) $f_{1}$

(b) $f_{2}$

(c) $f_{3}$

(d) $f_{4}$

Figure 6: Spectral density estimation (Gaussian kernel, Sheather-Jones
bandwidth). Kernel density estimator (blue) against spectral density (red) for
examples $f_{1},\dots,f_{4}$ with index of stability $\alpha=0.75$, sample
range $T=5000$ and sample size $n=10^{4}$.

We clearly see the improvements for the case $\alpha=1.5$ when comparing
Figures 3 and 5. On the other hand, a smaller index of stability $\alpha$
complicates the estimation, see Figures 5 and 6.

One thing we would like to mention is the the dip of the kernel density
estimates near $0$ in examples $f_{1}$ and $f_{4}$, see Figures 3, 5 and 6 (a)
and (d), respectively. These dips are a direct consequence of a commonly known
problem in frequency estimation, i.e. very low frequencies that are close to
$0$ are difficult to detect via the periodogram. This can be resolved by an
increase of the sampling range $T$ and sample size $n$. For more details we
refer to [13, Section 3.5].

##  6 Conclusion

Harmonizable $S\alpha S$ processes are one of the three main classes of
stationary $S\alpha S$ processes. Unlike stationary moving average $S\alpha S$
processes, the harmonizable case has not received much attention from a
statistical point of view. This is mainly due to the non-ergodicity of these
processes, which inhibits the application of standard empirical methods.

We considered the special case of stationary real harmonizable $S\alpha S$, in
which the circular control measure of the process is the product of the
uniform probability measure on the unit circle $(0,2\pi)$ and the so-called
control measure $m$. Assuming that the control measure $m$ has a symmetric
density $f$, the goal our work was to develop a consistent and efficient
statistical procedure to estimate $f$.

The series representation in Proposition 2 shows that a SRH $S\alpha S$
process $X$ is a conditional non-ergodic harmonic Gaussian process. In Theorem
2 we derived the non-ergodic limits of the empirical characteristic function
of $X$ and the lag process $\left\\{X(t+h)-X(t):t\in\mathbb{R}\right\\}$,
$h>0$. These limits can be explicitly given in terms of the Bessel function of
the first kind of order $0$ and the processes’ invariant sets.

Additionally, Theorem 2 also yields an equivalent series representation of $X$
in terms of amplitudes $\left\\{R_{k}\right\\}_{k=1}^{\infty}$, i.i.d. uniform
phases $\left\\{U_{k}\right\\}_{k=1}^{\infty}$ and frequencies
$\left\\{Z_{k}\right\\}_{k=1}^{\infty}$, see Equation (7). The frequencies
$Z_{k}$ are i.i.d. distributed with probability density function $f$. And
although these quantities are random, they are predetermined for a given path.

The periodogram is a fast and efficient tool for the estimation of the
absolute frequencies $|Z_{k}|$ as it relies on the fast Fourier algorithm. We
show that the frequency estimators are strongly consistent in Theorem 3.
Applying kernel density estimation yields an estimate of the spectral density.
Under minimal assumptions on the kernel function and the bandwidth, various
weak, strong as well as pointwise and uniform consistency results for our
spectral density estimator are proven in Theorem 4.

An extensive numerical analysis shows that our proposed estimation method
achieves satisfying results well on a variety of examples with index of
stability $\alpha\in(0,2)$. Results based on samples of size $n=1000$ alone
were demonstrated in Figures 3 and 4. The smaller the index of stability
$\alpha$ is the more difficult the estimation of the spectral density becomes
as peaks in the periodogram are harder to detect. But in general, good results
can still be achieved by increasing the sensitivity of the peak detecting
function. Moreover, an increase of the sampling range $T$ and sample size $n$
results in significant improvements of the estimation as shown in Figures 5
and 6. Further improvements might be achieved with the application of
different window functions as we have mentioned before.

Ultimately, aside from the consistency of our spectral density estimator and
its ease of computation, we would like to highlight that no other requirements
or prior knowledge on the process (such as e.g. the index of stability
$\alpha$) is needed for the estimation of the spectral density $f$.

##  Appendix A Periodogram computation

Recall the discrete Fourier transform of the sample $x(1),\dots,x(n)$ for
$\theta>0$ in Equation (12), i.e.

$\displaystyle F_{n}(\theta\delta)$
$\displaystyle=n^{-1}\sum\limits_{l=1}^{\infty}\Bigg{(}\frac{e^{iU_{l}}R_{l}}{2}\underbrace{\sum\limits_{j=1}^{n}e^{ij\delta\left(Z_{l}-\theta\right)}}_{=S_{l}^{(1)}(\theta,\delta)}+\frac{e^{-iU_{l}}R_{l}}{2}\underbrace{\sum\limits_{j=1}^{n}e^{-ij\delta\left(Z_{l}+\theta\right)}}_{=S_{l}^{(2)}(\theta,\delta)}\Bigg{)},$

where the sums $S_{l}^{(1)}$ and $S_{l}^{(2)}$ can be explicitly expressed by

$\displaystyle
S_{l}^{(1)}(\theta,\delta)=e^{i\delta\left(Z_{l}-\theta\right)}\frac{e^{in\delta\left(Z_{l}-\theta\right)}-1}{e^{i\delta\left(Z_{l}-\theta\right)}-1}\quad\text{and}\quad
S_{l}^{(2)}(\theta,\delta)=e^{-in\delta\left(Z_{l}-\theta\right)}\frac{e^{in\delta\left(Z_{l}+\theta\right)}-1}{e^{i\delta\left(Z_{l}+\theta\right)}-1}.$

The periodogram $I_{n}\left(\theta\right)$ is defined as the squared absolute
value of the discrete Fourier transform. Note that
$\left|x+y\right|^{2}=|x|^{2}+|y|^{2}+x^{\ast}y+xy^{\ast}=|x|^{2}+|y|^{2}+2Re(x^{\ast}y)$,
where $x^{\ast}$ denotes the complex conjugate of $x$. This yields

$\displaystyle I\left(\theta\right)$
$\displaystyle=\Bigg{|}n^{-1}\sum\limits_{l=1}^{\infty}\Bigg{(}\frac{e^{iU_{l}}R_{l}}{2}S_{l}^{(1)}(\theta,\delta)+\frac{e^{-iU_{l}}R_{l}}{2}S_{l}^{(2)}(\theta,\delta)\Bigg{)}\Bigg{|}^{2}$
$\displaystyle=n^{-2}\sum\limits_{l=1}^{\infty}\frac{R_{l}^{2}}{4}\left[\left|S_{l}^{(1)}(\theta,\delta)\right|^{2}+\left|S_{l}^{(2)}(\theta,\delta)\right|^{2}+2Re\left(\left(e^{iU_{l}}S_{l}^{(1)}(\theta,\delta)\right)^{\ast}e^{-iU_{l}}S_{l}^{(2)}(\theta,\delta)\right)\right]$
$\displaystyle+n^{-2}\sum\limits_{l=1}^{\infty}\sum\limits_{m>l}\frac{R_{l}R_{m}}{2}Re\left(\left(e^{iU_{l}}S_{l}^{(1)}(\theta,\delta)+e^{-iU_{l}}S_{l}^{(2)}(\theta,\delta)\right)^{\ast}\left(e^{iU_{m}}S_{m}^{(1)}(\theta,\delta)+e^{-iU_{m}}S_{m}^{(2)}(\theta,\delta)\right)\right)$
$\displaystyle=(\ast_{1})+(\ast_{2}).$

For ease of notation we assume $\delta=1$ and define the function
$S_{n}(x)=\sin\left(\frac{n}{2}x\right)/\sin\left(\frac{1}{2}x\right)$. In the
following we will not further denote the functions dependence on
$n\in\mathbb{N}$ with a subscript and just write $I=I_{n}$ and $S=S_{n}$
instead. The function $S$ is well-defined for all $x\in\mathbb{R}$, in
particular it holds that $S(0)=n$, hence it follows that
$\lim_{n\rightarrow\infty}n^{-1}S(x)=\pm 1$ for any $x=j\pi$,
$j\in\mathbb{Z}$, and $0$ otherwise.

Note that
$\left|e^{ix}-1\right|^{2}=\left(\cos(x)-1\right)^{2}+\sin^{2}(x)=2-2\cos(x)=4\sin^{2}\left(\frac{x}{2}\right)$
for all $x\in\mathbb{R}$, hence we can compute

$\displaystyle\left|S_{l}^{(1)}(\theta,1)\right|^{2}=\frac{\sin^{2}\left(\frac{n}{2}(Z_{l}-\theta)\right)}{\sin^{2}\left(\frac{1}{2}(Z_{l}-\theta)\right)}=S^{2}\left(Z_{l}-\theta\right),\quad\left|S_{l}^{(2)}(\theta,1)\right|^{2}=\frac{\sin^{2}\left(\frac{n}{2}(Z_{l}+\theta)\right)}{\sin^{2}\left(\frac{1}{2}(Z_{l}+\theta)\right)}=S^{2}\left(Z_{l}+\theta\right)$

for all $l\in\mathbb{N}$. Furthermore, it holds that

$\displaystyle\left(e^{iU_{l}}S_{l}^{(1)}(\theta,1)\right)^{\ast}\left(e^{-iU_{l}}S_{l}^{(2)}(\theta,1)\right)=e^{-i2U_{l}}\left(S_{l}^{(1)}\left(\theta,1\right)\right)^{\ast}S_{l}^{(2)}\left(\theta,1\right)$

with

$\displaystyle\left(S_{l}^{(1)}\left(\theta,1\right)\right)^{\ast}S_{l}^{(2)}\left(\theta,1\right)=S(Z_{l}-\theta)S(Z_{l}+\theta)\cos\left((n+1)Z_{l}\right)-iS(Z_{l}-\theta)S(Z_{l}+\theta)\sin\left((n+1)Z_{l}\right).$

Hence,

$\displaystyle
Re\left(e^{-i2U_{l}}\left(S_{l}^{(1)}\left(\theta,1\right)\right)^{\ast}S_{l}^{(2)}\left(\theta,1\right)\right)$
$\displaystyle=$ $\displaystyle
Re\left(e^{-i2U_{l}}\right)Re\left(\left(S_{l}^{(1)}\left(\theta,1\right)\right)^{\ast}S_{l}^{(2)}\left(\theta,1\right)\right)-Im\left(e^{-i2U_{l}}\right)Im\left(\left(S_{l}^{(1)}\left(\theta,1\right)\right)^{\ast}S_{l}^{(2)}\left(\theta,1\right)\right)$
$\displaystyle=$ $\displaystyle
S(Z_{l}-\theta)S(Z_{l}+\theta)\left(\cos(2U_{l})\cos\left((n+1)Z_{l}\right)-\sin(2U_{l})\sin\left((n+1)Z_{l}\right)\right)$
$\displaystyle=$ $\displaystyle
S(Z_{l}-\theta)S(Z_{l}+\theta)\cos\left(2U_{l}+(n+1)Z_{l}\right).$ (21)

Therefore,

$\displaystyle(\ast_{1})=n^{-2}\sum\limits_{l=1}^{\infty}\frac{R_{l}^{2}}{4}\left[S^{2}(Z_{l}-\theta)+S^{2}(Z_{l}+\theta)+2\cos\left(2U_{l}+(n+1)Z_{l}\right)S(Z_{l}-\theta)S(Z_{l}+\theta)\right].$

Multiplying out the expression in the real part of the sum $(\ast_{2})$ and
computing the real part of each summand similar to Equation (A) yields

$\displaystyle(\ast_{2})$
$\displaystyle=n^{-2}\sum\limits_{l=1}^{\infty}\sum\limits_{m>l}\frac{R_{l}R_{m}}{2}$
$\displaystyle\times\left[\left(S(Z_{l}-\theta)S(Z_{m}-\theta)+S(Z_{l}+\theta)S(Z_{m}+\theta)\right)\cos\left(U_{l}-U_{m}+\frac{1}{2}(n+1)(Z_{l}-Z_{m})\right)\right.$
$\displaystyle\left.+\left(S(Z_{l}-\theta)S(Z_{m}+\theta)+S(Z_{l}+\theta)S(Z_{m}-\theta)\right)\cos\left(U_{l}+U_{m}+\frac{1}{2}(n+1)(Z_{l}+Z_{m})\right)\right].$

Note that one can easily verify that the frequencies $Z_{l}$ and $Z_{m}$
inside the function $S$ can be replaced by their absolute values $|Z_{l}|$ and
$|Z_{m}|$. Ultimately, we get

$\displaystyle I(\theta)$
$\displaystyle=n^{-2}\sum\limits_{l=1}^{\infty}\frac{R_{l}^{2}}{4}\left[S^{2}(|Z_{l}|-\theta)+S^{2}(|Z_{l}|+\theta)+2\cos\left(2U_{l}+(n+1)Z_{l}\right)S(|Z_{l}|-\theta)S(|Z_{l}|+\theta)\right]$
$\displaystyle+n^{-2}\sum\limits_{l=1}^{\infty}\sum\limits_{m>l}\frac{R_{l}R_{m}}{2}$
$\displaystyle\times\left[\left(S(|Z_{l}|-\theta)S(|Z_{m}|-\theta)+S(|Z_{l}|+\theta)S(|Z_{m}|+\theta)\right)\cos\left(U_{l}-U_{m}+\frac{1}{2}(n+1)(Z_{l}-Z_{m})\right)\right.$
$\displaystyle+\left.\left(S(|Z_{l}|-\theta)S(|Z_{m}|+\theta)+S(|Z_{l}|+\theta)S(|Z_{m}|-\theta)\right)\cos\left(U_{l}+U_{m}+\frac{1}{2}(n+1)(Z_{l}+Z_{m})\right)\right].$
(22)

The periodogram’s behavior described in Equation (14), i.e. its convergence to
$0$ with rate $O(n^{-2})$ as $n\rightarrow\infty$ except for when
$\theta\rightarrow\infty|Z_{k}|$ for some $k\in\mathbb{N}$ can be readily
established by Lebesgue’s dominated convergence theorem and the fact that
$n^{-1}|S(x)|$ is bounded by $1$.

## References

  * Bowman [1984] Bowman, A.W..  An alternative method of cross-validation for the smoothing of kernel density estimates.  Biometrika 71 (2), pp. 353-360, 1984\.
  * Brockwell and Davis [2016] Brockwell, P.J., Davis, R.A..  Time series: Theory and methods.  Springer, 2016.
  * Cambanis et al. [1987] Cambanis, S., Hardin Jr., C.D., Weron, A..  Ergodic properties of stationary stable processes.  Stochastic Processes and their Applications 24 (1), pp. 1-18, 1987.
  * Doob [1991] Doob, J.L..  Stochastic processes.  Wiley, 1991.
  * Dym and McKean [1976] Dym, H., McKean, H.P..  Gaussian processes, function theory, and the inverse spectral problem.  Academic Press, 1976.
  * Fuller [1996] Fuller, W.A..  Introduction to statistical time series.  Wiley, 1996.
  * Hannan [1973] Hannan, E.J..  The estimation of frequency.  Journal of Applied Probability 10 (3), pp. 510-519, 1973.
  * Hoang and Spodarev [2021] Hoang, L.V., Spodarev, E..  Inversion of $\alpha$-sine and $\alpha$-cosine transforms on $\mathbb{R}$.  Inverse Problems 37 (8), 2021.
  * Koutrouvelis [1980] Koutrouvelis, I.A..  Regression-type estimation of the parameters of stable laws.  Journal of the American Statistical Association 75 (372), pp. 918-928, 1980.
  * Lindgren [2012] Lindgren, G..  Stationary stochastic processes: Theory and applications.  Chapmann and Hall, 2012.
  * McCulloch [1986] McCulloch, J.H..  Simple consistent estimators of stable distribution parameters.  Communications in Statistics - Simulation and Computation 15, pp. 1109-1136, 1986.
  * Oppenheim et al. [1999] Oppenheim, A.V., Buck, J.R., Schafer, R.W..  Discrete-time signal processing.  Prentice Hall, 1999.
  * Quinn and Hannan [2001] Quinn, B.G., Hannan, E.J..  The estimation and tracking of frequency.  Cambridge University Press, 2001\.
  * Rosinski [1995] Rosinski, J..  On the structure of stationary stable processes.  The Annals of Probability 23 (3), pp. 1163-1187, 1995\.
  * Rudeom [1982] Rudeom, M..  Empirical choice of histograms and kernel density estimators.  Scandinavian Journal of Statistics 9 (2), pp. 65-78, 1982.
  * Samorodnitsky [2005] Samorodnitsky, G..  Null flows, positive flows and the structure of stationary symmetric stable srocesses.  The Annals of Probability 33 (5), pp. 1782-1803, 2005\.
  * Samorodnitsky and Taqqu [1994] Samorodnitsky, G., Taqqu, M.S..  Stable non-Gaussian random processes: Stochastic models with infinite variance.  Chapman and Hall, 1994.
  * Sheather [2004] Sheather, S.J..  Density estimation.  Statistical Science 19 (4), pp. 588–597, 2004\.
  * Sheather and Jones [1984] Sheather, S.J., Jones, M.C..  A reliable data-based bandwidth selection method for kernel density estimation.  Journal of the Royal Statistical Society, Series B (Methodological) 53 (3), pp. 683-690, 1984.
  * Silverman [1986] Silverman, B.W..  Density esetimation for statistics and data analysis.  Chapman and Hall, 1986.
  * Ślezak [2017] Ślezak, J..  Asymptotic behaviour of time averages for non-ergodic Gaussian processes.  Annals of Physics 383, pp. 285-311, 2017\.
  * Wied and Weißbach [2012] Wied, D., Weißbach, R..  Consistency of the kernel density estimator - a survey.  Statistical Papers 53 (1), pp. 1-21, 2012\.