#  Periodic trawl processes: Simulation, statistical inference and
applications in energy markets

Almut E. D. Veraart
Department of Mathematics, Imperial College London
180 Queen’s Gate, London, SW7 2AZ, UK
<EMAIL_ADDRESS>
###### Abstract

This article introduces the class of periodic trawl processes, which are
continuous-time, infinitely divisible, stationary stochastic processes, that
allow for periodicity and flexible forms of their serial correlation,
including both short- and long-memory settings. We derive some of the key
probabilistic properties of periodic trawl processes and present relevant
examples. Moreover, we show how such processes can be simulated and establish
the asymptotic theory for their sample mean and sample autocovariances.
Consequently, we prove the asymptotic normality of a (generalised) method-of-
moments estimator for the model parameters. We illustrate the new model and
estimation methodology in an application to electricity prices.

MSC2020 subject classifications: Primary 62M15, 62F99, 60G10; secondary 62P99,
60F05, 60G57, 60E07.
Keywords: Periodic trawl process, Lévy basis, infinite divisibility,
simulation, asymptotic theory, weak dependence, generalised method of moments.

##  1 Introduction

Time series with periodic behaviours are ubiquitous and are present in many
applications such as climate science, meteorology, oceanography, hydrology,
economics, and communications networks. Our article is particularly motivated
by applications in energy, carbon, and temperature markets. Here, we need to
consider different types of seasonalities and periodicities when, for example,
modelling (renweable) energy production, demand and prices that are influenced
by the seasonal weather patterns, as well as different daily and weekly
fluctuations and economic cycles.

Traditionally, stochastic models for electricity prices focussed on mean-
reverting processes since spot prices can be regarded as equilibrium prices
balancing supply and demand. Barndorff-Nielsen et al. (2013) argued that, in
order to achieve a weaker version of mean reversion, one could model prices
directly by a stationary process. We follow the same idea in this article and
are interested in developing a flexible class of stationary stochastic
processes which can allow for periodic behaviour rather than employing widely-
used seasonal-trend decomposition models.

There is a broad literature on cyclostationary processes (and their
extensions), which are typically constructed by combining stationary processes
with periodic functions. The textbook by Hurd & Miamee (2007) contains
detailed historical background and Gardner et al. (2006) provides an excellent
review; the name _cyclostationary_ processes first appears in Bennett (1958),
whereas other authors referred to them as _periodically correlated_ processes,
see e.g. Gladyšev (1961), Hurd (1969), and Genton & Hall (2007), Das & Genton
(2021) for more recent work. There has been much interest in developing
suitable statistical methodologies for periodic and cyclostationary processes
and their extensions. For instance, Quinn & Thomson (1991) focuses on
estimating the frequency of a periodic function, whereas Hall et al. (2000)
presents a nonparametric estimator of a periodic function.

Moreover, from a stochastic modelling perspective, it is relevant to construct
models which can combine periodicity with either short- or long-memory
settings, with the latter being typically more challenging. However, there is
a long tradition in time series analysis tackling long-memory settings and
incorporating periodicities: Fractional differencing, introduced by Hosking
(1981) and further developed by Andel (1986) lead to long memory time series
models and fractionally integrated ARMA (FARIMA) processes, in particular. The
processes were further generalised by Gray et al. (1989), who introduced
Gegenbauer and Gegenbauer ARMA (GARMA) processes, see also Woodward et al.
(1998) for a $k$-factor GARMA long-memory model and Espejo et al. (2015) for
extensions to Gegenbauer random fields. The concept of cyclical long-memory
processes and the appropriate statistical methodology has been further
explored in Arteche & Robinson (2000), Giraitis et al. (2001), Ferrara &
Guégan (2001) and Arteche (2002). For example, Hidalgo & Soulier (2004) study
the estimation of the location and exponent of the spectral singularity of a
long-memory process. They maximise the periodogram to estimate the location of
the singularity, in the spirit of Yajima (2007), who proved consistency in the
Gaussian case and proved consistency under weaker (non-Gaussian) assumptions
Hidalgo & Soulier (2004). They used a GPH estimator, see Geweke & Porter-Hudak
(1983), for the long-memory parameter. Also, Hidalgo (2005) discusses the
semiparametric estimation of stationary processes whose spectra have an
unknown pole. Maddanu (2022) presents a harmonically weighted filter for
cyclical long-memory processes. In practice, one often estimates the
parameters of a cyclical/seasonal long-memory process in two steps: First, one
estimates the location of the singularity in the spectral density function,
second, with the location parameter now fixed, one estimates the long memory
parameter, see Ferrara & Guégan (2001) and Espejo et al. (2015). More
recently, a joint estimation of the periodicity and the long-memory parameter
have been developed in Alomari et al. (2020); Ayache et al. (2022).

While there is already extensive literature on discrete-time periodic
processes with a short or long memory, much less is known in the continuous-
time setting. To fill this gap, this article introduces the class of _periodic
trawl processes_ , which are continuous-time, infinitely divisible, stationary
stochastic processes, which allow for periodicity and flexible forms of their
serial correlation, including both short- and long-memory settings.

Periodic trawl processes extend the class of trawl processes first introduced
under this name in Barndorff-Nielsen (2011) but were earlier mentioned (under
a different name) in Wolpert & Taqqu (2005) and Wolpert & Brown (2011).
Motivated by applications in high-frequency financial data, Barndorff-Nielsen
et al. (2014) and Veraart (2019) introduced integer-valued trawl (IVT)
processes in a univariate and multivariate framework, respectively. The
statistical methodology for such processes was further developed in Bennedsen
et al. (2023); Shephard & Yang (2017, 2016). Moreover, it has been shown that
trawl processes can be applied in hierarchical models for extremes, see Noven
et al. (2018); Courgeau & Veraart (2022); Bacro et al. (2020). There is also a
growing interest in discrete-time versions of trawl processes, which were
first introduced by Doukhan et al. (2019), see also Doukhan et al. (2020), and
limit theorems for discrete- and continuous-time trawl processes, see for
instance Grahovac et al. (2018); Paulauskas (2020); Talarczyk & Treszczotko
(2020); Pakkanen et al. (2021).

###  1.1 Outline and main contributions of the article

The outline for the remainder of the article is as follows. Section 2 briefly
reviews mixed moving average processes (MMAPs) and then defines (periodic)
trawls as special cases of MMAPs. It presents some of the key probabilistic
properties of periodic trawl processes and relevant examples. Next, Section 3
describes how (periodic) trawl processes can be simulated and the
corresponding R code is available via the R package ambit on CRAN, see Veraart
(2022). Next, we tackle the question of how the parameters of periodic trawl
processes can be estimated. In order to develop a suitable inference
methodology for periodic trawl processes, we first develop a general
asymptotic theory for MMAPs in Section 4. Here we present central limit
theorems for the sample mean, sample autocovariance and sample
autocorrelations of MMAPs. Next, Section 5 tailors the asymptotic theory to
the case of (periodic) trawl processes and develops suitable moment-based
estimators. Our main focus here is on estimating the parameters of the kernel
function in the periodic trawl process, which determine the decay of the
autocorrelation function as well as the periodicity. This theory will be
developed under the assumption that the periodicity of the process is known.
Next, the more general setting is presented in Section 6, where the asymptotic
theory is developed for a generalised method-of-moment approach. The new
estimation methodology is illustrated in an empirical study of electricity
prices in Section 7; the code used in the empirical study is available on
GitHub and archived in Zenodo, see Veraart (2023). Section 8 concludes. The
proofs of the theoretical results, additional examples and detailed
discussions of the technical assumptions are relegated to the Appendix, see
Section A.

##  2 Mixed moving average and periodic trawl processes

Let us start off by reviewing the definition of the broader class of mixed
moving average processes (MMAPs), see for instance Surgailis et al. (1993);
Fasen (2005); Fuchs & Stelzer (2013); we will then introduce periodic trawl
processes as a special case of MMAPs.

###  2.1 Background

We briefly review the concept of Lévy bases. To this end, consider a Borel set
$S\subseteq\mathbb{R}$ and denote by $\mathcal{S}=\mathcal{B}(S)$ the Borel
$\sigma$-algebra on $S$ and by $\mathcal{B}_{b}(S)$ the bounded Borel sets of
$\mathcal{S}$.

######  Definition 1.

Let $L=\\{L(A)\,:\,A\in\mathcal{B}_{b}(S)\\}$ denote a collection of
$\mathbb{R}$-valued random variables. Then, $L$ is called an
$\mathbb{R}$-valued _Lévy basis_ on $S$ if the following three conditions
hold: 1) Random measure: For any sequence $A_{1},A_{2},\dots$ of disjoint
elements of $\mathcal{B}_{b}(S)$ satisfying
$\cup_{j=1}^{\infty}A_{j}\in\mathcal{B}_{b}(S)$ it holds that
$L\left(\cup_{j=1}^{\infty}A_{j}\right)=\sum_{j=1}^{\infty}L(A_{j})$ a.s.. 2)
Independent scatteredness: For any sequence $A_{1},A_{2},\dots$ of disjoint
elements of $\mathcal{B}_{b}(S)$, the random variables
$L(A_{1}),L(A_{2}),\dots$ are independent. 3) Infinite divisibility: For any
$A\in\mathcal{B}_{b}(S)$, the law of $L(A)$ is infinitely divisible (ID).

Recall that a Lévy basis admits a Lévy-Khintchine representation, see Rajput &
Rosiński (1989, Proposition 2.1 (a)).

######  Proposition 1.

Let $L$ denote a Lévy basis, and let $\theta\in\mathbb{R}$ and
$A\in\mathcal{B}_{b}(S)$, then the cumulant function of $L$ is given by
$C(\theta;L(A))=\mathrm{Log}\left(\mathbb{E}(\exp(i\theta L(A))\right)$, where

$\displaystyle\begin{split}C(\theta;L(A))&=i\theta\zeta^{*}(A)-\frac{1}{2}\theta^{2}a^{*}(A)+\int_{\mathbb{R}}\left(e^{i\theta\xi}-1-i\theta\xi\mathbb{I}_{[-1,1]}(\xi)\right)n(d\xi,A),\end{split}$
(1)

where $\zeta^{*}$ is a signed measure on $\mathcal{B}_{b}(S)$, $a^{*}$ is a
measure on $\mathcal{B}_{b}(S)$, and $n(\cdot,\cdot)$ is the generalised Lévy
measure, i.e. $n(d\xi,A)$ is a Lévy measure on $\mathbb{R}$ for fixed
$A\in\mathcal{B}_{b}(S)$ and a measure on $\mathcal{B}_{b}(S)$ for fixed $dx$.

Here $\mathrm{Log}$ denotes the _distinguished logarithm_ , see Sato (1999, p.
33).

In the following, we shall restrict our attention to homogeneous Lévy bases.

######  Definition 2 (Homogeneous Lévy basis).

A Lévy basis $L$ with Lévy-Khintchine representation (1) is called homogeneous
if we have the following representations for ${\bf z}\in\mathbb{R}^{2}$:

$\displaystyle\zeta^{*}(d{\bf z})=\zeta d{\bf z},$ $\displaystyle a^{*}(d{\bf
z})=ad{\bf z},$ $\displaystyle n(d\xi,d{\bf z})=\ell(d\xi)d{\bf z},$ (2)

for constants $\zeta\in\mathbb{R},a\geq 0$ and a Lévy measure $\ell$. We call
$(\zeta,a,\ell)$ the characteristic triplet associated with $L$.

In the case of a homogeneous Lévy basis with characteristic triplet
$(\zeta,a,\ell)$ the Lévy-Khintchine representation simplifies to

$\displaystyle\begin{split}C(\theta;L(A))&=\mathrm{Leb}(A)\left[i\theta\zeta-\frac{1}{2}\theta^{2}a+\int_{\mathbb{R}}\left(e^{i\theta\xi}-1-i\theta\xi\mathbb{I}_{[-1,1]}(\xi)\right)\ell(d\xi)\right],\end{split}$
(3)

for $\theta\in\mathbb{R}$ and $A\in\mathcal{B}_{b}(S)$. In the following, we
shall denote by $L^{\prime}$ the _Lévy seed_ associated with the Lévy basis
with characteristic triplet $(\zeta,a,\ell)$, i.e. a random variable with
cumulant function given by

$\displaystyle C(\theta;L^{\prime})=\mathrm{Log}\left(\mathbb{E}(\exp(i\theta
L^{\prime}))\right)=i\theta\zeta-\frac{1}{2}\theta^{2}a+\int_{\mathbb{R}}\left(e^{i\theta\xi}-1-i\theta\xi\mathbb{I}_{[-1,1]}(\xi)\right)\ell(d\xi),$
(4)

for $\theta\in\mathbb{R}$. We recall that a homogeneous Lévy basis admits a
Lévy-Itô representation, cf. Pedersen (2003, Proposition 4.5):

######  Proposition 2.

For a homogeneous Lévy basis $L$ with a characteristic triplet
$(\zeta,a,\ell)$ there exists a modification $L^{*}$ with the same
characteristic triplet that has the following Lévy-Itô decomposition. Let
$A\in\mathcal{B}_{b}(S)$, then

$\displaystyle L^{*}(A)$
$\displaystyle=\zeta\mathrm{Leb}(A)+W(A)+\int_{\\{|y|\leq
1\\}}y(N-n)(dy,A)+\int_{\\{|y|>1\\}}yN(dy,A),$ (5)

for a Gaussian basis $W$, with characteristic triplet $(0,a,0)$, i.e.
$W(A)\sim N(0,a\mathrm{Leb}(A))$, and a Poisson basis $N$ (independent of $W$)
with compensator $n(dy;A)=\mathbb{E}(N(dy;A))$, where $n(dx,d{\bf
z})=\ell(dx)d{\bf z}$.

###  2.2 Definition of a mixed moving average process

######  Definition 3 (Mixed moving average process).

Let $L$ denote a homogeneous Lévy basis with characteristic triplet
$(\zeta,a,\ell)$, $f:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ be a
deterministic function. The stochastic process $Y=(Y_{t})_{t\geq 0}$ with

$\displaystyle Y_{t}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}f(x,t-s)L(dx,ds),$ (6)

is called a mixed moving average process (MMAP).

An application of Rajput & Rosiński (1989, Theorem 2.7) leads to the following
result.

######  Proposition 3.

The MMA process defined in (8) is well-defined, if and only if the following
integrability conditions hold almost surely:

$\int_{\mathbb{R}\times\mathbb{R}}|V_{1}(f(x,t-s))|dxds<\infty,\qquad
a\int_{\mathbb{R}\times\mathbb{R}}|f(x,t-s)|^{2}dxds<\infty,\\\
\int_{\mathbb{R}\times\mathbb{R}}V_{2}(f(x,t-s))dxds<\infty,$ (7)

where for $\varrho(\xi):=\xi\mathbb{I}_{[-1,1]}(\xi)$, we define
$V_{1}(u):=u\zeta+\int_{\mathbb{R}}\left(\varrho(\xi
u)-u\varrho(\xi)\right)\ell(d\xi)$ and $V_{2}(u):=\int_{\mathbb{R}}\min(1,|\xi
u|^{2})\ell(d\xi)$. In that case, the cumulant function of the MMA process is
given by

$\displaystyle C(\theta;Y_{t})$
$\displaystyle=\mathrm{Log}(\mathbb{E}(\exp(i\theta
Y_{t})))=\int_{\mathbb{R}\times\mathbb{R}}C(\theta f(x,t-s);L^{\prime})dxds$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}C(\theta
f(x,u);L^{\prime})dxdu$
$\displaystyle=i\theta\zeta\int_{\mathbb{R}\times\mathbb{R}}f(x,u)dxdu-\frac{1}{2}\theta^{2}a\int_{\mathbb{R}\times\mathbb{R}}f^{2}(x,u)dxdu$
$\displaystyle+\int_{\mathbb{R}\times\mathbb{R}}\int_{\mathbb{R}}\left(e^{i\theta\xi
f(x,u)}-1-i\theta\xi f(x,u)\mathbb{I}_{[-1,1]}(\xi
f(x,u))\right)\ell(d\xi)dxdu,$

for $\theta\in\mathbb{R}$, which implies that the process is stationary and
infinitely divisible.

###  2.3 Periodic trawl processes

We will now show that trawl processes and the new class of periodic trawl
processes fall into the MMAP framework. First, recall the definition of a
periodic function.

######  Definition 4.

A function $p:[0,\infty)\to\mathbb{R}$ is called periodic with period $\tau>0$
if $p(x+\tau)=p(x)$ for all $x\geq 0$.

Now, we have all the ingredients for defining a periodic trawl process.

######  Definition 5 (Periodic trawl process).

Let $L$ denote a homogeneous Lévy basis with characteristic triplet
$(\zeta,a,\ell)$, $p:[0,\infty)\to\mathbb{R}$ a periodic function with period
$\tau>0$ and $g:[0,\infty)\to\mathbb{R}$ a continuous, monotonically
decreasing function. The stochastic process $Y=(Y_{t})_{t\geq 0}$ with

$\displaystyle Y_{t}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}f(x,t-s)L(dx,ds),$ (8)

with $f(x,t-s)=p(t-s)\mathbb{I}_{(0,g(t-s))}(x)\mathbb{I}_{[0,\infty)}(t-s)$
is called a $\tau$-periodic trawl process.

######  Remark 1.

In the case where $p\equiv 1$, we obtain a monotonic trawl process without
periodicity. Note that we use a slightly different notation than that used in
Barndorff-Nielsen et al. (2014). Here $g$ is defined for non-negative
arguments.

######  Remark 2.

Barndorff-Nielsen et al. (2014) proposed adding a periodic function as a
multiplicative factor to $g$ rather than as a kernel function as in (8).
However, we find that the version proposed above appears to be more
analytically tractable, as we shall discuss in more detail in the appendix;
see Remark 11. Note also that if one is interested in integer-valued trawl
processes, then they can be obtained by choosing $p$ to be integer-valued.

The following result is a special case of Proposition 3.

######  Proposition 4.

The $\tau$-periodic trawl process defined in Definition 5 is well defined, in
the sense that the stochastic integral exists, if and only if the
integrability conditions stated in (7) hold almost surely. In that case, the
cumulant function of the periodic trawl process is given by

$\displaystyle C(\theta;Y_{t})$
$\displaystyle=\mathrm{Log}(\mathbb{E}(\exp(i\theta
Y_{t})))=\int_{-\infty}^{t}\int_{\mathbb{R}}C(\theta
p(t-s)\mathbb{I}_{(0,g(t-s))}(x);L^{\prime})dxds$
$\displaystyle=\int_{0}^{\infty}\int_{\mathbb{R}}C(\theta
p(u)\mathbb{I}_{(0,g(u))}(x);L^{\prime})dxdu$
$\displaystyle=i\theta\zeta\int_{0}^{\infty}g(u)p(u)du-\frac{1}{2}\theta^{2}a\int_{0}^{\infty}g(u)p^{2}(u)du$
$\displaystyle+\int_{0}^{\infty}g(u)\int_{\mathbb{R}}\left(e^{i\theta\xi
p(u)}-1-i\theta\xi p(u)\mathbb{I}_{[-1,1]}(\xi p(u))\right)\ell(d\xi)du,$

for $\theta\in\mathbb{R}$, which implies that the process is stationary and
infinitely divisible.

######  Remark 3.

A sufficient condition for the integrability conditions to hold is that $p$ is
continuous and $\int_{0}^{\infty}g(u)du<\infty$. This is indeed the scenario
we typically consider for applications.

####  2.3.1 Second-order properties of periodic trawl processes

Let us now study the second-order properties of periodic trawl processes. The
proofs of our results are given in Appendix A.1.

######  Proposition 5.

Let $Y$ denote a $\tau$-periodic trawl process as defined in Definition 5. For
$t\geq 0$, we have
$\mathbb{E}(Y_{t})=\mathbb{E}(L^{\prime})\int_{0}^{\infty}p(u)g(u)du$,$\mathrm{Var}(Y_{t})=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}p^{2}(u)g(u)du$,
$\mathrm{Cov}(Y_{0},Y_{t})=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}p(u)p(t+u)g(t+u)du$
and

$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})=\frac{\int_{0}^{\infty}p(u)p(t+u)g(t+u)du}{\int_{0}^{\infty}p^{2}(u)g(u)du}.$

We get the following important result.

######  Proposition 6.

Let $Y$ denote a $\tau$-periodic trawl process as defined in Definition 5.
Suppose that the periodic function $p$ is continuous. Then there exists a
continuous $\tau$-periodic function $c:[0,\infty)\to\mathbb{R}$, which is
proportional to $p$ and satisfies $c(0)=1$, such that

$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})$
$\displaystyle=c(t)\frac{\int_{0}^{\infty}g(t+u)du}{\int_{0}^{\infty}g(u)du}.$

We note that the $\tau$-periodic function $c$ can be represented via its
corresponding Fourier series representation, which is useful if a parametric
model for $c$ is required.

######  Remark 4.

The corresponding results for a trawl process follow by setting $p\equiv
c\equiv 1$.

####  2.3.2 Examples

Let us study some examples of periodic trawl processes.

######  Example 1.

Consider an exponential trawl function given by $g(x)=\exp(-\lambda x)$ for
$\lambda>0,x\geq 0$ and the $\tau$-periodic function $p(x)=\sin(2\pi x/\tau)$.
Then

$\displaystyle\mathbb{E}(Y_{t})$
$\displaystyle=\mathbb{E}(L^{\prime})\frac{2\pi\tau}{\lambda^{2}\tau^{2}+4\pi^{2}},$
$\displaystyle\mathrm{Var}(Y_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})\frac{8\pi^{2}}{\lambda^{3}\tau^{2}+16\pi^{2}\lambda},$
$\displaystyle\mathrm{Cov}(Y_{0},Y_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})2\pi\frac{e^{-\lambda
t}\left[\sin\left(\frac{2\pi
t}{\tau}\right)\tau\lambda+4\pi\cos\left(\frac{2\pi
t}{\tau}\right)\right]}{\lambda(\lambda^{2}\tau^{2}+16\pi^{2})},$
$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})$ $\displaystyle=e^{-\lambda
t}c(t),\quad\mathrm{with}\quad c(t)=\frac{1}{4\pi}\left[\sin\left(\frac{2\pi
t}{\tau}\right)\tau\lambda+4\pi\cos\left(\frac{2\pi t}{\tau}\right)\right],$

where $c$ is a $\tau$-periodic function.

Next, we consider an example based on the supGamma trawl function, which can
allow for both short and long memory.

######  Example 2.

Consider the supGamma trawl function defined as
$g(x)=\left(1+\frac{x}{\alpha}\right)^{-H}$, for $\alpha>0,H>1,x\geq 0$. As
before, let $p(x)=\sin(2\pi x/\tau)$. Then

$\displaystyle\mathrm{Cov}(Y_{0},Y_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}p(u)p(t+u)g(t+u)du$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}\sin\left(\frac{2\pi
u}{\tau}\right)\sin\left(\frac{2\pi(t+u)}{\tau}\right)\left(1+\frac{(t+u)}{\alpha}\right)^{-H}du.$

We note that
$\int_{0}^{\infty}g(t+u)du=\int_{0}^{\infty}\left(1+\frac{(t+u)}{\alpha}\right)^{-H}du=\frac{\alpha^{H}}{H-1}(t+\alpha)^{1-H}$.
According to Proposition 6, there exists a $\tau$-periodic function $c$ such
that

$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})$
$\displaystyle=c(t)\frac{\int_{0}^{\infty}g(t+u)du}{\int_{0}^{\infty}g(u)du}=c(t)\left(1+\frac{t}{\alpha}\right)^{1-H}.$

For $H\in(1,2]$, we are in the long-memory case and for $H>2$ in the short-
memory case.

######  Example 3.

Let us briefly consider the general case of a superposition trawl. Let

$\displaystyle g(x)=\int_{0}^{\infty}e^{-\lambda
x}f_{\lambda}(\lambda)d\lambda,$

for a density $f_{\lambda}$. If

$\displaystyle
f_{\lambda}(\lambda)=\frac{1}{\Gamma(H)}\alpha^{H}\lambda^{H-1}e^{-\lambda\alpha},$

for $\alpha>0,H>1$. Then

$\displaystyle g(x)=\int_{0}^{\infty}e^{-\lambda
x}\frac{1}{\Gamma(H)}\alpha^{H}\lambda^{H-1}e^{-\lambda\alpha}d\lambda=\left(1+\frac{x}{\alpha}\right)^{-H},$

i.e. we are in the case of the previous example. Then

$\displaystyle\mathrm{Cov}(Y_{0},Y_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}p(u)p(t+u)g(t+u)du,$

where

$\displaystyle\int_{0}^{\infty}p(u)p(t+u)g(t+u)du$
$\displaystyle=\int_{0}^{\infty}p(u)p(t+u)\int_{0}^{\infty}e^{-\lambda(t+u)}f_{\lambda}(\lambda)d\lambda
du$ $\displaystyle=\int_{0}^{\infty}e^{-\lambda
t}\int_{0}^{\infty}p(u)p(t+u)e^{-\lambda u}duf_{\lambda}(\lambda)d\lambda$
$\displaystyle=\int_{0}^{\infty}2\pi\frac{e^{-\lambda
t}\left[\sin\left(\frac{2\pi
t}{\tau}\right)\tau\lambda+4\pi\cos\left(\frac{2\pi
t}{\tau}\right)\right]}{\lambda(\lambda^{2}\tau^{2}+16\pi^{2})}f_{\lambda}(\lambda)d\lambda$
$\displaystyle=2\pi\tau\sin\left(\frac{2\pi
t}{\tau}\right)\int_{0}^{\infty}\frac{e^{-\lambda
t}}{(\lambda^{2}\tau^{2}+16\pi^{2})}f_{\lambda}(\lambda)d\lambda+8\pi^{2}\cos\left(\frac{2\pi
t}{\tau}\right)\int_{0}^{\infty}\frac{e^{-\lambda
t}}{\lambda(\lambda^{2}\tau^{2}+16\pi^{2})}f_{\lambda}(\lambda)d\lambda.$

From Proposition 6, we deduce that there exists a $\tau$-periodic function $c$
such that

$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})=c(t)\frac{\int_{0}^{\infty}g(t+u)du}{\int_{0}^{\infty}g(u)du}.$

Note that

$\displaystyle\int_{0}^{\infty}g(t+u)du$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}e^{-\lambda(t+u)}f_{\lambda}(\lambda)d\lambda
du=\int_{0}^{\infty}e^{-\lambda
t}f_{\lambda}(\lambda)\int_{0}^{\infty}e^{-\lambda u}dud\lambda$
$\displaystyle=\int_{0}^{\infty}\frac{1}{\lambda}e^{-\lambda
t}f_{\lambda}(\lambda)d\lambda=\mathbb{E}(e^{-\Lambda t}/\Lambda),$

where $\Lambda$ is a continuous random variable with density $f_{\lambda}$.
Also,

$\displaystyle\int_{0}^{\infty}g(u)du$
$\displaystyle=\int_{0}^{\infty}\frac{1}{\lambda}f_{\lambda}(\lambda)d\lambda=\mathbb{E}(1/\Lambda).$

That is, we have

$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})=c(t)\frac{\mathbb{E}(e^{-\Lambda
t}/\Lambda)}{\mathbb{E}(1/\Lambda)}.$

##  3 Simulation of periodic trawl processes

We will now address the question of how a periodic trawl process can be
simulated efficiently. A natural first choice might be to develop a grid-based
method, using a suitable spatio-temporal grid. However, it turns out that we
can adapt a _slice-based_ simulation algorithm which was developed for trawl
processes, see Noven (2016) and also Leonte & Veraart (2022), and extend it to
allow for any additional temporal kernel function, which in our case, will be
a $\tau$-periodic function $p$.

The advantage of this slice-based approach is that we can have a coarser, but
exact, approximation in the spatial domain, and the simulation error only
appears through the discretisation in time, whereas a grid-based approach
would result in a discretisation error both in the spatial and the temporal
domain.

Suppose that we would like to simulate $Y$ on the grid $t_{0},\ldots,t_{n}$
with $t_{i}=i\Delta$ for $i=0,\ldots,n$ for $\Delta>0$ and $t_{n}=n\Delta=T$.

###  3.1 Slice-based simulation for trawl processes

The idea behind the slicing, see Noven (2016); Leonte & Veraart (2022), is
that we consider the partition of $\cup_{i=0}^{n}A_{t_{i}}$, where
$A_{t_{i}}=\\{(x,s):s\leq t,0\leq x\leq g(t-s)\\}$, which is obtained when
considering all the disjoint sets obtained from the intersections of the
various trawl sets. We illustrate this idea in Figure 1 for the case when
$n=4$ and spell out the mathematical details in Algorithm 1.

Figure 1: Slices for $n=4$.

We create a $(n+1)\times(n+1)$-dimensional slice matrix $S=(S_{ij})$. Capital
letters refer to slices themselves, which are Borel sets of finite Lebesgue
measure, and small letters $s=(s_{ij})$ to the Lebesgue measures of the
associated slices ($s_{ij}=Leb(S_{ij})$). We will then draw independent random
variables $L(S_{i,j})$ whose distribution has cumulant function
$s_{i,j}C(\theta;L^{\prime})$ and simulate periodic trawl processes of the
form
$Y_{t}=\int_{(-\infty,t]\times\mathbb{R}}p(t-s)\mathbb{I}_{(0,g(t-s))}(x)L(dx,ds)$,
using the following approximation, for $k=0,\ldots,n$,

$\displaystyle Y_{k\Delta_{n}}$
$\displaystyle=\sum_{j=1}^{k+1}p((k+2-j)\Delta_{n})\sum_{i=k+2-j}^{n+2-j}L(S_{i,j}).$

We will now describe the simulation algorithm for periodic trawl processes in
more detail by splitting the task at hand into three algorithms.

####  3.1.1 Computing the matrix of slices

We only consider slices $S_{ij}$ (or $s_{ij}$) with $j\leq n+2-i$ and store
them in the matrix of slices of the form

$\displaystyle\begin{pmatrix}S_{11}&S_{12}&\dots&S_{1n}&S_{1(n+1)}\\\
S_{21}&S_{22}&\dots&S_{2n}&\hrectangle\\\
\vdots&\vdots&\iddots&\iddots&\vdots\\\
S_{n1}&S_{n2}&\hrectangle\dots&\dots&\hrectangle\\\
S_{(n+1)1}&\hrectangle&\dots&\dots&\hrectangle\end{pmatrix},$

see also Algorithm 1.

Algorithm 1 Computing the slices $(s_{i,j})$, for $i\leq n+1,j\leq n+2-i$,
using the trawl function $g$, $n$, and the grid width $\Delta_{n}$ as inputs.

1:procedure Slices($g$, $n$, $\Delta_{n}$)

2: $b\leftarrow\mathrm{numeric}(n+1)$ $\triangleright$ Create four vectors of
zeros

3: $c\leftarrow\mathrm{numeric}(n)$

4: $d\leftarrow\mathrm{numeric}(n+1)$

5: $e\leftarrow\mathrm{numeric}(n)$

6: for k in 1:(n+1) do

7: $b[k]\leftarrow\int_{0}^{\Delta_{n}}g(k\Delta_{n}-x)dx$

8: $d[k]\leftarrow\int_{-\infty}^{0}g(x-(k-1)\Delta_{n})dx$

9: end for

10: for k in 1:n do

11: $c[k]\leftarrow b[k]-b[k+1]$

12: $e[k]=d[k]-d[k+1]$

13: end for

14: $s\leftarrow\mathrm{matrix}(0,n+1,n+1)$ $\triangleright$ Create slice
matrix of zeros

15: for k in 1:n do

16: $s[k,1:(n+1-k)]\leftarrow\mathrm{replicate}(c[k],n+1-k)$

17: $s[k,(n+1-k+1)]\leftarrow b[k]$

18: end for

19: $s[,1]\leftarrow c(e,d[n+1])$ $\triangleright$ Compute the column of first
slices

20: return $s$ $\triangleright$ Return the slice matrix

21:end procedure

####  3.1.2 Adding the weighted slices

Algorithm 2 Adding the weighted slices given a slicematrix $s=(s_{i,j})$, for
$i\leq n+1,j\leq n+2-i$ and an $n$-dimensional weight vector $w$.

1:procedure AddWeightedSlices($s,w$)

2: $n\leftarrow\mathrm{nrow}(s)-1$

3: $x\leftarrow\mathrm{numeric}(n+1)$$\triangleright$ Create vector of zeros

4: $tmp\leftarrow 0$

5: for k in 0:n do

6: $tmp\leftarrow 0$

7: for j in 1:(k+1) do

8: $tmp\leftarrow tmp+w[k+2-j]\cdot\mathrm{sum}(s[(k+2-j):(n+2-j),j])$

9: end for

10: $x[1+k]\leftarrow tmp$

11: end for

12: return $x$

13:end procedure

####  3.1.3 Simulating a periodic trawl process

Algorithm 3 Simulating a periodic trawl process on the grid
$0,\Delta_{n},\ldots,n\Delta_{n}$ given the functions $p,g$, the distribution
of $L^{\prime}$, and the grid width $\Delta_{n}$ and $n$.

1:procedure SimulatePeriodicTrawl($p,g,C(\cdot;L^{\prime}),\Delta_{n},n$)

2: $s\leftarrow\textsc{Slices}(g,n,\Delta_{n})$ $\triangleright$ Compute the
slice-matrix $(s_{i,j})$

3: $L\leftarrow\mathrm{matrix}(n+1,n+1)$ $\triangleright$ Create matrix of
r.v.s $L(S_{ij})$

4: for k in 1:n do

5: $L[k,1:(n+1-k)]\leftarrow$ vector of $(n+1-k)$ i.i.d. r.v.s $\sim
s[k,2]\cdot C(\cdot;L^{\prime})$

6: $L[k,(n+1-k+1)]\leftarrow$ r.v. $\sim s[k,(n+1-k+1)]\cdot
C(\cdot;L^{\prime})$

7: end for

8: for k in 1:(n+1) do

9: $L[k,1]\leftarrow$ r.v. $\sim s[k,1]\cdot C(\cdot;L^{\prime})$

10: end for

11: $w\leftarrow\mathrm{numeric}(n+1)$ $\triangleright$ Create weight vector

12: for k in 1:(n+1) do

13: $w[k]\leftarrow p(k\cdot\Delta_{n})$

14: end for

15: $y\leftarrow\textsc{AddWeightedSlices}(L,w)$

16: return $y$ $\triangleright$ Returns
$Y_{0},Y_{\Delta_{n}},\ldots,Y_{n\Delta_{n}}$

17:end procedure

We note that in the above algorithm, Algorithm 3, the notation ”r.v. $\sim
s[k,(n+1-k+1)]\cdot C(\cdot;L^{\prime})$” is a short-hand notation for a
random variable whose distribution is characterised by the cumulant function
$s[k,(n+1-k+1)]\cdot C(\cdot;L^{\prime})$.

We note that the approximation for the first column of the slice matrix is
rather rough. Hence, in practice, we should use a burn-in period in the
simulation to minimise the effect of the initial approximation error.

The simulation algorithm has been implemented in the R package ambit available
on CRAN, see Veraart (2022).

###  3.2 A note on stochastic versus deterministic seasonality

We note that periodic trawl processes are stationary with stochastic
seasonality, which is reflected in the time-invariant mean, variance, and
autocorrelation, where the autocorrelation function is a periodic function.

One might wonder how such periodic trawl processes compare to trawl processes
with deterministic seasonality. As before, let $L$ denote a homogeneous Lévy
basis with characteristic triplet $(\zeta,a,\ell)$,
$p:[0,\infty)\to\mathbb{R}$ a periodic function with period $\tau>0$ and
$g:[0,\infty)\to\mathbb{R}$ a continuous, monotonically decreasing function.

Consider the trawl process

$\displaystyle
X_{t}=\int_{(-\infty,t]\times\mathbb{R}}\mathbb{I}_{(0,g(t-s))}(x)L(dx,ds)$

and the periodic trawl process

$\displaystyle
Y_{t}=\int_{(-\infty,t]\times\mathbb{R}}p(t-s)\mathbb{I}_{(0,g(t-s))}(x)L(dx,ds).$

Let $q:[0,\infty)\to\mathbb{R}$ a periodic function with period $\tau>0$. We
consider trawl processes with additive and multiplicative seasonality.

$\displaystyle X^{a}(t):=q(t)+X_{t},$ $\displaystyle X^{m}(t):=q(t)X_{t}.$

Then $X^{a},X^{m}$ are not stationary and their second-order properties are
given by

$\displaystyle\mathbb{E}(X^{a}_{t})$
$\displaystyle=q(t)+\mathbb{E}(L^{\prime})\int_{0}^{\infty}g(u)du=q(t)+\mathbb{E}(X_{0})$
$\displaystyle\mathrm{Var}(X^{a}_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}g(u)du=\mathrm{Var}(X_{0}),$
$\displaystyle\mathrm{Cov}(X^{a}_{0},X^{a}_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}g(t+u)du=\mathrm{Cov}(X_{0},X_{t}),$
$\displaystyle\mathrm{Cor}(X^{a}_{0},X^{a}_{t})$
$\displaystyle=\frac{\int_{0}^{\infty}g(t+u)du}{\int_{0}^{\infty}g(u)du}=\mathrm{Cor}(X_{0},X_{t}),$

for $t\geq 0$. Also,

$\displaystyle\mathbb{E}(X^{m}_{t})$
$\displaystyle=q(t)\mathbb{E}(L^{\prime})\int_{0}^{\infty}g(u)du=q(t)\mathbb{E}(X_{0}),$
$\displaystyle\mathrm{Var}(X^{m}_{t})$
$\displaystyle=q^{2}(t)\mathrm{Var}(L^{\prime})\int_{0}^{\infty}g(u)du=q^{2}(t)\mathrm{Var}(X_{0}),$
$\displaystyle\mathrm{Cov}(X^{m}_{0},X^{m}_{t})$
$\displaystyle=q(0)q(t)\mathrm{Var}(L^{\prime})\int_{0}^{\infty}g(t+u)du=q(0)q(t)\mathrm{Cov}(X_{0},X_{t}),$
$\displaystyle\mathrm{Cor}(X^{m}_{0},X^{m}_{t})$
$\displaystyle=\frac{q(0)q(t)\int_{0}^{\infty}g(t+u)du}{q^{2}(t)\int_{0}^{\infty}g(u)du}=\frac{q(0)}{q(t)}\mathrm{Cor}(X_{0},X_{t}).$

######  Example 4.

Consider the case where $L^{\prime}\sim\mathrm{N}(0,1)$, $p(x)=q(x)=\sin(2\pi
x/\tau)$, for $\tau=3$, $g(x)=\exp(-0.5x)$. We simulate the processes $X$,
$X^{a}$, $X^{m}$, $Y$ and the function $q$ on the grid $t_{i}=i\Delta_{n}$,
for $i=0,\ldots,n=499$ and $\Delta_{n}=0.1$.

(a) Sample path of $X$ (b) Empirical ACF of $X$

(c) Sample path of $X^{a}$ (d) Empirical ACF of $X^{a}$

(e) Sample path of $X^{m}$ (f) Empirical ACF of $X^{m}$

(g) Sample path of $Y$ (h) Empirical ACF of $Y$

(i) Function $q()$ (j) Empirical ACF of $q()$

Figure 2: Comparing the sample paths and empirical autocorrelation functions
of the stationary trawl process $X$, the non-stationary processes $X^{a}$ and
$X^{m}$, the stationary periodic trawl process $Y$ and the seasonality
function $q$.

We visually compare the sample paths and empirical autocorrelation functions
of the stationary trawl process $X$, the nonstationary processes $X^{a}$ and
$X^{m}$, the stationary periodic trawl process $Y$ and the seasonality
function $q$ in Figure 2, see Veraart (2023) for the corresponding R code.
Here we consider Gaussian processes, and the stochastic noise terms in each
simulation of the various paths were kept the same to simplify the comparison.

We note that the sample paths and empirical autocorrelation functions of $Y$
and $X^{m}$ look very similar, which suggests that it might be hard to
distinguish between these two models in practice. Note however that the key
probabilistic difference between the two processes is that $Y$ is stationary
whereas $X^{m}$ is not.

##  4 Asymptotic theory for the sample mean, the sample autocovariances and
the sample autocorrelations of mixed moving average processes

Since MMAPs and hence also periodic trawl processes are mixing and, hence,
ergodic, see Fuchs & Stelzer (2013), we know that their corresponding moment
estimators are consistent. In order to develop a suitable estimation theory
for (periodic) trawl processes, we first derive the results in the more
general setting of MMAPs.

Consider an MMAP given by $Y=(Y_{t})_{t\in\mathbb{R}}$ for

$\displaystyle Y_{t}=\mu+\int_{\mathbb{R}\times\mathbb{R}}f(x,t-s)L(dx,ds),$

for $\mu\in\mathbb{R}$, where $f$ satisfies the integrability conditions
stated in (7) almost surely.

We are interested in the asymptotic behaviour of the sample mean and sample
autocovariances/-correlations for such processes.

Note that for all $t\in\mathbb{R}$,

$\displaystyle\mathbb{E}(Y_{t})=\mu+\mathbb{E}(L^{\prime})\int_{\mathbb{R}\times\mathbb{R}}f(x,u)dxdu,$

and the autocovariance function of $Y$ is denoted by

$\displaystyle\gamma(h)=\gamma_{f}(h)=\mathrm{Cov}(Y_{0},Y_{h})=\kappa_{2}\int_{\mathbb{R}\times\mathbb{R}}f(x,-s)f(x,h-s)dxds,$

for any $h\in\mathbb{R}$, where $\kappa_{2}=\mathrm{Var}(L^{\prime})$.

Suppose that the process $Y$ is sampled over a fixed-time grid of width
$\Delta>0$ at times $(n\Delta)_{n\in\mathbb{N}}$.

The proofs of the following results are presented in the Appendix, see A.2.

###  4.1 Asymptotic normality of the sample mean

We denote the sample mean by

$\displaystyle\overline{Y}_{n;\Delta}:=\frac{1}{n}\sum_{i=1}^{n}Y_{i\Delta}.$

By adapting Cohen & Lindner (2013, Theorem 2.1) to our more general setting,
we get the following asymptotic result for the sample mean.

######  Theorem 6.

Suppose that
$\mathbb{E}(L^{\prime})=0,\kappa_{2}=\mathrm{Var}(L^{\prime})<\infty,\mu\in\mathbb{R}$
and $\Delta>0$. Further, assume that

$\displaystyle\left(F_{\Delta}:\mathbb{R}\times[0,\Delta]\to[0,\infty],(x,u)\mapsto
F_{\Delta}(x,u)=\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)|\right)\in
L^{2}(\mathbb{R}\times[0,\Delta]).$ (9)

Then $\sum_{j=-\infty}^{\infty}|\gamma(\Delta j)|<\infty$,

$\displaystyle V_{\Delta}:=\sum_{j=-\infty}^{\infty}\gamma(\Delta
j)=\kappa_{2}\int_{\mathbb{R}\times[0,\Delta]}\left(\sum_{j=1}^{\infty}f(x,u+j\Delta)\right)^{2}dxdu,$
(10)

and the sample mean of $Y_{\Delta i}$, for $i=1,\ldots,n$, is asymptotically
Gaussian as $n\to\infty$, i.e.

$\displaystyle\sqrt{n}\left(\overline{Y}_{n;\Delta}-\mu\right)\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}\left(0,V_{\Delta}\right),\quad\mathrm{as}\,n\to\infty.$

######  Remark 5.

As in Cohen & Lindner (2013), we remark that the assumption that
$\mathbb{E}(L^{\prime})=0$ can be removed if the kernel function $f$ satisfies
$f\in L^{1}(\mathbb{R}\times\mathbb{R})\cap
L^{2}(\mathbb{R}\times\mathbb{R})$.

###  4.2 Asymptotic normality of the sample autocovariance

We define the sample autocovariance function as

$\displaystyle\widehat{\gamma}_{n;\Delta}(\Delta
h):=n^{-1}\sum_{j=1}^{n-h}(Y_{j\Delta}-\overline{Y}_{n,\Delta})(Y_{(j+h)\Delta}-\overline{Y}_{n,\Delta}),\quad
h\in\\{0,\ldots,n-1\\},$

and the sample autocorrelation function as

$\displaystyle\widehat{\rho}_{n;\Delta}(\Delta
h)=\frac{\widehat{\gamma}_{n;\Delta}(\Delta
h)}{\widehat{\gamma}_{n;\Delta}(0)},\quad h\in\\{0,\ldots,n-1\\}.$

In the case when $\mu=0$ (and we assumed that $\mathbb{E}(L^{\prime})=0$), one
could use the following simpler estimators

$\displaystyle\widehat{\gamma}_{n;\Delta}^{*}(\Delta
h):=n^{-1}\sum_{j=1}^{n}Y_{j\Delta}Y_{(j+h)\Delta},\quad
h\in\\{0,\ldots,n-1\\},$

and

$\displaystyle\widehat{\rho}_{n;\Delta}^{*}(\Delta
h)=\frac{\widehat{\gamma}_{n;\Delta}^{*}(\Delta
h)}{\widehat{\gamma}_{n;\Delta}^{*}(0)},\quad h\in\\{0,\ldots,n-1\\}.$

Consider the case where $\mu=0$, for which we will study the limiting
behaviour of $\mathrm{Cov}(\widehat{\gamma}_{n,\Delta}^{*}(\Delta
p),\widehat{\gamma}_{n,\Delta}^{*}(\Delta q))$. For this, we need a formula
for the fourth (joint) moments of the mixed moving average process, which we
derive in the following lemma.

######  Lemma 1.

Let $L^{\prime}$ be a non-zero Lévy seed with $\mathbb{E}(L^{\prime})=0$ and
$\mathbb{E}(L^{\prime 4})<\infty$. Let $\kappa_{2}:=\mathbb{E}(L^{\prime
2})=\mathrm{Var}(L^{\prime})$, $\eta:=\mathbb{E}(L^{\prime 4})/\kappa_{2}^{2}$
and $\kappa_{4}:=(\eta-3)\kappa_{2}^{2}$. Let $\Delta>0$, and assume that
$f\in L^{2}(\mathbb{R}\times\mathbb{R})\cap
L^{4}(\mathbb{R}\times\mathbb{R})$. Then, for
$t_{1},t_{2},t_{3},t_{4}\in\mathbb{R}$, we have

$\displaystyle\mathbb{E}(Y_{t_{1}}Y_{t_{2}}Y_{t_{3}}Y_{t_{4}})$
$\displaystyle=\gamma(t_{1}-t_{2})\gamma(t_{3}-t_{4})+\gamma(t_{1}-t_{3})\gamma(t_{2}-t_{4})+\gamma(t_{1}-t_{4})\gamma(t_{2}-t_{3})$
$\displaystyle\quad+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-t_{3}+s)f(x,t_{2}-t_{3})f(x,s)f(x,t_{4}-t_{3}+s)dxds.$

The following result is an extension of Cohen & Lindner (2013, Proposition
3.1).

######  Proposition 7.

Let $\mu=0$. Let $L^{\prime}$ be a non-zero Lévy seed with
$\mathbb{E}(L^{\prime})=0$ and $\mathbb{E}(L^{\prime 4})<\infty$. Let
$\kappa_{2}:=\mathbb{E}(L^{\prime 2})=\mathrm{Var}(L^{\prime})$ and
$\eta:=\mathbb{E}(L^{\prime 4})/\kappa_{2}^{2}$. Let $\Delta>0$, and assume
that $f\in L^{2}(\mathbb{R}\times\mathbb{R})\cap
L^{4}(\mathbb{R}\times\mathbb{R})$ and

$\displaystyle\left(\mathbb{R}\times[0,\Delta]\to\mathbb{R},(x,u)\mapsto\sum_{j=-\infty}^{\infty}f^{2}(x,u+j\Delta)\right)\in
L^{2}(\mathbb{R}\times[0,\Delta]).$ (11)

For $q\in\mathbb{Z}$, define the function

$\displaystyle\left(G_{q;\Delta}:\mathbb{R}\times[0,\Delta]\to\mathbb{R},(x,u)\mapsto
G_{q;\Delta}(x,u)=\sum_{j=-\infty}^{\infty}f(x,u+j\Delta)f(x,u+(j+q)\Delta)\right),$
(12)

which is in $L^{2}(\mathbb{R}\times[0,\Delta])$ due to (11).

Further, assume that

$\displaystyle\sum_{j=-\infty}^{\infty}|\gamma(\Delta j)|^{2}<\infty.$ (13)

Then, for each $p,q\in\mathbb{N}$, we have

$\displaystyle\begin{split}&\lim_{n\to\infty}n\mathrm{Cov}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p),\widehat{\gamma}_{n;\Delta}^{*}(\Delta q))=v_{pq;\Delta},\\\
v_{pq;\Delta}&:=(\eta-3)\kappa_{2}^{2}\int_{\mathbb{R}\times[0,\Delta]}G_{p;\Delta}(x,u)G_{q;\Delta}(x,u)dxdu\\\
&+\sum_{l=-\infty}^{\infty}[\gamma(l\Delta)\gamma((l+p-q)\Delta)+\gamma((l-q)\Delta)\gamma((l+p)\Delta)].\end{split}$
(14)

We now state a joint central limit theorem for the sample autocovariance,
autocorrelations and their counterparts when $\mu=0$.

######  Theorem 7.

  1. 1.

Suppose that the assumptions of Proposition 7 hold and that, in addition,

$\displaystyle\sum_{j=-\infty}^{\infty}\left(\int_{\mathbb{R}\times\mathbb{R}}|f(x,u)||f(x,u+j\Delta)|dxdu\right)^{2}<\infty.$
(15)

Then, for each $h\in\mathbb{N}$, we have

$\displaystyle\sqrt{n}(\gamma_{n;\Delta}^{*}(0)-\gamma(0),\ldots,\gamma_{n;\Delta}^{*}(h\Delta)-\gamma(h\Delta))^{\top}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,V_{\Delta}),\quad
n\to\infty,$

where the asymptotic covariance matrix is given by
$V_{\Delta}=(v_{pq,\Delta})_{p,q=0,\ldots,h}\in\mathbb{R}^{h+1,h+1}$ with
$v_{pq;\Delta}$ defined as in (LABEL:eq:sumT).

  2. 2.

Suppose that the same assumptions as in 1.) hold and further assume that the
function $(x,u)\mapsto\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)|$ is in
$L^{2}(\mathbb{R}\times[0,\Delta])$. Then, for each $h\in\mathbb{N}$, we have

$\displaystyle\sqrt{n}(\widehat{\gamma}_{n;\Delta}(0)-\gamma(0),\ldots,\widehat{\gamma}_{n;\Delta}(h\Delta)-\gamma(h\Delta))^{\top}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,V_{\Delta}),\quad
n\to\infty,$

where the asymptotic covariance matrix is given by
$V_{\Delta}=(v_{pq;\Delta})_{p,q=0,\ldots,h}\in\mathbb{R}^{h+1,h+1}$ with
$v_{pq;\Delta}$ defined as in (LABEL:eq:sumT).

  3. 3.

Suppose that the same assumptions as in 1.) hold and that $f$ is not almost
everywhere equal to zero. Then, for each $h\in\mathbb{N}$, we have

$\displaystyle\sqrt{n}(\rho^{*}_{n;\Delta}(\Delta)-\rho(\Delta),\ldots,\rho^{*}_{n;\Delta}(h\Delta)-\rho(h\Delta))^{\top}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,W_{\Delta}),\quad
n\to\infty,$

where the asymptotic covariance matrix is given by
$W_{\Delta}=(w_{pq;\Delta})_{p,q=0,\ldots,h}\in\mathbb{R}^{h,h}$ with

$\displaystyle w_{pq;\Delta}$
$\displaystyle=(v_{pq;\Delta}-\rho(p\Delta)v_{0q;\Delta}-\rho(q\Delta)v_{p0;\Delta}+\rho(p\Delta)\rho(q\Delta)v_{00;\Delta})/\gamma^{2}(0)$
$\displaystyle=\frac{(\eta-3)\kappa_{2}^{2}}{\gamma^{2}(0)}\int_{\mathbb{R}\times[0,\Delta]}(G_{p;\Delta}(x,u)-G_{0;\Delta}(x,u)\rho(p\Delta))$
$\displaystyle\quad\cdot(G_{q;\Delta}(x,u)-G_{0;\Delta}(x,u)\rho(q\Delta))dxdu$
$\displaystyle+\sum_{l=-\infty}^{\infty}\left[\rho((l+q)\Delta)\rho((l+p)\Delta)+\rho((l-q)\Delta)\rho((l+p)\Delta)\right.$
$\displaystyle-2\rho((l+q)\Delta)\rho(l\Delta)\rho(p\Delta)-2\rho(l\Delta)\rho((l+p)\Delta)\rho(q\Delta)$
$\displaystyle\left.+2\rho(p\Delta)\rho(q\Delta)\rho^{2}(l\Delta)\right].$

If, in addition, the function
$(x,u)\mapsto\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)|$ is in
$L^{2}(\mathbb{R}\times[0,\Delta])$, then we have

$\displaystyle\sqrt{n}(\widehat{\rho}_{n,\Delta}(\Delta)-\rho(\Delta),\ldots,\widehat{\rho}_{n,\Delta}(h\Delta)-\rho(h\Delta))^{\top}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,W_{\Delta}),\quad
n\to\infty,$

with $W_{\Delta}$ defined as above.

######  Remark 6.

Note that the assumptions used in the above theorems typically rule out long-
memory settings. We discuss the detailed implications in the context of
(periodic) trawl processes in the appendix, see Section A.4.

##  5 Inference for periodic trawl processes using methods of moments

Based on two working examples, we will now develop an estimation and inference
methodology for the parameters determining the serial correlation of
(periodic) trawl processes using a method-of-moments approach and present the
corresponding asymptotic theory. Our methodology can in principle be applied
to general trawl functions, but in order to simplify the exposition, we will
concentrate on two particularly relevant specifications, namely the
exponential trawl and the supGamma trawl function.

###  5.1 Exponential trawl function

Recall that, in the case of an exponential trawl function $g(x)=\exp(-\lambda
x),\lambda>0,x\geq 0$, the autocorrelation is, for all $t>0$, given by

$\displaystyle\rho(t)=\mathrm{Cor}(Y_{0},Y_{t})=c(t)\exp(-\lambda t),$

which is equivalent to

$\displaystyle\lambda=\log\left(\frac{\rho(t)}{c(t)}\right)/(-t).$

Setting $t=\Delta$ and assuming (for now) that $c(\Delta)$ is known (for a
(non-periodic) trawl process, we would have $c\equiv 1$), we get the following
estimator

$\displaystyle\hat{\lambda}=-\log\left(\frac{\rho_{n;\Delta}^{*}(\Delta)}{c(\Delta)}\right),$

where $\widehat{\rho}(\cdot)$ denotes the empirical autocorrelation function.

We now establish the corresponding limit theorem. From Theorem 7, assuming the
corresponding assumptions hold, we deduce that

$\displaystyle\sqrt{n}(\rho_{n;\Delta}^{*}(\Delta)-\rho(\Delta))\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,w_{11;\Delta}).$

Next, we define the function
$\phi:(0,\infty)\to\mathbb{R},\phi(x):=-\log(x/c(\Delta))$, which is
continuously differentiable with $\phi^{\prime}(x)=-\frac{1}{x}$ and set

$\displaystyle\widehat{\lambda}:=\phi(\rho_{n;\Delta}^{*}(\Delta)).$

An application of the delta method, see Brockwell & Davis (1987, Proposition
6.4.3), leads to

$\displaystyle\sqrt{n}(\widehat{\lambda}-\lambda)=\sqrt{n}(\phi(\rho_{n;\Delta}^{*}(\Delta))-\phi(\rho(\Delta)))\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,\Sigma_{\lambda;\Delta}),$

where

$\displaystyle\Sigma_{\lambda;\Delta}:=\phi^{\prime}(\rho(\Delta))w_{11;\Delta}\phi^{\prime}(\rho(\Delta))=\frac{w_{11;\Delta}}{\rho^{2}(\Delta)}=\frac{w_{11;\Delta}}{c^{2}(\Delta)\exp(-2\lambda\Delta)}=\frac{w_{11;\Delta}\exp(2\lambda\Delta)}{c^{2}(\Delta)}.$

######  Remark 7.

In the case of a nonperiodic trawl, we would have $c(\Delta)=1$, which
simplifies the asymptotic result. For periodic trawls, in the case of unknown
$c(\Delta)$, a plug-in estimator of the following form could be considered:
$\hat{\lambda}=-\log\left(\frac{\rho_{n;\Delta}^{*}(\Delta)}{\widehat{c(\Delta)}}\right)$.

Now, suppose that $c(\Delta)$ is not known, but the period $\tau$ of the
periodic function is known. I.e. we have that $c(x+\tau)=c(x)$, for all $x$.
Furthermore, assume that there exists a
$\widetilde{\tau}_{\Delta}\in\mathbb{N}$ such that
$\tau=\widetilde{\tau}_{\Delta}\Delta$. To simplify the notation, we shall
also define $\mathcal{T}_{\Delta}:=(1+\widetilde{\tau}_{\Delta})$. Then, we
have

$\displaystyle c(\mathcal{T}_{\Delta}\Delta)=c(\Delta).$

We can now construct an estimator for the memory parameter based on two
empirical autocorrelations. More precisely, note that

$\displaystyle\rho(\Delta)=e^{-\lambda\Delta}c(\Delta),$
$\displaystyle\rho(\mathcal{T}_{\Delta}\Delta)=e^{-\lambda\mathcal{T}_{\Delta}\Delta}c(\mathcal{T}_{\Delta}\Delta)=e^{-\lambda\mathcal{T}_{\Delta}\Delta}c(\Delta),$

which leaves us with two equations to estimate the two unknown parameters
$\lambda$ and $c(\Delta)$.

We can solve both equations for $c(\Delta)$ and obtain

$\displaystyle c(\Delta)=\rho(\Delta)e^{\lambda\Delta},$ $\displaystyle
c(\Delta)=\rho(\mathcal{T}_{\Delta}\Delta)e^{\lambda\mathcal{T}_{\Delta}\Delta}.$

Setting the two equations equal and solving for $\lambda$ leads to

$\displaystyle\lambda$
$\displaystyle=\frac{1}{\Delta(1-\mathcal{T}_{\Delta})}\log\left(\frac{\rho(\mathcal{T}_{\Delta}\Delta)}{\rho(\Delta)}\right).$

After $\lambda$ has been estimated, we can also estimate $c(l\Delta)$ for
$l=1,\ldots,\tilde{\tau}_{\Delta}$, using the relation

$\displaystyle c(l\Delta)=\rho(l\Delta)e^{\lambda l\Delta}.$

Altogether, we can consider the following estimator

$(\widehat{\lambda},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}\\\
:=\left(\frac{1}{\Delta(1-\mathcal{T}_{\Delta})}\log\left(\frac{\rho_{n;\Delta}^{*}(\mathcal{T}_{\Delta}\Delta)}{\rho_{n;\Delta}^{*}(\Delta)}\right),\left(\rho_{n;\Delta}^{*}(l\Delta)\exp\left\\{\frac{l}{(1-\mathcal{T}_{\Delta})}\log\left(\frac{\rho_{n;\Delta}^{*}(\mathcal{T}_{\Delta}\Delta)}{\rho_{n;\Delta}^{*}(\Delta)}\right)\right\\}\right)_{l=1,\ldots,\tilde{\tau}_{\Delta}}\right).$

We note that the vector
$(\widehat{\lambda},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}$
can be represented as a function of the empirical autocorrelation functions
$(\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}}$. As such,
we can apply the delta method to derive a joint central limit theorem as
follows. More precisely, we have

$\displaystyle(\widehat{\lambda},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}=F((\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}})$
$\displaystyle=(F_{1}((\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}}),\ldots,F_{\mathcal{T}_{\Delta}}((\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}}))^{\top},$

where the functions $F_{1},\ldots,F_{\mathcal{T}_{\Delta}}$ are given by

$\displaystyle F_{1}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\frac{1}{\Delta(1-\mathcal{T}_{\Delta})}\log\left(\frac{\rho_{\mathcal{T}_{\Delta}}}{\rho_{1}}\right),$
$\displaystyle F_{2}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{2}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\rho_{1}\exp(\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle
F_{3}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{3}(\rho_{1},\rho_{2},\rho_{\mathcal{T}_{\Delta}})=\rho_{2}\exp(2\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle
F_{4}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{4}(\rho_{1},\rho_{3},\rho_{\mathcal{T}_{\Delta}})=\rho_{3}\exp(3\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\vdots$
$\displaystyle
F_{\mathcal{T}_{\Delta}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{\mathcal{T}_{\Delta}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}-1},\rho_{\mathcal{T}_{\Delta}})=\rho_{\mathcal{T}_{\Delta}-1}\exp((\mathcal{T}_{\Delta}-1)\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})).$

Under the assumptions of Theorem 7,

$\displaystyle\sqrt{n}((\widehat{\lambda},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}-(\lambda,c(\Delta),\ldots,c(\mathcal{T}_{\Delta}\Delta))$

converges to a multivariate normal random vector with zero mean and variance
given by $D{W}_{\Delta}D^{\prime}$, where ${W}_{\Delta}$ is the
$\mathcal{T}_{\Delta}\times\mathcal{T}_{\Delta}$-matrix whose elements are the
ones defined in Theorem 7. The matrix $D$ is a
$\mathcal{T}_{\Delta}\times\mathcal{T}_{\Delta}$-matrix $[(\partial
F_{i}/\partial\rho_{j})(\rho(\Delta),\ldots,\rho(\mathcal{T}_{\Delta}\Delta))^{\top}]$,
where

$\displaystyle\frac{\partial
F_{1}}{\partial\rho_{1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{1}}{\partial\rho_{1}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\frac{1}{\Delta(\mathcal{T}_{\Delta}-1)\rho_{1}},$
$\displaystyle\frac{\partial
F_{1}}{\partial\rho_{2}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\cdots=\frac{\partial
F_{1}}{\partial\rho_{\mathcal{T}_{\Delta}-1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})=0,$
$\displaystyle\frac{\partial
F_{1}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{1}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\frac{1}{\Delta(1-\mathcal{T}_{\Delta})\rho_{\mathcal{T}_{\Delta}}},$
$\displaystyle\frac{\partial
F_{2}}{\partial\rho_{1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{2}}{\partial\rho_{1}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\frac{\mathcal{T}_{\Delta}}{\mathcal{T}_{\Delta}-1}\exp(\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\frac{\partial
F_{2}}{\partial\rho_{2}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\cdots=\frac{\partial
F_{2}}{\partial\rho_{\mathcal{T}_{\Delta}-1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})=0,$
$\displaystyle\frac{\partial
F_{2}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{2}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\frac{1}{1-\mathcal{T}_{\Delta}}\frac{\rho_{1}}{\rho_{\mathcal{T}_{\Delta}}}\exp(\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\frac{\partial
F_{3}}{\partial\rho_{1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{3}}{\partial\rho_{1}}(\rho_{1},\rho_{2},\rho_{\mathcal{T}_{\Delta}})=\frac{2}{\mathcal{T}_{\Delta}-1}\frac{\rho_{2}}{\rho_{1}}\exp(2\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\frac{\partial
F_{3}}{\partial\rho_{2}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{3}}{\partial\rho_{2}}(\rho_{1},\rho_{2},\rho_{\mathcal{T}_{\Delta}})=\exp(2\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\frac{\partial
F_{3}}{\partial\rho_{3}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\cdots=\frac{\partial
F_{3}}{\partial\rho_{\mathcal{T}_{\Delta}-1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})=0,$
$\displaystyle\frac{\partial
F_{3}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{3}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\rho_{2},\rho_{\mathcal{T}_{\Delta}})=\frac{2}{1-\mathcal{T}_{\Delta}}\frac{\rho_{2}}{\rho_{\mathcal{T}_{\Delta}}}\exp(2\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\vdots$
$\displaystyle\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{1}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}-1}\rho_{\mathcal{T}_{\Delta}})=\frac{\rho_{\mathcal{T}_{\Delta}-1}}{\rho_{1}}\exp((\mathcal{T}_{\Delta}-1)\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{2}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\cdots=\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{\mathcal{T}_{\Delta}-2}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})=0,$
$\displaystyle\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{\mathcal{T}_{\Delta}-1}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{\mathcal{T}_{\Delta}-1}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}-1},\rho_{\mathcal{T}_{\Delta}})=\exp((\mathcal{T}_{\Delta}-1)\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})),$ $\displaystyle\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=\frac{\partial
F_{\mathcal{T}_{\Delta}}}{\partial\rho_{\mathcal{T}_{\Delta}}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}-1},\rho_{\mathcal{T}_{\Delta}})=-\frac{\rho_{\mathcal{T}_{\Delta}-1}}{\rho_{\mathcal{T}_{\Delta}}}\exp((\mathcal{T}_{\Delta}-1)\Delta
F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})).$

###  5.2 SupGamma trawl function

Let us now focus on the case of a supGamma trawl function that allows for both
short- and long-memory settings.

We note that the inference for the long-memory case appears numerically
unstable when approximating the corresponding integrals appearing in the
autocorrelation function (in simulation experiments not reported here). Hence
we suggest, as before, using the mean value theorem result for inference on
the memory parameter. Let $g(x)=\left(1+\frac{x}{\alpha}\right)^{-H}$, for
$\alpha>0,H>1$. As before, there exists a $\tau$-periodic function $c$ such
that, for $t\geq 0$,

$\displaystyle\rho(t)=\mathrm{Cor}(Y_{0},Y_{t})$
$\displaystyle=c(t)\left(1+\frac{t}{\alpha}\right)^{1-H}.$

For $H\in(1,2]$, we are in the long-memory case and for $H>2$ in the short-
memory case.

We are interested in estimating the Hurst index $H$ and therefore in the
following assume that the parameter $\alpha$ is known.

Using the analogous procedure as in the short-memory case, we can write that,
for all $t\geq 0$,

$\displaystyle\frac{\rho(t)}{c(t)}=\exp((1-H)\log(1+t/\alpha))\Leftrightarrow\log\left(\frac{\rho(t)}{c(t)}\right)=(1-H)\log(1+t/\alpha).$

This suggests the following estimator of the Hurst index:

$\displaystyle\widehat{H}=1-\log\left(\frac{\rho_{n;\Delta}^{*}(\Delta)}{c(\Delta)}\right)/\log(1+\Delta/\alpha),$

assuming that $c(\Delta)$ and $\alpha$ are known.

Now, define the function
$\phi:(0,\infty)\to\mathbb{R},\phi(x):=1-\log\left(\frac{x}{c(\Delta)}\right)/\log(1+\Delta/\alpha)$,
which is continuously differentiable with
$\phi^{\prime}(x)=-\frac{1}{\log(1+\Delta/\alpha)}\frac{1}{x}$ and set

$\displaystyle\widehat{H}:=\phi(\rho_{n;\Delta}^{*}(\Delta)).$

An application of the delta method, see Brockwell & Davis (1987, Proposition
6.4.3), leads to

$\displaystyle\sqrt{n}(\widehat{H}-H)=\sqrt{n}(\phi(\rho_{n;\Delta}^{*}(\Delta))-\phi(\rho(\Delta)))\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}(0,\Sigma_{H;\Delta}),$

where

$\displaystyle\Sigma_{H;\Delta}$
$\displaystyle:=\phi^{\prime}(\rho(\Delta))w_{11;\Delta}\phi^{\prime}(\rho(\Delta))=\frac{w_{11;\Delta}}{(\log(1+\Delta/\alpha))^{2}\rho^{2}(\Delta)}$
$\displaystyle=\frac{w_{11;\Delta}}{(\log(1+\Delta/\alpha))^{2}c^{2}(\Delta)\left(1+\frac{\Delta}{\alpha}\right)^{2-2H}}=\frac{w_{11;\Delta}\left(1+\frac{\Delta}{\alpha}\right)^{2H-2}}{(\log(1+\Delta/\alpha))^{2}c^{2}(\Delta)}.$

Now suppose that $\tau$ is known and proceed as in the exponential case:
Assume that there exists a $\widetilde{\tau}_{\Delta}\in\mathbb{N}$ such that
$\tau=\widetilde{\tau}_{\Delta}\Delta$ and set
$\mathcal{T}_{\Delta}:=(1+\widetilde{\tau}_{\Delta})$. Then, we have

$\displaystyle c(\mathcal{T}_{\Delta}\Delta)=c(\Delta).$

We can now construct an estimator for the memory parameter based on two
empirical autocorrelations. More precisely, note that

$\displaystyle\rho(\Delta)=(1+\Delta/\alpha)^{H-1}c(\Delta),$
$\displaystyle\rho(\mathcal{T}_{\Delta}\Delta)=(1+\mathcal{T}_{\Delta}\Delta/\alpha)^{H-1}c(\mathcal{T}_{\Delta}\Delta)=(1+\mathcal{T}_{\Delta}\Delta/\alpha)^{H-1}c(\Delta),$

which leaves us with two equations to estimate the two unknown parameters $H$
and $c(\Delta)$, assuming that $\alpha$ is known.

We can solve both equations for $c(\Delta)$ and obtain

$\displaystyle c(\Delta)=\rho(\Delta)(1+\Delta/\alpha)^{1-H},$ $\displaystyle
c(\Delta)=\rho(\mathcal{T}_{\Delta}\Delta)(1+\mathcal{T}_{\Delta}\Delta/\alpha)^{1-H}.$

Setting the two equations equal and solving for $H$ leads to

$\displaystyle
H=1+\log\left(\frac{\rho(\mathcal{T}_{\Delta}\Delta)}{\rho(\Delta)}\right)/\log\left(\frac{\alpha+\Delta}{\alpha+\mathcal{T}_{\Delta}\Delta}\right).$

After $H$ has been estimated, we can also estimate $c(l\Delta)$ for
$l=1,\ldots,\tilde{\tau}_{\Delta}$, using the relation

$\displaystyle c(l\Delta)=\rho(l\Delta)(1+l\Delta/\alpha)^{1-H}.$

Altogether, we can consider the following estimator

$(\widehat{H},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}\\\
:=\left(1+\log\left(\frac{\rho_{n;\Delta}^{*}(\mathcal{T}_{\Delta}\Delta)}{\rho_{n;\Delta}^{*}(\Delta)}\right)/\log\left(\frac{\alpha+\Delta}{\alpha+\mathcal{T}_{\Delta}\Delta}\right),\right.\\\
\left.\left(\rho_{n;\Delta}^{*}(l\Delta)(1+l\Delta/\alpha)^{\log\left(\frac{\rho_{n;\Delta}^{*}(\mathcal{T}_{\Delta}\Delta)}{\rho_{n;\Delta}^{*}(\Delta)}\right)/\log\left(\frac{\alpha+\Delta}{\alpha+\mathcal{T}_{\Delta}\Delta}\right)}\right)_{l=1,\ldots,\tilde{\tau}_{\Delta}}\right).$

We note that the vector
$(\widehat{H},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}$
can be represented as a function of the empirical autocorrelation functions
$(\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}}$. That is,
we have

$\displaystyle(\widehat{H},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}=F((\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}})$
$\displaystyle=(F_{1}((\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}}),\ldots,F_{\mathcal{T}_{\Delta}}((\rho_{n;\Delta}^{*}(l\Delta))_{l=1,\ldots,\mathcal{T}_{\Delta}}))^{\top},$

where the functions $F_{1},\ldots,F_{\mathcal{T}_{\Delta}}$ are given by

$\displaystyle F_{1}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=1+\log\left(\frac{\rho_{\mathcal{T}_{\Delta}}}{\rho_{1}}\right)/\log\left(\frac{\alpha+\Delta}{\alpha+\mathcal{T}_{\Delta}\Delta}\right),$
$\displaystyle F_{2}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{2}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})=\rho_{1}(1+\Delta/\alpha)^{1-F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})}$
$\displaystyle F_{3}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{3}(\rho_{1},\rho_{2},\rho_{\mathcal{T}_{\Delta}})=\rho_{2}(1+2\Delta/\alpha)^{1-F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})},$
$\displaystyle F_{4}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{4}(\rho_{1},\rho_{3},\rho_{\mathcal{T}_{\Delta}})=\rho_{3}(1+3\Delta/\alpha)^{1-F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})},$
$\displaystyle\vdots$ $\displaystyle
F_{\mathcal{T}_{\Delta}}(\rho_{1},\ldots,\rho_{\mathcal{T}_{\Delta}})$
$\displaystyle=F_{\mathcal{T}_{\Delta}}(\rho_{1},\rho_{\mathcal{T}_{\Delta}-1},\rho_{\mathcal{T}_{\Delta}})=\rho_{\mathcal{T}_{\Delta}-1}(1+(\mathcal{T}_{\Delta}-1)\Delta/\alpha)^{1-F_{1}(\rho_{1},\rho_{\mathcal{T}_{\Delta}})}.$

Under the assumptions of Theorem 7

$\displaystyle\sqrt{n}((\widehat{H},\widehat{c(\Delta)},\ldots,\widehat{c(\tilde{\tau}_{\Delta}\Delta)})^{\top}-(H,c(\Delta),\ldots,c(\mathcal{T}_{\Delta}\Delta))$

converges to a multivariate normal random vector with zero mean and variance
given by $D{W}_{\Delta}D^{\prime}$, where ${W}_{\Delta}$ is the
$\mathcal{T}_{\Delta}\times\mathcal{T}_{\Delta}$-matrix whose elements are the
ones defined in Theorem 7. The matrix $D$ is a
$\mathcal{T}_{\Delta}\times\mathcal{T}_{\Delta}$-matrix $[(\partial
F_{i}/\partial\rho_{j})(\rho(\Delta),\ldots,\rho(\mathcal{T}_{\Delta}\Delta))^{\top}]$.

######  Remark 8.

Note that the assumptions of Theorem 7 are not satisfied in the case of a
long-memory supGamma trawl function, i.e. when $H\in(1,2]$, see Subsection A.4
in the Appendix.

##  6 Inference for periodic trawl processes using a generalised-method-of-
moments approach

In the previous section, we showed how the memory parameter and periodic
function of a periodic trawl process can be inferred if the period is known.
More generally, if we would like to estimate all model parameters, including
the parameters of the Lévy basis, we can proceed by using the generalised
method of moments (GMM).

Hence, we will now develop the asymptotic theory for estimating the parameters
of a periodic trawl process using a generalised method of moments (GMM). This
extends the work presented in Bennedsen et al. (2023) for (integer-valued)
trawl processes to the case of periodic trawl processes.

Recall that we denote the periodic trawl process by

$\displaystyle Y_{t}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}p(t-s)\mathbb{I}_{(0,g(t-s))}(x)\mathbb{I}_{[0,\infty)}(t-s)L(dx,ds),$

and we denote the corresponding trawl process, where $p\equiv 1$, by

$\displaystyle X_{t}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}\mathbb{I}_{(0,g(t-s))}(x)\mathbb{I}_{[0,\infty)}(t-s)L(dx,ds).$
(16)

###  6.1 Weak dependence

We will first show that periodic trawl processes are $\theta$-weakly
dependent, see Curato & Stelzer (2019, Definition 3.2).

We note that a periodic trawl process is a special case of a causal mixed
moving average process as defined in Curato & Stelzer (2019, Definition 3.3).
Hence, under the assumption that $\int_{|\xi|>1}|\xi|^{2}\ell(d\xi)<\infty$,
we can deduce from Curato & Stelzer (2019, Corollary 3.4) the periodic trawl
process is $\theta$-weakly dependent, see Curato & Stelzer (2019, Definition
3.2) with coefficient, for $r\geq 0$,

$\displaystyle\theta_{Y}(r)$
$\displaystyle=\left(\mathrm{Var}(L^{\prime})\int_{(-\infty,-r)\times\mathbb{R}}p^{2}(-s)\mathbb{I}_{(0,g(-s))}^{2}(x)\mathbb{I}_{[0,\infty)}^{2}(-s)dxds\right.$
$\displaystyle\left.+\left|\mathbb{E}(L^{\prime})\int_{(-\infty,-r)\times\mathbb{R}}p(-s)\mathbb{I}_{(0,g(-s))}(x)\mathbb{I}_{[0,\infty)}(-s)dxds\right|^{2}\right)^{1/2}$
$\displaystyle=\left(\mathrm{Var}(L^{\prime})\int_{-\infty}^{-r}p^{2}(-s)g(-s)ds+(\mathbb{E}(L^{\prime}))^{2}\left(\int_{-\infty}^{-r}p(-s)g(-s)ds\right)^{2}\right)^{1/2}$
$\displaystyle=\left(\mathrm{Var}(L^{\prime})\int_{r}^{\infty}p^{2}(s)g(s)ds+(\mathbb{E}(L^{\prime}))^{2}\left(\int_{r}^{\infty}p(s)g(s)ds\right)^{2}\right)^{1/2}.$

Since the periodic function is continuous and bounded, we note that the weak-
dependence coefficient can be bounded by a function $\widetilde{\theta}_{Y}$
such that $\theta_{Y}(r)\leq\widetilde{\theta}_{Y}(r)$ for all $r\geq 0$,
where

$\displaystyle\widetilde{\theta}_{Y}(r)$
$\displaystyle=\left(c_{1}\mathrm{Cov}(X_{0},X_{r})+c_{2}(\mathbb{E}(L^{\prime}))^{2}(\mathrm{Cov}(X_{0},X_{r}))^{2}\right)^{1/2},$

where $\mathbb{E}(L^{\prime})=\gamma+\int_{|\xi|>1}\ell(d\xi)$, and
$\mathrm{Var}(L^{\prime})=a+\int_{\mathbb{R}}\xi^{2}\ell(d\xi)$ and
$c_{1},c_{2}>0$ are constants.

Let us briefly consider the special case when $L^{\prime}$ is of finite
variation, i.e. when the characteristic triplet is given by $(\gamma,0,\ell)$
with $\int_{\mathbb{R}}|\xi|\ell(\xi)<\infty$. This case is of interest since
it includes the well-known class of integer-valued trawl processes. Here, the
weak dependence coefficient is, for $r\geq 0$, given by

$\displaystyle\theta_{Y}(r)$
$\displaystyle=\int_{(-\infty,-r)\times\mathbb{R}}\int_{\mathbb{R}}|p(-s)\mathbb{I}_{(0,g(-s))}(x)\mathbb{I}_{[0,\infty)}(-s)\xi|\ell(d\xi)dxds$
$\displaystyle+\int_{(-\infty,-r)\times\mathbb{R}}|p(-s)\mathbb{I}_{(0,g(-s))}(x)\mathbb{I}_{[0,\infty)}(-s)\gamma_{0}|dxds$
$\displaystyle=\int_{\mathbb{R}}|\xi|\ell(d\xi)\int_{-\infty}^{-r}|p(-s)|g(-s)ds+|\gamma_{0}|\int_{-\infty}^{-r}|p(-s)|g(-s)ds$
$\displaystyle=\left(\int_{\mathbb{R}}|\xi|\ell(d\xi)+|\gamma_{0}|\right)\int_{r}^{\infty}|p(s)|g(s)ds.$

As before, we note that the weak-dependence coefficient can be bounded by a
function $\widetilde{\theta}_{Y}$ such that
$\theta_{Y}(r)\leq\widetilde{\theta}_{Y}(r)$ for all $r\geq 0$, where

$\displaystyle\widetilde{\theta}_{Y}(r)$
$\displaystyle=c\mathrm{Cov}(X_{0},X_{r}),$

for a positive constant $c>0$ and $\gamma_{0}=\gamma-\int_{|\xi|\leq
1}\xi\ell(d\xi)$.

######  Remark 9.

We note that the $\theta$-dependence coefficient of the periodic trawl process
$Y$ can be related to the $\theta$-dependence coefficient of $X$ via

$\displaystyle\theta_{Y}(r)\leq\widetilde{\theta}_{Y}(r)=O(\theta_{X}(r)).$

###  6.2 GMM estimation for periodic trawl processes

In this section, we will describe how the parameters of a periodic trawl
process can be estimated via the generalised method of moments (GMM).

As before, we consider the equidistantly sampled process
$Y_{\Delta},Y_{2\Delta},\ldots,Y_{n\Delta}$, for
$\Delta=T/n>0,T>0,n\in\mathbb{N}$. We will consider a GMM estimator, which is
based on the sample mean, sample variance and sample autocovariances up to lag
$m\geq 2$.

Consider the vector

$Y_{t}^{(m)}=(Y_{t\Delta},Y_{(t+1)\Delta},\ldots,Y_{(t+m)\Delta}),$

for $t=1,\ldots,n-m$. We denote by $\Theta$ the parameter space of the
periodic trawl process and set $\mu:=\mu(\theta)=\mathbb{E}(Y_{0})$ and
$D(k):=D(k,\theta):=\mathbb{E}(Y_{0}Y_{k\Delta})$, for $k=0,\ldots,m$. So, as
soon as we specify a parametric model for $Y$, then $D(k)$ is a function of
the model parameter(s) $\theta$.

Next, we define the measurable function
$h:\mathbb{R}^{m+1}\times\Theta\to\mathbb{R}^{m+2}$ by

$\displaystyle h(Y_{t}^{(m)},\theta)$
$\displaystyle=\left(\begin{array}[]{c}h_{E}(Y_{t}^{(m)},\theta)\\\
h_{0}(Y_{t}^{(m)},\theta)\\\ h_{1}(Y_{t}^{(m)},\theta)\\\ \vdots\\\
h_{m}(Y_{t}^{(m)},\theta)\\\
\end{array}\right)=\left(\begin{array}[]{c}Y_{t\Delta}-\mu(\theta)\\\
Y_{t\Delta}^{2}-D(0,\theta)\\\ Y_{t\Delta}Y_{(t+1)\Delta}-D(1,\theta)\\\
\vdots\\\ Y_{t\Delta}Y_{(t+m)\Delta}-D(m,\theta)\end{array}\right).$

Moreover, define the corresponding sample moments as

$\displaystyle g_{n,m}(\theta)$
$\displaystyle=\frac{1}{n-m}\sum_{t=1}^{n-m}h(Y_{t}^{(m)},\theta)=\left(\begin{array}[]{c}\frac{1}{n-m}\sum_{t=1}^{n-m}h_{E}(Y_{t}^{(m)},\theta)\\\
\frac{1}{n-m}\sum_{t=1}^{n-m}h_{0}(Y_{t}^{(m)},\theta)\\\ \vdots\\\
\frac{1}{n-m}\sum_{t=1}^{n-m}h_{m}(Y_{t}^{(m)},\theta)\end{array}\right).$

Suppose that the true parameter (vector) is denoted by $\theta_{0}$, say,
which can be estimated by minimising the objective function of the GMM, i.e.
the GMM estimator is given by

$\displaystyle\widehat{\theta}_{0,\mathrm{GMM}}^{n,m}=\mathrm{argmin}g_{n,m}(\theta)^{\top}A_{n,m}g_{n,m}(\theta),$
(17)

for a positive-definite weight matrix $A_{n,m}$.

We now derive the weak consistency and asymptotic normality of the GMM
estimator under suitable (standard) assumptions.

First, we present the assumptions which guarantee the weak consistency of the
estimator, cf. Mátyás (1999, Assumptions 1.1-1.3)

######  Assumption 1.

  1. (i)

Suppose that the expectation $\mathbb{E}(h(Y_{t}^{(m)},\theta))$ exists and is
finite for all $\theta\in\Theta$ and for all $t$.

  2. (ii)

Set $h_{t}^{(m)}(\theta)=\mathbb{E}(h(Y_{t}^{(m)},\theta))$. There exists a
$\theta_{0}\in\Theta$ such that $h_{t}^{(m)}(\theta)=0$ for all $t$ if and
only if $\theta=\theta_{0}$.

######  Assumption 2.

Let $h^{(m)}(\theta)=\sum_{t=1}^{n-m}h_{t}^{(m)}(\theta)$ and denote the $j$th
component of the $m+2$-dimensional vectors $h^{(m)}(\theta)$ and
$g_{n,m}(\theta)$ by $h^{(m)}_{j}(\theta)$ and $g_{n,m;j}(\theta)$,
respectively. Suppose that, for $j=1,\ldots,m+2$, as $n\to\infty$,

$\sup_{\theta\in\Theta}|h^{(m)}_{j}(\theta)-g_{n,m;j}(\theta)|\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0,$

######  Assumption 3.

There exists a sequence of non-random, positive definite matrices
$\overline{A}_{n,m}$ such that, as $n\to\infty$,
$|A_{n,m}-\overline{A}_{n,m}|\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0$

From Mátyás (1999, Theorem 1.1), we deduce the following result.

######  Theorem 8.

Assume that Assumptions 1, 2, 3 hold. Then the GMM estimator
$\widehat{\theta}_{0,\mathrm{GMM}}^{n,m}$ defined in (17) is weakly
consistent.

Next, we formulate the assumptions needed for the central limit theorem.

######  Assumption 4.

$\Theta$ is a compact parameter space which includes the true parameter
$\theta_{0}$.

######  Assumption 5.

The weight matrix $A_{n,m}$ converges in probability to a positive definite
matrix $A$.

######  Assumption 6.

The covariance matrix $\Sigma_{a}$ defined in (18) below is positive definite.

######  Theorem 9.

Consider a periodic trawl process $Y$ with characteristic triplet
$(\zeta,a,\ell)$ and suppose that
$\int_{|\xi|>1}|\xi|^{4+\delta}\ell(d\xi)<\infty$, for some $\delta>0$ and and
suppose that the $\theta$-weak dependence coefficient of the periodic trawl
process satisfies $\theta_{Y}(r)\leq O(r^{-\alpha})$, for
$\alpha>\left(1+\frac{1}{\delta}\right)\left(1+\frac{1}{2+\delta}\right)$.
Suppose that Assumptions 1, 4, 5, 6 hold. Then, as $n\to\infty$,

$\sqrt{n}(\widehat{\theta}_{0,\mathrm{GMM}}^{n,m}-\theta_{0})\stackrel{{\scriptstyle
d}}{{\to}}\mathrm{N}(0,M\Sigma_{a}M^{\top}),$

where

$\displaystyle\Sigma_{a}$
$\displaystyle=\sum_{l\in\mathbb{Z}}\mathrm{Cov}(h(Y_{0}^{(m)},\theta_{0}),h(Y_{l}^{(m)},\theta_{0})),$
(18) $\displaystyle M$
$\displaystyle=(G_{0}^{\top}AG_{0})^{-1}G_{0}^{\top}A,\quad\mathrm{where}\;\;\;G_{0}=\mathbb{E}\left[\frac{\partial
h(Y_{t}^{(m)},\theta)}{\partial\theta^{\top}}\right]_{\theta=\theta_{0}}.$

###### Proof of Theorem 9.

The proof follows, with very minor modifications, the steps of the proof
presented for trawl processes in Bennedsen et al. (2023), which is based on
the arguments of the proofs of Mátyás (1999, Theorem 1.2), see also Curato &
Stelzer (2019, Proof of Theorem 6.2) for the case of a supOU process. ∎

######  Remark 10.

We note that the assumption on the $\theta$-weak dependence index rule out
long-memory settings, see Bennedsen et al. (2023) for a more detailed
discussion in the case of trawl processes.

##  7 Empirical illustration on electricity day-ahead prices

In this section, we will illustrate how the proposed methodology for
estimating the kernel function of periodic trawl processes developed in
Section 5 can be used in practice, see Veraart (2023) for the corresponding
data and R code.

We consider day-ahead electricity baseload prices for Germany and Luxembourg
from 1st October 2018 to 1st January 2023 recorded in EUR/MWh. The data have
been downloaded from the website https://www.smard.de/en, which is maintained
by the Bundesnetzagentur in Germany. The time series is depicted in Figure 5.

Given the turmoil experienced by the electricity market in recent years, a
single stationary stochastic process model would not fit such data
appropriately. Hence, we rather split the dataset into two parts, representing
the relatively calm period from 01.10.2018–31.12.2020 (823 observations,
referred to as TS1) and the volatile period from 01.01.2021–01.01.2023 (731
observations, referred to as TS2). We remark that more sophisticated methods,
e.g. from change-point-detection, could be applied to split the dataset.
Figure 4 depicts the time series, the empirical autocorrelation functions and
the empirical densities of both time series. We note a rather slowly decaying
serial correlation with a distinct weakly periodic pattern. The marginal
distributions are rather different for the calm and the volatile regime: In
the calm regime, we observe the well-known fact of a slightly skewed
distribution, which has slightly heavier tails than the normal distribution
and includes negative prices. For the volatile regime, the distribution
appears to be skewed, multimodal, heavy-tailed and with significantly larger
empirical mean and variance than in the calm regime.

Figure 3: Time series plot of the day-ahead electricity baseload prices for
Germany and Luxembourg from 1st October 2018 to 1st January 2023 recorded in
EUR/MWh.

(a) TS1: Daily prices from 01.10.2018–31.12.2020 (b) TS2: Daily prices from
01.01.2021–01.01.2023

(c) Empirical ACF of TS1 (d) Empirical ACF of TS2

(e) Empirical density and scaled histogram of TS1  (f) Empirical density and
scaled histogram of TS2

Figure 4: Exploratory data analysis of the day-ahead electricity baseload
prices for Germany and Luxembourg from 1st October 2018 to 1st January 2023
recorded in EUR/MWh, which has been split into two times series (TS1 and TS2).

We will now focus on estimating the kernel function of a periodic trawl
process. First, we consider the relatively calm period from
01.10.2018–31.12.2020 and estimate an exponential trawl function with period
$\tau=7$. We set $\Delta=1$, representing one day. We obtain the following
estimates: $\widehat{\lambda}=0.055$ and
$(\widehat{c(0)},\widehat{c(1)},\ldots,\widehat{c(6)})=(1.000,0.705,0.496,0.451,0.445,0.490,0.642)$.
The empirical and fitted autocorrelation functions are depicted in Figure 5a.
We note that the exponential decay appears to be too fast to capture the
empirical autocorrelation function well. Hence we also fit a supGamma periodic
trawl function. We first ran a GMM estimation on the empirical autocorrelation
function to estimate $\alpha$ and then set $\alpha=1.2$. We then employ our
method of moment estimator and obtain $\widehat{H}=1.264$, indicating a long-
memory regime, and
$(\widehat{c(0)},\widehat{c(1)},\ldots,\widehat{c(6)})=(1.000,0.789,0.582,0.539,0.534,0.583,0.752)$.
We note that the model fit looks very good, see Figure 5b. Note that the
sensitivity to the particular choice of $\alpha$ appeared low, but this needs
to be investigated in more detail in the future.

We repeat the analysis on the more volatile time period from
01.01.2021–01.01.2023. For the exponential-periodic trawl process we obtain
the estimates $\widehat{\lambda}=0.032$ and
$(\widehat{c(0)},\widehat{c(1)},\ldots,\widehat{c(6)})=(1.000,0.689,0.473,0.421,0.406,0.437,0.558)$,
which results in a decay which is too fast, see Figure 5c. Again, the fit of
the supGamma-periodic trawl model appears better, where we set $\alpha=5.297$
(based on a GMM-estimation) and then find that $\widehat{H}=1.298$, which as
before falls into the long-memory setting, and
$(\widehat{c(0)},\widehat{c(1)},\ldots,\widehat{c(6)})=(1.000,0.958,0.901,0.887,0.888,0.912,0.957)$,
see Figure 5d.

(a) Exponential-periodic trawl for TS1  (b) SupGamma-periodic trawl for TS1

(c) Exponential-periodic trawl for TS2 (d) SupGamma-periodic trawl for TS2

Figure 5: Empirical and fitted exponential-periodic and supGamma periodic
trawl for the time series TS1 and TS2.

This small illustration demonstrates how easily the kernel function of a
periodic trawl process can be fitted to periodic autocorrelation functions.

In future work, it will be interesting to further extend and fine-tune the
methodology, in particular also the GMM methodology derived in Section 6.2 to
also estimate the parameters of the driving Lévy noise (under suitable
identifiability conditions) in practical settings. Here, a multistep approach,
where the parameters of the kernel functions are inferred first (as in this
article), followed by the parameters of the Lévy basis, might be a useful
direction to investigate further. We also remark that in cases when the period
is not known, preliminary simulation experiments not reported here indicate
that the period can be estimated using a smoothed empirical periodogram, which
is in line with existing literature.

It will also be worthwhile to explore extensions beyond the stationary
setting, by for instance considering a linear combination of an additive
seasonal function and a periodic trawl process to allow for non-stationarities
in the mean.

##  8 Conclusion and outlook

This article introduced an alternative definition of periodic trawl processes
compared to the one given in Barndorff-Nielsen et al. (2014), which appears to
have slightly higher analytical tractability. We derived some of the key
probabilistic properties of periodic trawl processes and have studied relevant
examples including both short- and long-memory settings. We showed that a
slice method can be used to simulate periodic trawl processes effectively.
Under suitable technical conditions, which currently only cover short-memory
scenarios, we have proved the asymptotic normality of the sample mean, sample
autocovariances, sample autocorrelations and the GMM-estimator. We have
implemented the proposed methodology in R and showed that the serial
correlation of periodic trawl processes describes the empirical one of
electricity spot prices very well. These are promising results, which suggest
that it will be worthwhile to further advance the theory and statistical
methodology for periodic trawl processes. For instance, for the GMM approach
to be applicable in practice, we need to specify suitable identifiability
conditions for fully parametric settings. It will also be interesting to
account for multiple periods, for instance, including weekly and yearly
periodicities, and tailor methods from spectral analysis to such settings.
Moreover, model selection tools are needed for choosing the appropriate Lévy
basis. For the particular application considered here, we note that our
framework allows for a smooth (Gaussian) noise terms as well as for jumps,
where the latter can model the positive and even more dominant negative jumps,
as in Veraart & Veraart (2014), but also the more volatile behaviour observed
in the most recent electricity prices. That said, in future research, we
should also investigate whether including a stochastic volatility term, either
via temporal or spatial stochastic scaling of the periodic trawl process is
needed to describe the data even better, which would bring us into the general
framework of ambit fields and processes, see Barndorff-Nielsen et al. (2018).

From a theoretical point of view, it will be interesting to develop an
asymptotic theory for periodic trawl processes which can allow for long memory
settings. We would then need to check how well such an asymptotic theory works
in finite samples to decide whether confidence bounds for the obtained
parameters can be constructed based on asymptotic guarantees or whether
bootstrap approaches might be more useful in practice.

##  Appendix A Appendix

The appendix contains the proofs of all the technical results presented in the
main paper, additional examples and a discussion of when the technical
assumptions needed in our main theorems hold for periodic trawl processes.

###  A.1 Proof of the second order properties

First, we derive the joint characteristic/cumulant function.

######  Proposition 8.

Let $t_{1}<t_{2}$ and $\theta_{1},\theta_{2}\in\mathbb{R}$. Then

$\displaystyle\mathrm{Log}(\mathbb{E}(\exp(i(\theta_{1},\theta_{2})(Y_{t_{1}},Y_{t_{2}})^{\top})))=\mathrm{Log}(\mathbb{E}(\exp(i\theta_{1}Y_{t_{1}}+i\theta_{2}Y_{t_{2}})))$
$\displaystyle=\int_{(-\infty,t_{1}]\times\mathbb{R}}C_{L^{\prime}}(\theta_{1}p(t_{1}-s)\mathbb{I}_{(0,g(t_{1}-s))}(x)+\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x))dxds$
$\displaystyle+\int_{(t_{1},t_{2}]\times\mathbb{R}}C_{L^{\prime}}(\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x))dxds,$

where $C_{L^{\prime}}$ denotes the cumulant function of the Lévy seed
$L^{\prime}$ associated with the Lévy basis $L$.

###### Proof.

Let $t_{1}<t_{2}$ and $\theta_{1},\theta_{2}\in\mathbb{R}$. Then the joint
characteristic function is given by

$\displaystyle\mathbb{E}(\exp(i(\theta_{1},\theta_{2})(Y_{t_{1}},Y_{t_{2}})^{\top}))=\mathbb{E}(\exp(i\theta_{1}Y_{t_{1}}+i\theta_{2}Y_{t_{2}}))$
$\displaystyle=\mathbb{E}\left[\exp\left(i\theta_{1}\int_{(-\infty,t_{1}]\times\mathbb{R}}p(t_{1}-s)\mathbb{I}_{(0,g(t_{1}-s))}(x)L(dx,ds)\right.\right.$
$\displaystyle\left.\left.+i\theta_{2}\int_{(-\infty,t_{1}]\times\mathbb{R}}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x)L(dx,ds)+i\theta_{2}\int_{(t_{1},t_{2}]\times\mathbb{R}}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x)L(dx,ds)\right)\right]$
$\displaystyle=\mathbb{E}\left[\exp\left(i\int_{(-\infty,t_{1}]\times\mathbb{R}}\\{\theta_{1}p(t_{1}-s)\mathbb{I}_{(0,g(t_{1}-s))}(x)+\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x)\\}L(dx,ds)\right.\right.$
$\displaystyle\left.\left.+i\theta_{2}\int_{(t_{1},t_{2}]\times\mathbb{R}}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x)L(dx,ds)\right)\right]$
$\displaystyle=\exp\left(\int_{(-\infty,t_{1}]\times\mathbb{R}}C(\theta_{1}p(t_{1}-s)\mathbb{I}_{(0,g(t_{1}-s))}(x)+\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x);L^{\prime})dxds\right)$
$\displaystyle\quad\cdot\exp\left(\int_{(t_{1},t_{2}]\times\mathbb{R}}C(\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x);L^{\prime})dxds\right).$

I.e.

$\displaystyle\log(\mathbb{E}(\exp(i\theta_{1}Y_{t_{1}}+i\theta_{2}Y_{t_{2}})))$
$\displaystyle=\int_{(-\infty,t_{1}]\times\mathbb{R}}C(\theta_{1}p(t_{1}-s)\mathbb{I}_{(0,g(t_{1}-s))}(x)+\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x);L^{\prime})dxds$
$\displaystyle+\int_{(t_{1},t_{2}]\times\mathbb{R}}C(\theta_{2}p(t_{2}-s)\mathbb{I}_{(0,g(t_{2}-s))}(x);L^{\prime})dxds.$

∎

We can now easily derive the second-order properties of the periodic trawl
process:

###### Proof of Proposition 5.

For $t,t_{1},t_{2}\in\mathbb{R},t_{1}<t_{2}$, we have

$\displaystyle\mathbb{E}(Y_{t})$
$\displaystyle=\mathbb{E}(L^{\prime})\int_{-\infty}^{t}p(t-s)g(t-s)ds=\mathbb{E}(L^{\prime})\int_{0}^{\infty}p(u)g(u)du,$
$\displaystyle\mathrm{Var}(Y_{t})$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}p^{2}(u)g(u)du,$
$\displaystyle\mathrm{Cov}(Y_{t_{1}},Y_{t_{2}})$
$\displaystyle=-\left.\frac{\partial^{2}}{\partial\theta_{1}\partial\theta_{2}}\log(\mathbb{E}(\exp(i\theta_{1}Y_{t_{1}}+i\theta_{2}Y_{t_{2}})))\right|_{\theta_{1}=\theta_{2}=0}$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{-\infty}^{t_{1}}p(t_{1}-s)p(t_{2}-s)\min(g(t_{1}-s),g(t_{2}-s))ds,$
$\displaystyle\mathrm{Cor}(Y_{t_{1}},Y_{t_{2}})$
$\displaystyle=\frac{\int_{-\infty}^{t_{1}}p(t_{1}-s)p(t_{2}-s)\min(g(t_{1}-s),g(t_{2}-s))ds}{\int_{0}^{\infty}p^{2}(u)g(u)du}.$

Recall that we assume that $g$ is monotonically decreasing, i.e. if $x\leq y$,
then $g(x)\geq g(y)$.

Since $t_{1}<t_{2}$, we have $t_{1}-s<t_{2}-s$ for $s<t_{1}$ and

$\displaystyle\min(g(t_{1}-s),g(t_{2}-s))=g(t_{2}-s),$

hence the above expressions simplify to

$\displaystyle\mathrm{Cov}(Y_{t_{1}},Y_{t_{2}})$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{-\infty}^{t_{1}}p(t_{1}-s)p(t_{2}-s)g(t_{2}-s)ds$
$\displaystyle=\mathrm{Var}(L^{\prime})\int_{0}^{\infty}p(u)p(t_{2}-t_{1}+u)g(t_{2}-t_{1}+u)du,$
$\displaystyle\mathrm{Cor}(Y_{t_{1}},Y_{t_{2}})$
$\displaystyle=\frac{\int_{0}^{\infty}p(u)p(t_{2}-t_{1}+u)g(t_{2}-t_{1}+u)du}{\int_{0}^{\infty}p^{2}(u)g(u)du}.$

∎

###### Proof of Proposition 6.

Recall that

$\displaystyle\mathrm{Cor}(Y_{0},Y_{t})=\frac{\int_{0}^{\infty}p(u)p(t+u)g(t+u)du}{\int_{0}^{\infty}p^{2}(u)g(u)du}.$

We consider a constant $M>\tau$. Since $p$ is periodic with period $\tau$,
there exist $\xi_{1},\xi_{2}\in[0,\tau]$ such that,

$\displaystyle\int_{0}^{M}p(u)p(t+u)g(t+u)du$
$\displaystyle=p(\xi_{1})p(\xi_{1}+t)\int_{0}^{M}g(t+u)du,$
$\displaystyle\int_{0}^{M}(p(u))^{2}g(u)du$
$\displaystyle=(p(\xi_{2}))^{2}\int_{0}^{M}g(t+u)du,$

by the mean value theorem. We note that $p$ is assumed to be continuous, and
since it is also periodic, it is bounded. Also, the integrability conditions
in (7) guarantee the existence of the integrals when taking the limit as
$M\to\infty$. Taking the limit and setting
$c(t)=p(\xi_{1})p(\xi_{1}+t)/(p(\xi_{2}))^{2}$ leads the result; since
$Cor(Y_{0},Y_{0})=1$, we deduce that $c(0)=1$. Also, we observe that $c$ is
proportional to the $\tau$-periodic function $p$ and is hence $\tau$-periodic
itself.

∎

######  Remark 11.

As mentioned in Remark 2, Barndorff-Nielsen et al. (2014) proposed adding a
periodic function as a multiplicative factor to $g$ rather than as kernel
function as in (8), which results in a process $(Z_{t})_{t\geq 0}$ with

$\displaystyle Z_{t}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}\mathbb{I}_{(0,p(t-s)g(t-s))}(x)\mathbb{I}_{[0,\infty)}(t-s)L(dx,ds),$
(19)

compared to our earlier definition of $(Y_{t})_{t\geq 0}$ with

$\displaystyle Y_{t}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}p(t-s)\mathbb{I}_{(0,g(t-s))}(x)\mathbb{I}_{[0,\infty)}(t-s)L(dx,ds).$

The autocorrelation function of the process $Z$ is of the form, for
$t_{1}<t_{2}$,

$\displaystyle\mathrm{Cor}(Z_{t_{1}},Z_{t_{2}})$
$\displaystyle=\frac{\int_{-\infty}^{t_{1}}\min(p(t_{1}-s)g(t_{1}-s),p(t_{2}-s)g(t_{2}-s))ds}{\int_{0}^{\infty}p(u)g(u)du},$

which is potentially slightly more difficult to deal with than the
autocorrelation function of our proposed periodic trawl process $Y$.

###  A.2 Proofs of the asymptotic theory

The following proofs extend the ideas presented in the work by Cohen & Lindner
(2013). Alternatively, we could have deduced the results from the more recent
work by Curato & Stelzer (2019).

###### Proof of Theorem 6.

The proof is a straightforward extension of the arguments provided in the
proof of Theorem 2.1 in Cohen & Lindner (2013). For the convenience of the
reader and to keep this article self-contained, we will present the steps to
extend the proof by Cohen & Lindner (2013) to our more general setting of
mixed moving average processes driven by homogeneous Lévy bases.

First of all, we continue the function $F_{\Delta}$ periodically on
$\mathbb{R}$ by setting

$\displaystyle F_{\Delta}(x,u)=\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)|,\quad
u\in\mathbb{R},$

where $F_{\Delta}(x,u)=F_{\Delta}(x,u+j\Delta)$ for all $j\in\mathbb{Z}$,
$u,x\in\mathbb{R}$.

We note that the autocovariance function of $Y$ satisfies

$\displaystyle|\gamma_{f}(j\Delta)|\leq\kappa_{2}\int_{\mathbb{R}\times\mathbb{R}}|f(x,-s)||f(x,j\Delta-s)|dxds,$

for any $j\in\mathbb{Z}$ and

$\displaystyle\begin{split}\frac{1}{\kappa_{2}}\sum_{j=-\infty}^{\infty}|\gamma_{f}(j\Delta)|&\leq\sum_{j=-\infty}^{\infty}\int_{\mathbb{R}\times\mathbb{R}}|f(x,-s)||f(x,j\Delta-s)|dxds\\\
&\leq\int_{\mathbb{R}\times\mathbb{R}}|f(x,-s)|\sum_{j=-\infty}^{\infty}|f(x,j\Delta-s)|dxds\\\
&=\int_{\mathbb{R}\times\mathbb{R}}|f(x,-s)|F_{\Delta}(x,-s)dxds\\\
&=\int_{\mathbb{R}\times\mathbb{R}}|f(x,s)|F_{\Delta}(x,s)dxds\\\
&=\sum_{j=-\infty}^{\infty}\int_{\mathbb{R}\times[0,\Delta]}|f(x,j\Delta+s)|F_{\Delta}(x,s)dxds\\\
&=\int_{\mathbb{R}\times[0,\Delta]}\sum_{j=-\infty}^{\infty}|f(x,j\Delta+s)|F_{\Delta}(x,s)dxds\\\
&=\int_{\mathbb{R}\times[0,\Delta]}F_{\Delta}^{2}(x,s)dxds<\infty.\end{split}$
(20)

The above computations can be repeated without the modulus, which implies that
$\sum_{j=-\infty}^{\infty}\gamma_{f}(j\Delta)=V_{\Delta}$.

To simplify the exposition, we shall now assume that $\mu=0$. We proceed as in
Cohen & Lindner (2013). Define the function
$f_{m;\Delta}(x,s):=f(x,s)\mathbb{I}_{(-m\Delta,m\Delta)}(s)$, for
$m\in\mathbb{N},x,s\in\mathbb{R}$, and set

$\displaystyle Y_{j;\Delta}^{m}$
$\displaystyle:=\int_{\mathbb{R}\times\mathbb{R}}f_{m;\Delta}(x,s)L(dx,ds)=\int_{\mathbb{R}\times(-m\Delta,m\Delta)}f(x,s)L(dx,ds)$
$\displaystyle=\int_{\mathbb{R}\times((-m+j)\Delta,(m+j)\Delta)}f(x,j\Delta-s)L(dx,ds).$

Since $L$ is independently scattered, we can deduce that
$(Y_{j;\Delta}^{(m)})_{j\in\mathbb{Z}}$ is a $(2m-1)$-dependent sequence,
which is also strictly stationary. Hence, by Brockwell & Davis (1987, Theorem
6.4.2), we know that

$\displaystyle\sqrt{n}\;\overline{Y}_{n;\Delta}^{(m)}=n^{-1/2}\sum_{j=1}^{n}Y_{j;\Delta}^{(m)}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}Z^{(m)}_{\Delta},\quad\mathrm{as}\;n\to\infty,$

where the random variable $Z^{(m)}_{\Delta}$ satisfies
$Z^{(m)}_{\Delta}\stackrel{{\scriptstyle\mathrm{d}}}{{=}}\mathrm{N}(0,V^{(m)})$,
where

$\displaystyle V^{(m)}_{\Delta}=\sum_{j=-2m}^{2m}\gamma_{f_{m}}(j\Delta),$

for

$\displaystyle\gamma_{f_{m}}(j\Delta)$
$\displaystyle=\mathrm{Cov}(Y_{0;\Delta}^{(m)},Y_{j;\Delta}^{(m)})=\kappa_{2}\int_{\mathbb{R}\times\mathbb{R}}f_{m;\Delta}(x,-s)f_{m;\Delta}(x,j\Delta-s)dxds$
$\displaystyle=\int_{\mathbb{R}\times((-m+j)\Delta,(m+j)\Delta)}f(x,-s)f(x,j\Delta-s)dxds.$

We observe that $\lim_{m\to\infty}\gamma_{f_{m}}(j\Delta)=\gamma_{f}(j\Delta)$
for all $j\in\mathbb{Z}$; also

$\displaystyle|\gamma_{f_{m}}(j\Delta)|\leq\kappa_{2}\int_{\mathbb{R}\times\mathbb{R}}|f(x,-s)||f(x,j\Delta-s)|dxds,$

and
$\sum_{j=-\infty}^{\infty}\int_{\mathbb{R}\times\mathbb{R}}|f(x,-s)||f(x,j\Delta-s)|<\infty$
by the computations in (20). Hence, Lebesgue’s Dominated Convergence Theorem
implies that $\lim_{m\to\infty}V^{(m)}_{\Delta}=V_{\Delta}$ and we get that
$Z_{\Delta}^{(m)}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}Z_{\Delta}$, where
$Z\stackrel{{\scriptstyle\mathrm{d}}}{{=}}\mathrm{N}(0,V_{\Delta})$.

It remains to control the difference
$n^{1/2}(\overline{Y}_{n;\Delta}-\overline{Y}_{n;\Delta}^{(m)})$. We argue as
follows. Using similar arguments as above, we note that
$\lim_{m\to\infty}\sum_{j=-\infty}^{\infty}\gamma_{f-f_{m;\Delta}}(j\Delta)=0$.
Hence, we have that

$\displaystyle\lim_{m\to\infty}\lim_{n\to\infty}\mathrm{Var}(n^{1/2}(\overline{Y}_{n;\Delta}-\overline{Y}_{n;\Delta}^{(m)}))$
$\displaystyle=\lim_{m\to\infty}\lim_{n\to\infty}n\mathrm{Var}\left(n^{-1}\sum_{j=1}^{n}\int_{\mathbb{R}\times\mathbb{R}}(f(x,j\Delta-s)-f_{m;\Delta}(x,j\Delta-s))L(dx,ds)\right)$
$\displaystyle\stackrel{{\scriptstyle(\star)}}{{=}}\lim_{m\to\infty}\sum_{j=-\infty}^{\infty}\gamma_{f-f_{m;\Delta}}(j\Delta)=0,$

where the equality $(\star)$ follows from Brockwell & Davis (1987, Theorem
7.1.1). Chebychef’s inequality allows us to conclude that, for any
$\epsilon>0$,

$\displaystyle\lim_{m\to\infty}\limsup_{n\to\infty}\mathbb{P}(n^{1/2}|\overline{Y}_{n;\Delta}-\overline{Y}_{n;\Delta}^{(m)}|>\epsilon)=0.$

As stated in Cohen & Lindner (2013), the final step of the proof consists of
an application of a Slutsky-type theorem as presented in Brockwell & Davis
(1987, Proposition 6.3.9). ∎

###### Proof of Lemma 1.

For $t_{1},t_{2},t_{3},t_{4}\in\mathbb{R}$, we have, for any
$a_{1},a_{2},a_{3},a_{4}\in\mathbb{R}$, the following expression for the joint
characteristic function

$\displaystyle\psi((a_{1},a_{2},a_{3},a_{4});(Y_{t_{1}},Y_{t_{2}},Y_{t_{3}},Y_{t_{4}})):=\mathbb{E}(\exp(i(a_{1}Y_{t_{1}}+a_{2}Y_{t_{2}}+a_{3}Y_{t_{3}}+a_{4}Y_{t_{4}}))$
$\displaystyle=\mathbb{E}\left[\exp\left(i\int_{\mathbb{R}\times\mathbb{R}}\left(\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)\right)L(ds,dx)\right)\right]$
$\displaystyle=\exp\left[\int_{\mathbb{R}\times\mathbb{R}}C\left(\sum_{j=1}^{4}a_{j}f(x,t_{j}-s);L^{\prime}\right)dxds\right],$

where $C(\cdot;L^{\prime})$ denotes the cumulant function of the Lévy seed
$L^{\prime}$, which we will present next.

Suppose $L^{\prime}$ has characteristic triplet $(c,A,\nu)$ w.r.t. the
truncation function $\tau(y)=\mathbb{I}_{[-1,1]}(y)$. I.e. we have the
following representation for its characteristic function, for any
$\theta\in\mathbb{R}$,

$\displaystyle\mathbb{E}(\exp(i\theta
L^{\prime}))=\exp\left(ic\theta-\frac{1}{2}A\theta^{2}+\int_{\mathbb{R}}(e^{iy\theta}-1-i\theta
y\tau(y))\nu(dy)\right).$

We recall that
$\mathbb{E}(L^{\prime})=c+\int_{\mathbb{R}}y(1-\tau(y))\nu(dy)$. Since we are
assuming that $\mathbb{E}(L^{\prime})=0$, we get that
$c=-\int_{\mathbb{R}}y(1-\tau(y))\nu(dy)=-\int_{\mathbb{R}}y\mathbb{I}_{[-1,1]^{c}}(y)\nu(dy)$
and, hence,

$\displaystyle\mathbb{E}(\exp(i\theta
L^{\prime}))=\exp\left(-\frac{1}{2}A\theta^{2}+\int_{\mathbb{R}}(e^{iy\theta}-1-i\theta
y)\nu(dy)\right).$

I.e. the corresponding cumulant function is given by

$\displaystyle
C(\theta;L^{\prime})=-\frac{1}{2}A\theta^{2}+\int_{\mathbb{R}}(e^{iy\theta}-1-i\theta
y)\nu(dy).$

Moreover,

$\displaystyle
C((a_{1},a_{2},a_{3},a_{4});(Y_{t_{1}},Y_{t_{2}},Y_{t_{3}},Y_{t_{4}})):=\int_{\mathbb{R}\times\mathbb{R}}C\left(\sum_{j=1}^{4}a_{j}f(x,t_{j}-s);L^{\prime}\right)dxds$
$\displaystyle=-\frac{1}{2}A\int_{\mathbb{R}\times\mathbb{R}}\left(\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)\right)^{2}dxds$
$\displaystyle\quad+\int_{\mathbb{R}\times\mathbb{R}}\int_{\mathbb{R}}\left(e^{iy\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)}-1-i\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)y\right)\nu(dy)dxds$
$\displaystyle=-\frac{1}{2}A\sum_{j,k=1}^{4}a_{j}a_{k}\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{j}-s)f(x,t_{k}-s)dxds$
$\displaystyle\quad+\int_{\mathbb{R}\times\mathbb{R}}\int_{\mathbb{R}}\left(e^{iy\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)}-1-i\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)y\right)\nu(dy)dxds$
$\displaystyle=-\frac{1}{2}A\sum_{j=1}^{4}a_{j}^{2}\int_{\mathbb{R}\times\mathbb{R}}f^{2}(x,t_{j}-s)dxds-\frac{1}{2}A\sum_{j,k=1,j\not=k}^{4}a_{j}a_{k}\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{j}-s)f(x,t_{k}-s)dxds$
$\displaystyle\quad+\int_{\mathbb{R}\times\mathbb{R}}\int_{\mathbb{R}}\left(e^{iy\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)}-1-i\sum_{j=1}^{4}a_{j}f(x,t_{j}-s)y\right)\nu(dy)dxds.$

Next, we compute the fourth moments, where we recall that

$\displaystyle\left.\frac{\partial^{4}}{\partial a_{1}\partial a_{2}\partial
a_{3}\partial
a_{4}}\psi((a_{1},a_{2},a_{3},a_{4});(Y_{t_{1}},Y_{t_{2}},Y_{t_{3}},Y_{t_{4}}))\right|_{a_{1}=a_{2}=a_{3}=a_{4}=0}=\mathbb{E}(Y_{t_{1}}Y_{t_{2}}Y_{t_{3}}Y_{t_{4}}).$

We now abbreviate the functions to $\psi$ and $C$ without stating their
arguments and a subscript denotes the corresponding partial derivative, e.g.
$C_{a_{1}}=\frac{\partial}{\partial
a_{1}}C((a_{1},a_{2},a_{3},a_{4});(Y_{t_{1}},Y_{t_{2}},Y_{t_{3}},Y_{t_{4}})$
and similarly for higher order partial derivatives. Since $\psi=\exp(C)$, we
have

$\displaystyle\psi_{a1}$ $\displaystyle=\psi C_{a_{1}},$
$\displaystyle\psi_{a_{1},a_{2}}$
$\displaystyle=\psi[C_{a_{1},a_{2}}+C_{a_{1}}C_{a_{2}}],$
$\displaystyle\psi_{a_{1},a_{2},a_{3}}$
$\displaystyle=\psi[(C_{a_{1},a_{2}}+C_{a_{1}}C_{a_{2}})C_{a_{3}}+C_{a_{1},a_{2},a_{3}}+C_{a_{1},a_{3}}C_{a_{2}}+C_{a_{1}}C_{a_{2},a_{3}}]$
$\displaystyle=\psi[C_{a_{1}}C_{a_{2}}C_{a_{3}}+C_{a_{1}}C_{a_{2},a_{3}}+C_{a_{2}}C_{a_{1},a_{3}}+C_{a_{3}}C_{a_{1},a_{2}}++C_{a_{1},a_{2},a_{3}}],$
$\displaystyle\psi_{a_{1},a_{2},a_{3},a_{4}}$
$\displaystyle=\psi[C_{a_{1}}C_{a_{2}}C_{a_{3}}C_{a_{4}}+C_{a_{1}}C_{a_{2},a_{3}}C_{a_{4}}+C_{a_{2}}C_{a_{1},a_{3}}C_{a_{4}}+C_{a_{3}}C_{a_{1},a_{2}}C_{a_{4}}+C_{a_{1},a_{2},a_{3}}C_{a_{4}}$
$\displaystyle+C_{a_{1},a_{4}}C_{a_{2}}C_{a_{3}}+C_{a_{1}}C_{a_{2},a_{4}}C_{a_{3}}+C_{a_{1}}C_{a_{2}}C_{a_{3},a_{4}}+C_{a_{1},a_{4}}C_{a_{2},a_{3}}+C_{a_{1}}C_{a_{2},a_{3},a_{4}}$
$\displaystyle+C_{a_{2},a_{4}}C_{a_{1},a_{3}}+C_{a_{2}}C_{a_{1},a_{3},a_{4}}+C_{a_{3},a_{4}}C_{a_{1},a_{2}}+C_{a_{3}}C_{a_{1},a_{2},a_{4}}+C_{a_{1},a_{2},a_{3},a_{4}}].$

Here we have

$\displaystyle C_{a_{i}}|_{a_{1}=a_{2}=a_{3}=a_{4}=0}$
$\displaystyle=0,\quad\mathrm{for}\;i=1,2,3,4,$ $\displaystyle
C_{a_{i},a_{j}}|_{a_{1}=a_{2}=a_{3}=a_{4}=0}$
$\displaystyle=-\left(A+\int_{\mathbb{R}}y^{2}\nu(dy)\right)\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{i}-s)f(x,t_{j}-s)dxds,\quad\mathrm{for}\,i,j=1,2,3,4,$
$\displaystyle C_{a_{1},a_{2},a_{3},a_{4}}$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}\prod_{j=1}^{4}f(x,t_{j}-s)dxds\int_{\mathbb{R}}y^{4}\nu(dy).$

The above results imply that

$\displaystyle\mathbb{E}(Y_{t_{1}}Y_{t_{2}}Y_{t_{3}}Y_{t_{4}})$
$\displaystyle=\left(A+\int_{\mathbb{R}}y^{2}\nu(dy)\right)^{2}$
$\displaystyle\left(\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-s)f(x,t_{2}-s)dxds\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{3}-s)f(x,t_{4}-s)dxds\right.$
$\displaystyle+\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-s)f(x,t_{3}-s)dxds\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{2}-s)f(x,t_{4}-s)dxds$
$\displaystyle\left.+\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-s)f(x,t_{4}-s)dxds\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{2}-s)f(x,t_{3}-s)dxds\right)$
$\displaystyle+\int_{\mathbb{R}}y^{4}\nu(dy)\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-s)f(x,t_{2}-s)f(x,t_{3}-s)f(x,t_{4}-s)dxds.$

We note that
$\kappa_{4}:=\int_{\mathbb{R}}y^{4}\nu(dy)=(\eta-3)\kappa_{2}^{2}$ and
$\kappa_{2}=A+\int_{\mathbb{R}}y^{2}\nu(dy)$

We can further simplify the above formula as follows:

$\displaystyle\mathbb{E}(Y_{t_{1}}Y_{t_{2}}Y_{t_{3}}Y_{t_{4}})$
$\displaystyle=\left(A+\int_{\mathbb{R}}y^{2}\nu(dy)\right)^{2}$
$\displaystyle\left(\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-t_{2}+s)f(x,s)dxds\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{3}-t_{4}+s)f(x,s)dxds\right.$
$\displaystyle+\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-t_{3}+s)f(x,s)dxds\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{2}-t_{4}+s)f(x,s)dxds$
$\displaystyle\left.+\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-t_{4}+s)f(x,s)dxds\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{2}-t_{3}+s)f(x,s)dxds\right)$
$\displaystyle+\int_{\mathbb{R}}y^{4}\nu(dy)\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-t_{3}+s)f(x,t_{2}-t_{3})f(x,s)f(x,t_{4}-t_{3}+s)dxds$
$\displaystyle=\gamma(t_{1}-t_{2})\gamma(t_{3}-t_{4})+\gamma(t_{1}-t_{3})\gamma(t_{2}-t_{4})+\gamma(t_{1}-t_{4})\gamma(t_{2}-t_{3})$
$\displaystyle+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,t_{1}-t_{3}+s)f(x,t_{2}-t_{3})f(x,s)f(x,t_{4}-t_{3}+s)dxds.$

∎

###### Proof of Proposition 7.

We first expand the covariance of the sample autocovariances as follows

$\displaystyle\mathrm{Cov}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p),\widehat{\gamma}_{n,\Delta}^{*}(\Delta
q))=\mathbb{E}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p)\widehat{\gamma}_{n;\Delta}^{*}(\Delta
q))-\mathbb{E}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p))\mathbb{E}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta q)),$

where

$\displaystyle\mathbb{E}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta p))$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n}\mathbb{E}(Y_{j\Delta}Y_{(j+p)\Delta})=\frac{1}{n}\sum_{j=1}^{n}\gamma(p\Delta)=\gamma(p\Delta),$
$\displaystyle\mathbb{E}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta q))$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{n}\mathbb{E}(Y_{k\Delta}Y_{(k+q)\Delta})=\frac{1}{n}\sum_{k=1}^{n}\gamma(q\Delta)=\gamma(q\Delta).$

Also,

$\displaystyle\mathbb{E}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p)\widehat{\gamma}_{n;\Delta}^{*}(\Delta
q))=\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}\mathbb{E}(Y_{j\Delta}Y_{(j+p)\Delta}Y_{k\Delta}Y_{(k+q)\Delta}),$

where

$\displaystyle\mathbb{E}(Y_{j\Delta}Y_{(j+p)\Delta}Y_{k\Delta}Y_{(k+q)\Delta})$
$\displaystyle=\gamma(p\Delta)\gamma(q\Delta)+\gamma((k-j)\Delta)\gamma((j+p-k-q)\Delta)+\gamma((j-k-q)\Delta)\gamma((j+p-k)\Delta)$
$\displaystyle\quad+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,(j-k)\Delta+s)f(x,(j-k+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds$
$\displaystyle=\gamma(p\Delta)\gamma(q\Delta)+\gamma((j-k)\Delta)\gamma((j-k+p-q)\Delta)+\gamma((j-k-q)\Delta)\gamma((j-k+p)\Delta)$
$\displaystyle\quad+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,(j-k)\Delta+s)f(x,(j-k+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds$
$\displaystyle\stackrel{{\scriptstyle
l=j-k}}{{=}}\gamma(p\Delta)\gamma(q\Delta)+\gamma(l\Delta)\gamma((l+p-q)\Delta)+\gamma((l-q)\Delta)\gamma((l+p)\Delta)$
$\displaystyle\quad+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,l\Delta+s)f(x,(l+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds.$

Now we subtract $\gamma(p\Delta)\gamma(q\Delta)$, we set $l=j-k$, interchange
the order of summation and use the stationarity to obtain

$\displaystyle\mathrm{Cov}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p),\widehat{\gamma}_{n;\Delta}^{*}(\Delta q))$
$\displaystyle=\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}\mathbb{E}(Y_{j\Delta}Y_{(j+p)\Delta}Y_{k\Delta}Y_{(k+q)\Delta})-\gamma(p\Delta)\gamma(q\Delta)$
$\displaystyle=\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}\left(\gamma((j-k)\Delta)\gamma((j-k+p-q)\Delta)+\gamma((j-k-q)\Delta)\gamma((j-k+p)\Delta)\right.$
$\displaystyle\quad+\left.\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,(j-k)\Delta+s)f(x,(j-k+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds\right)$
$\displaystyle=\frac{1}{n^{2}}\sum_{|l|<n}\sum_{k=1}^{n-|l|}\left(\gamma(l\Delta)\gamma((l+p-q)\Delta)+\gamma((l-q)\Delta)\gamma((l+p)\Delta)\right.$
$\displaystyle\quad+\left.\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,l\Delta+s)f(x,(l+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds\right)$
$\displaystyle=\frac{1}{n^{2}}\sum_{|l|<n}(n-|l|)T_{l,p,q;\Delta}=\frac{1}{n}\sum_{|l|<n}\left(1-\frac{|l|}{n}\right)T_{l,p,q;\Delta},$

where

$\displaystyle T_{l,p,q;\Delta}$
$\displaystyle:=\gamma(l\Delta)\gamma((l+p-q)\Delta)+\gamma((l-q)\Delta)\gamma((l+p)\Delta)$
$\displaystyle\quad+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}f(x,l\Delta+s)f(x,(l+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds.$

Hence, we have

$\displaystyle\lim_{n\to\infty}n\mathrm{Cov}(\widehat{\gamma}_{n;\Delta}^{*}(\Delta
p),\widehat{\gamma}_{n;\Delta}^{*}(\Delta
q))=\lim_{n\to\infty}\sum_{|l|<n}\left(1-\frac{|l|}{n}\right)T_{l,p,q;\Delta}=\sum_{l=-\infty}^{\infty}T_{l,p,q;\Delta},$

where

$\displaystyle\sum_{l=-\infty}^{\infty}T_{l,p,q;\Delta}$
$\displaystyle=\sum_{l=-\infty}^{\infty}[\gamma(l\Delta)\gamma((l+p-q)\Delta)+\gamma((l-q)\Delta)\gamma((l+p)\Delta)]$
$\displaystyle\quad+\kappa_{4}\int_{\mathbb{R}\times\mathbb{R}}\sum_{l=-\infty}^{\infty}f(x,l\Delta+s)f(x,(l+p)\Delta+s)f(x,s)f(x,q\Delta+s)dxds,$

by the Dominated Convergence Theorem since (13) holds. More precisely, let us
justify why $\sum_{l=-\infty}^{\infty}|T_{l,p,q;\Delta}|<\infty$. The
finiteness of
$\sum_{l=-\infty}^{\infty}|\gamma(l\Delta)\gamma((l+p-q)\Delta)+\gamma((l-q)\Delta)\gamma((l+p)\Delta)|$
follows from (13). For the second term, for $q\in\mathbb{Z}$, define

$\displaystyle\left(\widetilde{G}_{q;\Delta}:\mathbb{R}\times[0,\Delta]\to\mathbb{R},(x,u)\mapsto\widetilde{G}_{q;\Delta}(x,u)=\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)||f(x,u+(j+q)\Delta)|\right),$

which is in $L^{2}(\mathbb{R}\times[0,\Delta])$ due to (11). We consider the
periodic continuation of $\widetilde{G}_{q;\Delta}$ and set

$\displaystyle\left(\widetilde{G}_{q;\Delta}:\mathbb{R}\times\mathbb{R}\to\mathbb{R},(x,u)\mapsto\widetilde{G}_{q;\Delta}(x,u)=\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)||f(x,u+(j+q)\Delta)|\right).$

Since $\widetilde{G}_{q;\Delta}$ is periodic and, restricted to
$\mathbb{R}\times[0,\Delta]$ square-integrable, we have

$\displaystyle\int_{\mathbb{R}\times\mathbb{R}}\sum_{l=-\infty}^{\infty}|f(x,l\Delta+s)f(x,(l+p)\Delta+s)||f(x,s)f(x,q\Delta+s)|dxds$
$\displaystyle=\int_{\mathbb{R}\times\mathbb{R}}\widetilde{G}_{p}(s)|f(x,s)f(x,q\Delta+s)|dxds=\sum_{j=-\infty}^{\infty}\int_{\mathbb{R}\times[j\Delta,(j+1)\Delta]}\widetilde{G}_{p}(s)|f(x,s)f(x,q\Delta+s)|dxds$
$\displaystyle=\sum_{j=-\infty}^{\infty}\int_{\mathbb{R}\times[0,\Delta]}\widetilde{G}_{p}(s+j\Delta)|f(x,s+j\Delta)f(x,(q+j)\Delta+s)|dxds$
$\displaystyle=\sum_{j=-\infty}^{\infty}\int_{\mathbb{R}\times[0,\Delta]}\widetilde{G}_{p}(s)|f(x,s+j\Delta)f(x,(q+j)\Delta+s)|dxds$
$\displaystyle=\int_{\mathbb{R}\times[0,\Delta]}\widetilde{G}_{p}(s)\sum_{j=-\infty}^{\infty}|f(x,s+j\Delta)f(x,(q+j)\Delta+s)|dxds$
$\displaystyle=\int_{\mathbb{R}\times[0,\Delta]}\widetilde{G}_{p}(s)\widetilde{G}_{q}(s)dxds<\infty.$

Equation (LABEL:eq:sumT) follows from the same calculations as above without
the modulus sign in the definition of $\widetilde{G}$.

∎

###### Proof of Theorem 7.

  1. 1.

For a function $f$ with compact support, the result can be deduced as in
Brockwell & Davis (1987, Proposition 7.3.2). The general case can be handled
as follows, where we adapt the proof of Cohen & Lindner (2013, Theorem 3.5) to
our more general setting. As in the proof for the sample mean, define the
function $f_{m;\Delta}(x,s):=f(x,s)\mathbb{I}_{(-m\Delta,m\Delta)}(s)$, for
$m\in\mathbb{N},x,s\in\mathbb{R}$, and set

$\displaystyle Y_{j;\Delta}^{m}$
$\displaystyle:=\int_{\mathbb{R}\times\mathbb{R}}f_{m;\Delta}(x,s)L(dx,ds)=\int_{\mathbb{R}\times((-m+j)\Delta,(m+j)\Delta)}f(x,j\Delta-s)L(dx,ds).$

We denote by $\gamma_{m}$ the autocovariance function of the process
$(Y_{j;\Delta}^{m})_{j\in\mathbb{Z}}$. We set

$\displaystyle\gamma_{n;\Delta;m}^{*}(p\Delta)=\sum_{j=1}^{n}Y_{j;\Delta}^{m}Y_{j+p;\Delta}^{m},\quad
p=0,\ldots,h.$

Then, we have

$\displaystyle\sqrt{n}(\gamma_{n;\Delta;m}^{*}(0)-\gamma_{m}(0),\ldots,\gamma_{n;\Delta;m}^{*}(h\Delta)-\gamma_{m}(h\Delta))^{\top}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}Z_{\Delta;m}\sim\mathrm{N}(0,V_{\Delta;m}),\quad
n\to\infty,$

where the asymptotic covariance matrix is given by
$V_{\Delta;m}=(v_{pq;\Delta;m})_{p,q=0,\ldots,h}\in\mathbb{R}^{h+1,h+1}$ with
$v_{pq;\Delta;m}$ defined as

$\displaystyle v_{pq;\Delta;m}$
$\displaystyle:=(\eta-3)\kappa_{2}^{2}\int_{\mathbb{R}\times[0,\Delta]}G_{p;\Delta;m}(x,u)G_{q;\Delta;m}(x,u)dxdu$
$\displaystyle+\sum_{l=-\infty}^{\infty}[\gamma_{m}(l\Delta)\gamma_{m}((l+p-q)\Delta)+\gamma_{m}((l-q)\Delta)\gamma_{m}((l+p)\Delta),]$
$\displaystyle G_{q;\Delta;m}(x,u)$
$\displaystyle:=\sum_{j=-\infty}^{\infty}f_{m}(x,u+j\Delta)f_{m}(x,u+(j+q)\Delta);\quad
u\in[0,\Delta].$

We would like to show that $\lim_{m\to\infty}V_{\Delta;m}=V_{\Delta}$. For
this, we note that

$\displaystyle G_{q;\Delta;m}(x,u)$
$\displaystyle=\sum_{j=-\infty}^{\infty}f_{m}(x,u+j\Delta)f_{m}(x,u+(j+q)\Delta)$
$\displaystyle\to
G_{q;\Delta}(x,u)=\sum_{j=-\infty}^{\infty}f(x,u+j\Delta)f(x,u+(j+q)\Delta),$

uniformly in $u\in[0,\Delta]$, as $m\to\infty$, by Lebesgue’s Dominated
Convergence Theorem, since the function
$(x,u)\mapsto\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)||f(x,u+(j+q)\Delta)|$ is
in $L^{2}(\mathbb{R}\times[0,\Delta])$ by (11) and is therefore almost surely
finite. Moreover, we note that

$\displaystyle|G_{q;\Delta;m}(x,u)|$
$\displaystyle\leq\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)||f(x,u+(j+q)\Delta)|,$

uniformly in $u$ and $m$. Hence, an application of the Dominated Convergence
Theorem leads to $G_{q;\Delta;m}\to G_{q;\Delta}$ in
$L^{2}(\mathbb{R}\times[0,\Delta])$ as $m\to\infty$. We also note that

$\displaystyle|\gamma_{m}(j\Delta)|$
$\displaystyle\leq\int_{\mathbb{R}\times\mathbb{R}}|f(x,u)||f(x,u+j\Delta)|dxdu,\quad\forall
m\in\mathbb{N},\forall j\in\mathbb{Z}.$

Moreover, $\lim_{m\to\infty}\gamma_{m}(j\Delta)=\gamma(j\Delta)$ for all
$j\in\mathbb{Z}$. Assumption (15) together with the Dominated Convergence
Theorem allows us to conclude that $(\gamma_{m}(j\Delta))_{j\in\mathbb{Z}}$
converges in $L^{2}(\mathbb{Z})$ to $(\gamma(j\Delta))_{j\in\mathbb{Z}}$.
Combining this result with our earlier finding of the convergence of
$G_{q;\Delta;m}$ implies that $\lim_{m\to\infty}V_{\Delta;m}=V_{\Delta}$ and

$\displaystyle
Z_{\Delta;m}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}Z_{\Delta},\quad
m\to\infty,$

where
$Z_{\Delta}\stackrel{{\scriptstyle\mathrm{d}}}{{=}}\mathrm{N}(0,V_{\Delta})$.

Now, using the same arguments as in the proof of Brockwell & Davis (1987,
Equation (7.3.9)), we can show that

$\displaystyle\lim_{m\to\infty}\limsup_{n\to\infty}\mathbb{P}(n^{1/2}|\gamma_{n;\Delta;m}^{*}(q\Delta)-\gamma_{m}(q\Delta)-\gamma^{*}(q\Delta)+\gamma(q\Delta)|>\epsilon)=0,$

for all $\epsilon>0,q\in\\{0,\ldots,h\\}$. An application of a variant of
Slutsky’s theorem, see Brockwell & Davis (1987, Proposition 6.3.9) completes
the proof.

  2. 2.

This part can be proven in similar way as the proof of Brockwell & Davis
(1987, Proposition 7.3.4). Also, as in the proof of Cohen & Lindner (2013,
Theorem 3.5 b)), we observe that $\sqrt{n}\overline{Y}_{n;\Delta}$ converges
to a Gaussian random variable as $n\to\infty$ due to Theorem 6 and
$\overline{Y}_{n;\Delta}$ converges to 0 in probability as $n\to\infty$ (since
we assume here that $\mu=0).$

  3. 3.

For the final part of the theorem, we can argue as in the proof of Brockwell &
Davis (1987, Theorem 7.2.1), where the $w_{pq;\Delta}$ are obtained via the
Bartlett formula

$\displaystyle
w_{pq;\Delta}=(v_{pq;\Delta}-\rho(p\Delta)v_{0q;\Delta}-\rho(q\Delta)v_{p0;\Delta}+\rho(p\Delta)\rho(q\Delta)v_{00;\Delta})/\gamma^{2}(0).$

We can simplify the above formula and write

$\displaystyle w_{pq;\Delta}=w_{pq;\Delta}^{(1)}+w_{pq;\Delta}^{(2)},$

where

$\displaystyle w_{pq;\Delta}^{(1)}$
$\displaystyle:=\frac{(\eta-3)\kappa_{2}^{2}}{\gamma^{2}(0)}\int_{\mathbb{R}\times[0,\Delta]}(G_{p;\Delta}(x,u)-G_{0;\Delta}(x,u)\rho(p\Delta))$
$\displaystyle\quad\cdot(G_{q;\Delta}(x,u)-G_{0;\Delta}(x,u)\rho(q\Delta))dxdu,$

and

$\displaystyle w_{pq;\Delta}^{(2)}$
$\displaystyle:=\sum_{l=-\infty}^{\infty}\left[\rho(l\Delta)\rho((l+p-q)\Delta)+\rho((l-q)\Delta)\rho((l+p)\Delta)\right.$
$\displaystyle-2\rho(l\Delta)\rho((l-q)\Delta)\rho(p\Delta)-2\rho(l\Delta)\rho((l+p)\Delta)\rho(q\Delta)$
$\displaystyle\left.+2\rho(p\Delta)\rho(q\Delta)\rho^{2}(l\Delta)\right].$

Note that

$\displaystyle\sum_{l=-\infty}^{\infty}\rho(l\Delta)\rho((l+p-q)\Delta)$
$\displaystyle=\sum_{l=-\infty}^{\infty}\rho((l+q)\Delta)\rho((l+p)\Delta),$
$\displaystyle\sum_{l=-\infty}^{\infty}\rho(l\Delta)\rho((l-q)\Delta)\rho(p\Delta)$
$\displaystyle=\sum_{l=-\infty}^{\infty}\rho((l+q)\Delta)\rho(l\Delta)\rho(p\Delta).$

Hence,

$\displaystyle w_{pq;\Delta}^{(2)}$
$\displaystyle=\sum_{l=-\infty}^{\infty}\left[\rho((l+q)\Delta)\rho((l+p)\Delta)+\rho((l-q)\Delta)\rho((l+p)\Delta)\right.$
$\displaystyle-2\rho((l+q)\Delta)\rho(l\Delta)\rho(p\Delta)-2\rho(l\Delta)\rho((l+p)\Delta)\rho(q\Delta)$
$\displaystyle\left.+2\rho(p\Delta)\rho(q\Delta)\rho^{2}(l\Delta)\right].$

Hence,

$\displaystyle w_{pq;\Delta}$
$\displaystyle=\frac{(\eta-3)\kappa_{2}^{2}}{\gamma^{2}(0)}\int_{\mathbb{R}\times[0,\Delta]}(G_{p;\Delta}(x,u)-G_{0;\Delta}(x,u)\rho(p\Delta))$
$\displaystyle\quad\cdot(G_{q;\Delta}(x,u)-G_{0;\Delta}(x,u)\rho(q\Delta))dxdu$
$\displaystyle+\sum_{l=-\infty}^{\infty}\left[\rho((l+q)\Delta)\rho((l+p)\Delta)+\rho((l-q)\Delta)\rho((l+p)\Delta)\right.$
$\displaystyle-2\rho((l+q)\Delta)\rho(l\Delta)\rho(p\Delta)-2\rho(l\Delta)\rho((l+p)\Delta)\rho(q\Delta)$
$\displaystyle\left.+2\rho(p\Delta)\rho(q\Delta)\rho^{2}(l\Delta)\right].$

∎

###  A.3 Examples

In this subsection, we show how the asymptotic variances appearing in the
asymptotic theory for the sample mean can be computed for trawl processes with
either an exponential or a supGamma trawl function.

Note that, in the case when $p\equiv 1$, i.e. for the (non-periodic) trawl
process, we get

$\displaystyle\gamma(h)=\mathrm{Cov}(Y_{t},Y_{t+h})=\int_{|h|}^{\infty}g(u)du,$

for $t,h\in\mathbb{R}$.

####  A.3.1 Exponential trawl

Consider the case of an exponential trawl function with $g(x)=\exp(-\lambda
x)$, for $x\geq 0$. The autocovariance function is given by

$\displaystyle\gamma(t)=\mathrm{Cov}(Y_{t},Y_{0})=\frac{1}{\lambda}e^{-\lambda|t|},$

for $t\in\mathbb{R}$, and the autocorrelation function is given by

$\displaystyle\rho(t)=\mathrm{Cor}(Y_{t},Y_{0})=e^{-\lambda|t|},$

for $t\in\mathbb{R}$.

For the sample mean, we have the following result. Suppose that
$\mathbb{E}(L^{\prime})=0,\kappa_{2}=\mathrm{Var}(L^{\prime})<\infty,\mu\in\mathbb{R}$
and $\Delta>0$. Then

$\displaystyle\left(F_{\Delta}:\mathbb{R}\times[0,\Delta]\to[0,\infty],(x,u)\mapsto
F_{\Delta}(x,u)=\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)|\right)\in
L^{2}(\mathbb{R}\times[0,\Delta])$

since

$\displaystyle
F_{\Delta}(x,u)=\sum_{j=-\infty}^{\infty}|f(x,u+j\Delta)|=\sum_{j=-\infty}^{\infty}\mathbb{I}_{(0,g(u+j\Delta)}(x),$

and we have that

$\displaystyle\int_{\mathbb{R}\times[0,\Delta]}|F_{\Delta}(x,u)|^{2}dxdu$
$\displaystyle=\int_{\mathbb{R}\times[0,\Delta]}\sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty}\mathbb{I}_{(0,g(u+j\Delta)}(x)\mathbb{I}_{(0,g(u+k\Delta)}(x)dxdu$
$\displaystyle=\sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty}\int_{\mathbb{R}\times[0,\Delta]}\mathbb{I}_{(0,g(u+j\Delta)}(x)\mathbb{I}_{(0,g(u+k\Delta)}(x)dxdu$
$\displaystyle=\sum_{j=0}^{\infty}\sum_{k=0}^{\infty}\int_{[0,\Delta]}\min\\{g(u+j\Delta),g(u+k\Delta)\\}du$
$\displaystyle=\sum_{j=0}^{\infty}\int_{[0,\Delta]}g(u+j\Delta)du+2\sum_{j=0}^{\infty}\sum_{k=j+1}^{\infty}\int_{[0,\Delta]}g(u+k\Delta)du$

where we applied Tonelli’s theorem.

Then $\sum_{j=-\infty}^{\infty}|\gamma(\Delta j)|<\infty$,

$\displaystyle V_{\Delta}:=\sum_{j=-\infty}^{\infty}\gamma(\Delta
j)=\kappa_{2}\int_{\mathbb{R}\times[0,\Delta]}\left(\sum_{j=1}^{\infty}f(x,u+j\Delta)\right)^{2}dxdu,$
(21)

and the sample mean of $Y_{\Delta i}$, for $i=1,\ldots,n$, is asymptotically
Gaussian as $n\to\infty$, i.e.

$\displaystyle\sqrt{n}\left(\overline{Y}_{n;\Delta}-\mu\right)\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\mathrm{N}\left(0,V_{\Delta}\right),\quad\mathrm{as}\,n\to\infty.$

For the case of an exponential trawl function, we get

$\displaystyle V_{\Delta}$
$\displaystyle=\sum_{j=-\infty}^{\infty}\gamma(\Delta
j)=\sum_{j=-\infty}^{\infty}\frac{1}{\lambda}e^{-\lambda|j|\Delta}=\frac{1}{\lambda}\left(1+2\sum_{j=1}^{\infty}e^{-\lambda
j\Delta}\right)=\frac{1}{\lambda}\left(1+2\frac{e^{-\lambda\Delta}}{1-e^{-\lambda\Delta}}\right)$
$\displaystyle=\frac{1+e^{-\lambda\Delta}}{\lambda(1-e^{-\lambda\Delta})}=\frac{1+e^{\lambda\Delta}}{\lambda(e^{\lambda\Delta}-1)}.$

####  A.3.2 SupGamma trawl

In the case when $g(x)=(1+x/\alpha)^{-H}$, $\alpha>0,H>2$, i.e. we require a
short-memory setting, $x\geq 0$, we have

$\displaystyle\gamma(h)=\int_{|h|}^{\infty}g(x)dx=\frac{\alpha}{H-1}\left(1+\frac{|h|}{\alpha}\right)^{1-H}.$

Then

$\displaystyle V_{\Delta}$
$\displaystyle=\sum_{j=-\infty}^{\infty}\gamma(j\Delta)=\frac{\alpha}{H-1}\sum_{j=-\infty}^{\infty}\left(1+\frac{|j|\Delta}{\alpha}\right)^{1-H}$
$\displaystyle=\frac{\alpha}{H-1}\left(\frac{\Delta}{\alpha}\right)^{1-H}\sum_{j=-\infty}^{\infty}\left(\frac{\alpha}{\Delta}+|j|\right)^{1-H}$
$\displaystyle=\frac{\alpha}{H-1}\left(\frac{\Delta}{\alpha}\right)^{1-H}\left(\frac{\alpha}{\Delta}\right)^{1-H}\left[2\left(\frac{\alpha}{\Delta}\right)^{H-1}\zeta(H-1,\alpha/\Delta)-1\right]$
$\displaystyle=\frac{\alpha}{H-1}\left[2\left(\frac{\alpha}{\Delta}\right)^{H-1}\zeta(H-1,\alpha/\Delta)-1\right],$

where $\zeta$ denotes the Hurwitz Zeta function defined by
$\zeta(s,a)=\sum_{k=0}^{\infty}\frac{1}{(k+a)^{s}}$, for $\mathrm{Re}(s)>1$.

###  A.4 Verifying the assumptions of Theorem 7 for selected periodic trawl
processes

For the applications discussed in Section 5, we need to verify the condition
(11) from Proposition 7 and Assumption (15) from Theorem 7 assuming that the
corresponding moment assumptions for the Lévy seed hold.

For both conditions, it is sufficient to check that a (non-periodic) trawl
process satisfies the stated conditions since the periodic function is
bounded. Hence, in the following, we shall set $p\equiv 1$.

####  A.4.1 Verifying Condition (11) from Proposition 7

We need to check that

$\displaystyle\left(\mathbb{R}\times[0,\Delta]\to\mathbb{R},(x,u)\mapsto\sum_{j=-\infty}^{\infty}f^{2}(x,u+j\Delta)\right)\in
L^{2}(\mathbb{R}\times[0,\Delta]).$

This condition holds for trawl processes if, for
$(x,u)\in\mathbb{R}\times[0,\Delta]$,

$\displaystyle\sum_{j=-\infty}^{\infty}f^{2}(x,u+j\Delta)=\sum_{j=-\infty}^{\infty}f(x,u+j\Delta)=\sum_{j=-\infty}^{\infty}\mathbb{I}_{(0,g(u+j\Delta))}(x)\mathbb{I}_{[0,\infty)}(u+j\Delta)=F_{\Delta}(x,u)\in
L^{2}(\mathbb{R}\times[0,\Delta]).$

This is equivalent to checking that

$\displaystyle\begin{split}&\int_{\mathbb{R}\times[0,\Delta]}|F_{\Delta}(x,u)|^{2}dxdu\\\
&=\int_{\mathbb{R}\times[0,\Delta]}\sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty}\mathbb{I}_{(0,g(u+j\Delta))}(x)\mathbb{I}_{[0,\infty)}(u+j\Delta)\mathbb{I}_{(0,g(u+k\Delta))}(x)\mathbb{I}_{[0,\infty)}(u+k\Delta)dxdu\\\
&\propto\sum_{j=-\infty}^{\infty}\gamma(j\Delta)<\infty.\end{split}$ (22)

This condition is satisfied both for an exponential trawl function and also
for a supGamma trawl function with short memory. In the latter case, we have
that $\gamma(x)\propto(1+|x|/\alpha)^{1-H}$ for $\alpha>0,H>2$. Then, the
finiteness of (LABEL:fin) follows using the $\zeta$-function representation.

####  A.4.2 Verifying Assumption (15) from Theorem 7

We need to verify

$\displaystyle\sum_{j=-\infty}^{\infty}\left(\int_{\mathbb{R}\times\mathbb{R}}|f(x,u)||f(x,u+j\Delta)|dxdu\right)^{2}<\infty.$

Using very similar computations as before, we find that the above condition is
equivalent to

$\displaystyle\sum_{j=-\infty}^{\infty}\left(\int_{\mathbb{R}\times\mathbb{R}}|f(x,u)||f(x,u+j\Delta)|dxdu\right)^{2}$
$\displaystyle=\sum_{j=-\infty}^{\infty}\left(\int_{\mathbb{R}\times\mathbb{R}}\mathbb{I}_{(0,g(u)}(x)\mathbb{I}_{[0,\infty)}(u)\mathbb{I}_{(0,g(u+j\Delta))}(x)\mathbb{I}_{[0,\infty)}(u+j\Delta)dxdu\right)^{2}\propto\sum_{j=-\infty}^{\infty}\gamma^{2}(j\Delta)<\infty,$

which is satisfied by the exponential trawl function and the supGamma trawl
functions with $H>3/2$, which includes some long-memory settings.

## Acknowledgement

I would like to thank Paul Doukhan for suggesting a study of periodic trawl
processes and for helpful discussions, as well as Michele Nguyen for
commenting on an earlier version of this article.

## References

  * (1)
  * Alomari et al. (2020) Alomari, H. M., Ayache, A., Fradon, M. & Olenko, A. (2020), ‘Estimation of cyclic long-memory parameters’, Scandinavian Journal of Statistics. Theory and Applications 47(1), 104–133.
https://doi.org/10.1111/sjos.12404

  * Andel (1986) Andel, J. (1986), ‘Long memory time series models’, Kybernetika 22, 105–123.
  * Arteche (2002) Arteche, J. (2002), ‘Semiparametric robust tests on seasonal or cyclical long memory time series’, Journal of Time Series Analysis 23(3), 251–285.
https://onlinelibrary.wiley.com/doi/abs/10.1111/1467-9892.00264

  * Arteche & Robinson (2000) Arteche, J. & Robinson, P. M. (2000), ‘Semiparametric inference in seasonal and cyclical long memory processes’, Journal of Time Series Analysis 21(1), 1–25.
https://onlinelibrary.wiley.com/doi/abs/10.1111/1467-9892.00170

  * Ayache et al. (2022) Ayache, A., Fradon, M., Nanayakkara, R. & Olenko, A. (2022), ‘Asymptotic normality of simultaneous estimators of cyclic long-memory processes’, Electronic Journal of Statistics 16(1), 84–115.
https://doi.org/10.1214/21-ejs1953

  * Bacro et al. (2020) Bacro, J., Gaetan, C., Opitz, T. & Toulemonde, G. (2020), ‘Hierarchical space-time modeling of asymptotically independent exceedances with an application to precipitation data’, Journal of the American Statistical Association 115(530), 555–569.
https://doi.org/10.1080/01621459.2019.1617152

  * Barndorff-Nielsen (2011) Barndorff-Nielsen, O. E. (2011), ‘Stationary infinitely divisible processes’, Brazilian Journal of Probability and Statistics 25(3), 294–322.
https://doi.org/10.1214/11-BJPS140

  * Barndorff-Nielsen et al. (2013) Barndorff-Nielsen, O. E., Benth, F. E. & Veraart, A. E. D. (2013), ‘Modelling energy spot prices by volatility modulated Lévy-driven Volterra processes’, Bernoulli 19(3), 803–845.
http://www.jstor.org/stable/23525714

  * Barndorff-Nielsen et al. (2018) Barndorff-Nielsen, O. E., Benth, F. E. & Veraart, A. E. D. (2018), Ambit stochastics, Vol. 88 of Probability Theory and Stochastic Modelling, Springer, Cham.
https://doi.org/10.1007/978-3-319-94129-5

  * Barndorff-Nielsen et al. (2014) Barndorff-Nielsen, O. E., Lunde, A., Shephard, N. & Veraart, A. E. D. (2014), ‘Integer-valued trawl processes: A class of stationary infinitely divisible processes’, Scandinavian Journal of Statistics 41(3), 693–724.
  * Bennedsen et al. (2023) Bennedsen, M., Lunde, A., Shephard, N. & Veraart, A. E. D. (2023), ‘Inference and forecasting for continuous-time integer-valued trawl processes’, Journal of Econometrics 236(2), 105476.
https://www.sciencedirect.com/science/article/pii/S0304407623001926

  * Bennett (1958) Bennett, W. R. (1958), ‘Statistics of regenerative digital transmission’, The Bell System Technical Journal 37, 1501–1542.
https://doi.org/10.1002/j.1538-7305.1958.tb01560.x

  * Brockwell & Davis (1987) Brockwell, P. J. & Davis, R. A. (1987), Time series: theory and methods, Springer Series in Statistics, Springer-Verlag, New York.
https://doi.org/10.1007/978-1-4899-0004-3

  * Cohen & Lindner (2013) Cohen, S. & Lindner, A. (2013), ‘A central limit theorem for the sample autocorrelations of a Lévy driven continuous time moving average process’, Journal of Statistical Planning and Inference 143(8), 1295–1306.
https://www.sciencedirect.com/science/article/pii/S0378375813000670

  * Courgeau & Veraart (2022) Courgeau, V. & Veraart, A. E. (2022), ‘Asymptotic theory for the inference of the latent trawl model for extreme values’, Scandinavian Journal of Statistics 49(4), 1448–1495.
https://onlinelibrary.wiley.com/doi/abs/10.1111/sjos.12563

  * Curato & Stelzer (2019) Curato, I. V. & Stelzer, R. (2019), ‘Weak dependence and GMM estimation of supOU and mixed moving average processes’, Electronic Journal of Statistics 13(1), 310 – 360.
https://doi.org/10.1214/18-EJS1523

  * Das & Genton (2021) Das, S. & Genton, M. G. (2021), ‘Cyclostationary processes with evolving periods and amplitudes’, IEEE Transactions on Signal Processing 69, 1579–1590.
  * Doukhan et al. (2019) Doukhan, P., Jakubowski, A., Lopes, S. & Surgailis, D. (2019), ‘Discrete-time trawl processes’, Stochastic Processes and their Applications 129(4), 1326–1348.
https://www.sciencedirect.com/science/article/pii/S0304414918301571

  * Doukhan et al. (2020) Doukhan, P., Roueff, F. & Rynkiewicz, J. (2020), ‘Spectral estimation for non-linear long range dependent discrete time trawl processes’, Electronic Journal of Statistics 14(2), 3157 – 3191.
https://doi.org/10.1214/20-EJS1742

  * Espejo et al. (2015) Espejo, R. M., Leonenko, N. N., Olenko, A. & Ruiz-Medina, M. D. (2015), ‘On a class of minimum contrast estimators for Gegenbauer random fields’, TEST 24, 657–680.
https://doi.org/10.1007/s11749-015-0428-4

  * Fasen (2005) Fasen, V. (2005), ‘Extremes of regularly varying Lévy-driven mixed moving average processes’, Advances in Applied Probability 37(4), 993–1014.
  * Ferrara & Guégan (2001) Ferrara, L. & Guégan, D. (2001), Comparison of parameter estimation methods in cyclical long memory time series, in C. Junis, J. Moody & A. Timmermann, eds, ‘Development in Forecast Combination and Portfolio Choice’, Wiley, New York, USA.
  * Fuchs & Stelzer (2013) Fuchs, F. & Stelzer, R. (2013), ‘Mixing conditions for multivariate infinitely divisible processes with an application to mixed moving averages and the supOU stochastic volatility model’, ESAIM: Probability and Statistics 17, 455–471.
  * Gardner et al. (2006) Gardner, W. A., Napolitano, A. & Paura, L. (2006), ‘Cyclostationarity: Half a century of research’, Signal Processing 86(4), 639–697.
https://www.sciencedirect.com/science/article/pii/S0165168405002409

  * Genton & Hall (2007) Genton, M. G. & Hall, P. (2007), ‘Statistical inference for evolving periodic functions’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69(4), 643–657.
https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9868.2007.00604.x

  * Geweke & Porter-Hudak (1983) Geweke, J. & Porter-Hudak, S. (1983), ‘The estimation and application of long memory time series models’, Journal of Time Series Analysis 4(4), 221–238.
https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9892.1983.tb00371.x

  * Giraitis et al. (2001) Giraitis, L., Hidalgo, J. & Robinson, P. M. (2001), ‘Gaussian estimation of parametric spectral density with unknown pole’, The Annals of Statistics 29(4), 987–1023.
http://www.jstor.org/stable/2674066

  * Gladyšev (1961) Gladyšev, E. G. (1961), ‘Periodically correlated random sequences’, Doklady Akademii Nauk SSSR 137, 1026–1029.
  * Grahovac et al. (2018) Grahovac, D., Leonenko, N. N. & Taqqu, M. S. (2018), ‘Intermittency of trawl processes’, Statistics and Probability Letters 137, 235 – 242.
https://www.sciencedirect.com/science/article/pii/S0167715218300415

  * Gray et al. (1989) Gray, H. L., Zhang, N.-F. & Woodward, W. A. (1989), ‘On generalized fractional processes’, Journal of Time Series Analysis 10(3), 233–257.
https://doi.org/10.1111/j.1467-9892.1989.tb00026.x

  * Hall et al. (2000) Hall, P., Reimann, J. & Rice, J. (2000), ‘Nonparametric estimation of a periodic function’, Biometrika 87(3), 545–557.
http://www.jstor.org/stable/2673629

  * Hidalgo (2005) Hidalgo, J. (2005), ‘Semiparametric estimation for stationary processes whose spectra have an unknown pole’, The Annals of Statistics 33(4), 1843 – 1889.
https://doi.org/10.1214/009053605000000318

  * Hidalgo & Soulier (2004) Hidalgo, J. & Soulier, P. (2004), ‘Estimation of the location and exponent of the spectral singularity of a long memory process’, Journal of Time Series Analysis 25(1), 55–81.
https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9892.2004.00337.x

  * Hosking (1981) Hosking, J. R. M. (1981), ‘Fractional differencing’, Biometrika 68(1), 165–176.
http://www.jstor.org/stable/2335817

  * Hurd (1969) Hurd, H. L. (1969), An investigation of periodically correlated stochastic processes, PhD thesis, Duke University, Department of Electrical Engineering.
  * Hurd & Miamee (2007) Hurd, H. L. & Miamee, A. (2007), Periodically correlated random sequences, Wiley Series in Probability and Statistics, Wiley-Interscience [John Wiley & Sons], Hoboken, NJ.  Spectral theory and practice.
https://doi.org/10.1002/9780470182833

  * Leonte & Veraart (2022) Leonte, D. & Veraart, A. E. D. (2022), ‘Simulation methods and error analysis for trawl processes and ambit fields’.
https://arxiv.org/abs/2208.08784

  * Maddanu (2022) Maddanu, F. (2022), ‘A harmonically weighted filter for cyclical long memory processes’, AStA Advances in Statistical Analysis 106, 49–78.
https://doi.org/10.1007/s10182-021-00394-9

  * Mátyás (1999) Mátyás, L., ed. (1999), Generalized method of moments estimation, Cambridge University Press, Cambridge.