# Bootstrapping Whittle Estimators

Jens-Peter Kreiss  Technische Universität Braunschweig, Institut für
Mathematische Stochastik, Universitätsplatz 2, D–38106 Braunschweig, Germany.
and  Efstathios Paparoditis  University of Cyprus, Department of Mathematics
and Statistics, 1678 Nicosia, Cyprus.

###### Abstract.

Fitting parametric models by optimizing frequency domain objective functions
is an attractive approach of parameter estimation in time series analysis.
Whittle estimators are a prominent example in this context. Under weak
conditions and the (realistic) assumption that the true spectral density of
the underlying process does not necessarily belong to the parametric class of
spectral densities fitted, the distribution of Whittle estimators typically
depends on difficult to estimate characteristics of the underlying process.
This makes the implementation of asymptotic results for the construction of
confidence intervals or for assessing the variability of estimators, difficult
in practice. This paper proposes a frequency domain bootstrap method to
estimate the distribution of Whittle estimators which is asymptotically valid
under assumptions that not only allow for (possible) model misspecification
but also for weak dependence conditions which are satisfied by a wide range of
stationary stochastic processes. Adaptions of the bootstrap procedure
developed to incorporate different modifications of Whittle estimators
proposed in the literature, like for instance, tapered, de-biased or boundary
extended Whittle estimators, are also considered. Simulations demonstrate the
capabilities of the bootstrap method proposed and its good finite sample
performance. A real-life data analysis also is presented.

###### Key words and phrases:

Periodogram, Bootstrap, Linear process, Nonparametric kernel estimation,
Spectral means.

###### 2000 Mathematics Subject Classification:

Primary 62M10, 62M15; secondary 62G09

##  1\. Introduction

Optimizing frequency domain objective functions is an attractive approach of
fitting parametric models to observed time series. Let
$X_{1},X_{2},\ldots,X_{n}$ be a time series stemming from a stationary process
$\\{X_{t},t\in{\mathbb{Z}}\\}$ and assume that this process possesses a
spectral density $f$. Let ${\mathcal{F}}_{\theta}$ be a family of parametric
spectral densities, where $f_{\theta}\in{\mathcal{F}}_{\theta}$ is determined
by a $m$-dimensional parameter vector $\theta$ and suppose that we are
interested in fitting to $X_{1},X_{2},\ldots,X_{n}$ a model from the class
${\mathcal{F}}_{\theta}$. Notice that we do not assume
$f\in{\mathcal{F}}_{\theta}$, that is, we allow for the practical important
case of model misspecification where the parametric model class considered
does not necessarily contain the true spectral density $f$.

Several approaches for selecting a frequency domain objective function and
consequently for developing a frequency domain procedure to fit parametric
models exist; we refer here to Taniguchi (1987) and to Dahlhaus and Wefelmeyer
(1996) for examples. In this context, Whittle estimators play an important
role. This is due to the fact that Whittle estimators are computationally fast
and they are obtained via minimizing a frequency domain approximation of the
(Gaussian) log likelihood function. Moreover, for a variety of models and
under different assumptions, Whittle estimators are asymptotically normal,
asymptotically equivalent to the exact maximum likelihood estimators and
asymptotically efficient (in Fisher-sense) if $f\in{\mathcal{F}}_{\theta}$;
see Sykulski et al. (2019) and Subba Rao and Yang (2020) for a recent
discussion of the related literature. Furthermore and even if
$f\notin{\mathcal{F}}_{\theta}$, Whittle estimators retain certain efficiency
properties; Dahlhaus and Wefelmeyer (1996). Despite these nice properties,
however, the limiting distribution of Whittle’s estimators is affected by
characteristics of the process which make the implementation of asymptotic
results for assessing their variability or for constructing confidence
intervals, difficult in practice. To elaborate, allowing for possible model
misspecification and avoiding restrictive structural assumptions for the
underlying process class, like for instance linearity assumptions, the
limiting distribution of Whittle estimators typically depends on the
parametric spectral density from the class ${\mathcal{F}}_{\theta}$ which
“best fits” the data, say $f_{\theta_{0}}$, the true spectral density $f$ as
well as the entire fourth order cumulant structure of the underlying process
$\\{X_{t};t\in{\mathbb{Z}}\\}$. Estimation of the last mentioned quantity is a
rather difficult problem.

In situations like the above, bootstrapping may offer an alternative to
classical large sample approximations. Bootstrapping Whittle estimators has
been discussed in the literature under certain structural assumptions on the
underlying process $\\{X_{t},t\in{\mathbb{Z}}\\}$. It is typically assumed
that $\\{X_{t},t\in{\mathbb{Z}}\\}$ is a linear process driven by i.i.d.
innovations and that Kolmogorov’s formula holds true; see Dahlhaus and Janas
(1996) and Kim and Nordman (2013). Such structural assumptions lead to a
simplification of the limiting distribution of Whittle estimators and
consequently of the features of the underlying linear process that the
bootstrap procedure has to appropriately mimic in order to be consistent. To
elaborate, notice first that, under linearity assumptions, the problem of
estimating the variance $\sigma^{2}$ of the i.i.d. innovations driving the
linear process can be separated from the problem of estimating the remaining
coefficients of the parameter vector, denoted by $\tau\in{\mathbb{R}}^{m-1}$,
i.e., $\theta=(\sigma^{2},\tau)$. Second and more importantly, under the
linearity assumption, the limiting distribution of the Whittle estimator of
the parameter part $\tau$, does not depend on the fourth order cumulants of
the process; see Section 2 for details. Therefore, if one is solely interested
in estimating the distribution of the Whittle estimator of the innovation free
part $\tau$, then standard frequency domain bootstrap procedures which
generate pseudo periodogram ordinates that are independent across frequencies,
can successfully be applied. In this context, the multiplicative bootstrap,
see Hurvich an Zeger (1987), Franke and Härdle (1992) and Dahlhaus and Janas
(1996), or the local periodogram bootstrap, Paparoditis and Politis (1999),
are consistent. However, the pseudo periodogram ordinates generated by the
aforementioned bootstrap procedures are independent across frequencies.
Therefore, these bootstrap procedures are not able to imitate the covariance
structure of the periodogram ordinates which is responsible for the fact that
the fourth structure of the process shows up in the limiting distribution of
Whittle estimators. As a consequence, these procedures fail in all cases where
the fourth order cumulants of the process affect the distribution of interest.

In this paper we present a frequency domain, hybrid bootstrap procedure for
Whittle estimators which is valid under weak assumptions on the underlying
process $\\{X_{t},t\in{\mathbb{Z}}\\}$ and which also covers the practical
important case where the true spectral density $f$ does not necessarily belong
to the parametric class ${\mathcal{F}}_{\theta}$. The procedure consists of
two main parts: A multiplicative frequency domain bootstrap part and a part
based on the convolution of resampled periodograms of subsamples. The latter
is an adaption to the frequency domain of the convolved subsampling idea
proposed in Tewes et al. (2019). The two parts contribute differently and
complementary in estimating the distribution of interest. The multiplicative
part of the bootstrap procedure is used to estimate all features of the
distribution of Whittle estimators including the parts of the limiting
covariance matrix that depend on the second order structure of the underlying
process and of the parametric model fitted to the time series at hand.
However, the components of the covariance matrix of these estimators that
depend on the fourth order structure of the underlying process
$\\{X_{t},t\in{\mathbb{Z}}\\}$, are estimated using the part of the bootstrap
procedure which is based on the convolution of resampled periodograms of
subsamples. Putting the two parts together in an appropriate way, leads to a
bootstrap procedure which is asymptotically valid under conditions on the
dependence structure of the process $\\{X_{t},t\in{\mathbb{Z}}\\}$ which go
far beyond linearity and at the same time appropriately captures the effects
of (possible) model misspecification on the distribution of interest.

The frequency domain, hybrid approach proposed in this paper and which uses
convolution of resampled periodograms of subsamples together with the
multiplicative periodogram bootstrap, is related to the proposal of Meyer et
al. (2020). However and additional to differences in the technical tools used
to establish bootstrap consistency, the merging of the multiplicatve and of
the convolved part of the bootstrap procedure presented here, is different,
more involved and tailormade for Whittle estimators. Furthermore, we show, how
the described bootstrap procedure can appropriately be modified, respectively,
extended to incorporate several modifications of Whittle estimators which have
been proposed in order to improve the finite (small) sample bias of these
estimators. This concerns tapered Whittle estimators, Dahlhaus (1988), so-
called de-biased Whittle estimators, Sykulski et al. (2019), and improvements
of Whittle’s quasi Gaussian likelihood approximation based on boundary
corrected periodograms; see Subba Rao and Yang (2020). The corresponding
extensions of the frequency domain bootstrap proposed in this context are
novel and of interest on their own.

The paper is organized as follows. In Section 2 we review some basic results
on Whittle estimators which are important for our subsequent development of
the bootstrap. The basic frequency domain bootstrap procedure introduced is
presented in Section 3. Section 4 is devoted to the derivation of theoretical
results and establishes consistency of the bootstrap. Section 5 deals with the
extensions of the basic bootstrap procedure in order to incorporate the
aforementioned modifications of standard Whittle estimators. Section 6
discusses some issues related to the practical implementation of the bootstrap
algorithm, presents simulations which investigate the finite sample
performance of the new bootstrap method and make comparisons with the
asymptotic Gaussian approximation, respectively, the multiplicative
periodogram bootstrap. A real-life data application also is discussed.
Auxiliary lemmas as well as proofs of the main results are deferred to Section
7.

##  2\. Whittle Estimators

In his PhD thesis Whittle introduced a frequency domain approximation of the
log-likelihood function of a stationary Gaussian time series (cf. Whittle
(1951) and Whittle (1953)). This approximation can be written as

$\displaystyle l_{n}(\theta)$
$\displaystyle=2n\log(2\pi)+\sum_{j\in{\mathcal{F}}_{n}}\Big{\\{}\log
f(\lambda_{j,n})+\frac{\displaystyle I_{n}(\lambda_{j,n})}{\displaystyle
f(\lambda_{j,n})}\Big{\\}},$

where $I_{n}(\lambda_{j,n})$, denotes the periodogram of the time series
$X_{1},X_{2},\ldots,X_{n}$ evaluated at the Fourier frequency
$\lambda_{j,n}=2\pi j/n\in{\mathcal{F}}_{n}$ and
${\mathcal{F}}_{n}=\\{-[(n-1)/2],\ldots,[n/2]\\}$ is the set of Fourier
frequencies. Ignoring the first additive term and approximating the integral
over the set ${\mathcal{G}}(n)=\\{-N,-N+1,\ldots,-1,1,\ldots,N\\}$, where
$N=[n/2]$, Whittle’s approximation to the log-likelihood function used in this
paper, is given by

(2.1) $D_{n}(\theta,I_{n})=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{\\{}\log
f_{\theta}(\lambda_{j,n})+\frac{\displaystyle
I_{n}(\lambda_{j,n})}{\displaystyle f_{\theta}(\lambda_{j,n})}\Big{\\}}.$

A minimizer $\widehat{\theta}_{n}$ of

(2.2) $\theta\mapsto D_{n}(\theta,I_{n}),$

is called a Whittle estimator of $\theta$. Note that (2.1) can be considered
as a Riemann sum approximation of

(2.3) $D(\theta,I_{n})=\frac{1}{2\pi}\int_{-\pi}^{\pi}\Big{\\{}\log
f_{\theta}(\lambda)+\frac{\displaystyle I_{n}(\lambda)}{\displaystyle
f_{\theta}(\lambda)}\Big{\\}}d\lambda.$

Let

(2.4) $D(\theta,f)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\Big{\\{}\log
f_{\theta}(\lambda)+\frac{\displaystyle f(\lambda)}{\displaystyle
f_{\theta}(\lambda)}\Big{\\}}d\lambda$

and assume that

(2.5) $\theta_{0}=\arg\min_{\theta\in\Theta}D(\theta,f),$

exists and is unique. $\widehat{\theta}_{n}$ is then an estimator of
$\theta_{0}$ and $f_{\theta_{0}}$ denotes the spectral density from the
parametric family ${\mathcal{F}}_{\theta}$ which best fits the spectral
density $f$ of the underlying process $\\{X_{t},t\in{\mathbb{Z}}\\}$ in the
sense of minimizing the divergence measure (2.4). Note that if
$f\notin{\mathcal{F}}_{\theta}$, then the best approximating parametric
spectral density $f_{\theta_{0}}$ from the class ${\mathcal{F}}_{\theta}$,
clearly depends on the particular divergence measure $D(\theta,f)$ associated
with Whittle’s approximation of the log-likelihood function. Observe that
$\theta_{0}$ which minimizes $D(\theta,f)$ is the same as the one which
minimizes the so-called Kullback-Leibler information divergence. The later is
given for Gaussian processes by

$\frac{1}{2\pi}\int_{-\pi}^{\pi}\Big{\\{}\log\frac{f_{\theta}(\lambda)}{f(\lambda)}+\frac{\displaystyle
f(\lambda)}{\displaystyle f_{\theta}(\lambda)}-1\Big{\\}}d\lambda;$

see Dahlhaus and Wefelmeyer (1996). Different modifications of Whittle’s
likelihood approximation (2.3) have been proposed in the literature in order
to improve the finite sample behavior and more specifically the bias of the
estimator $\widehat{\theta}_{n}$. We mention here the tapered Whittle
likelihood proposed by Dahlhaus (1988), the debiased Whittle likelihood
proposed by Sykulski et al. (2019), the boundary corrected and the hybrid
Whittle likelihood approximation proposed by Subba Rao and Yang (2020). For
the sake of a better presentation, however, we will first focus on the
standard Whittle likelihood approximation (2.1), respectively, (2.3). Later
on, we will elaborate on how to appropriately modify the basic bootstrap
procedure proposed in order to take into account the aforementioned
modifications/extensions of Whittle estimators.

Let us briefly review the main ideas involved in deriving the limiting
distribution of the estimator $\widehat{\theta}_{n}$ and which are important
for our subsequent discussion of the bootstrap. Consider (2.3) and (2.4 ),
assume that $f_{\theta}$ is sufficiently smooth with respect to $\theta$ and
recall that $\widehat{\theta}_{n}$ and $\theta_{0}$ satisfy the score
equations

(2.6)
$\frac{\displaystyle\partial}{\displaystyle\partial\theta}D_{n}(\theta,I_{n})\Big{|}_{\theta=\widehat{\theta}_{n}}=0\
\ \ \mbox{and}\ \
\frac{\displaystyle\partial}{\displaystyle\partial\theta}D(\theta,f)\Big{|}_{\theta=\theta_{0}}=0,$

respectively. Using a linear approximation of
$(\partial/\partial\theta)D_{n}(\theta,I_{n})\big{|}_{\theta=\widehat{\theta}_{n}}$
around $\theta_{0}$ we get, taking into account (2.6), that

$0=\frac{\displaystyle\partial}{\displaystyle\partial\theta}D_{n}(\theta_{0},I_{n})+\frac{\displaystyle\partial^{2}}{\displaystyle\partial\theta\partial\theta^{\top}}D_{n}(\theta_{0},I_{n})(\widehat{\theta}_{n}-\theta_{0})+R_{n}.$

Notice that for simplicity, the notation
$(\partial/\partial\theta)D_{n}(\theta_{0},I_{n})$ for
$(\partial/\partial\theta)D_{n}(\theta,I_{n})\big{|}_{\theta=\theta_{0}}$ has
been used with an analogue notation for the matrix of second order partial
derivatives
$\partial^{2}\big{/}\partial\theta\partial\theta^{\top}D_{n}(\theta_{0},I_{n})$.
Provided that the remainder $R_{n}$ is $o_{P}(n^{-1/2})$ and that
$(\partial^{2}/\partial\theta\partial\theta^{\top})D_{n}(\theta_{0},I_{n})$ is
invertible, the following basic expression is then obtained,

(2.7)
$\sqrt{n}\big{(}\widehat{\theta}_{n}-\theta_{0}\big{)}=\Big{(}\frac{\displaystyle\partial^{2}}{\displaystyle\partial\theta\partial\theta^{\top}}D_{n}(\theta_{0},I_{n})\Big{)}^{-1}\sqrt{n}\int_{-\pi}^{\pi}g_{\theta_{0}}(\lambda)\big{(}I_{n}(\lambda)-f(\lambda)\big{)}d\lambda+o_{P}(1).$

Here $g_{\theta_{0}}(\lambda)$ is a $m$-dimensional vector of scores, the
$j$th element of which is given by

$g_{j,\theta_{0}}(\lambda)=\frac{1}{2\pi
f_{\theta_{0}}(\lambda)}\frac{\partial}{\partial\theta_{j}}\log
f_{\theta_{0}}(\lambda)=-\frac{1}{2\pi}\frac{\partial}{\partial\theta_{j}}f^{-1}_{\theta_{0}}(\lambda),$

$j=1,2,\ldots,m$ and $f^{-1}_{\theta}=1/f_{\theta}$. Equation (2.7) suggests
that the distribution of Whittle’s estimator $\widehat{\theta}_{n}$, can be
well approximated by the product of the inverse of the $m\times m$ random
matrix
$W_{n}=\big{(}\partial^{2}/\partial\theta_{j}\partial\theta_{k}D_{n}(\theta_{0},I_{n})\big{)}_{j,k=1,2,\ldots,m}$
with the $m$-dimensional vector of integrated periodograms
$\big{(}\sqrt{n}\int_{-\pi}^{\pi}g_{j,\theta_{0}}(\lambda)\big{(}I_{n}(\lambda)-f(\lambda)\big{)}d\lambda,j=1,2,\ldots,m\big{)}^{\top}$.

Based on expression (2.7), asymptotic theory for Whittle estimators has been
developed in the literature under a variety of assumptions on the dependence
structure of the underlying process. In early papers Walker (1964) and Hannan
(1973) derived asymptotic normality under linearity assumptions while Hosoya
(1979) also allowed for long range dependence. For a more recent review of the
related literature as well as derivations of asymptotic normality based on a
physical dependence measure, we refer to Shao (2010). It is typically shown,
that,

(2.8)
$\sqrt{n}\big{(}\widehat{\theta}_{n}-\theta_{0}\big{)}\stackrel{{\scriptstyle
D}}{{\rightarrow}}{\mathcal{N}}\big{(}0,W^{-1}(V_{1}+V_{2})W^{-1}\big{)},$

as $n\rightarrow\infty$, where the $m\times m$ matrices $W$, $V_{1}$ and
$V_{2}$ are given by

$W=\Big{(}\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}D(\theta_{0},f)\Big{)}_{j,k=1,2,\ldots,m},$
$V_{1}=\Big{(}4\pi\int_{-\pi}^{\pi}g_{j,\theta_{0}}(\lambda)g_{k,\theta_{0}}(\lambda)f^{2}(\lambda)d\lambda\Big{)}_{j,k=1,2,\ldots,m}$

and

$V_{2}=\Big{(}2\pi\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}g_{j,\theta_{0}}(\lambda_{1})g_{k,\theta_{0}}(\lambda_{2})f_{4}(\lambda_{1},\lambda_{2},-\lambda_{2})d\lambda_{1}d\lambda_{2}\Big{)}_{j,k=1,2,\ldots,m}.$

Observe that $W$ and $V_{1}$ only depend on the parametric and the true
spectral densities, that is on $f_{\theta_{0}}$ and $f$, while the matrix
$V_{2}$ depends on $f_{\theta_{0}}$ and on the fourth order cumulant spectral
density $f_{4}$ of the underlying process $\\{X_{t},t\in{\mathbb{Z}}\\}$. The
latter is defined as

$f_{4}(\lambda_{1},\lambda_{2},\lambda_{3})=\frac{1}{(2\pi)^{3}}\sum_{h_{1},h_{2},h_{3}\in{\mathbb{Z}}}cum(X_{0},X_{h_{1}},X_{h_{2}},X_{h_{3}})\exp^{-i(h_{1}\lambda_{1}+h_{2}\lambda_{2}+h_{3}\lambda_{3})},$

where $cum(X_{0},X_{h_{1}},X_{h_{2}},X_{h_{3}})$ denotes the fourth order
cumulant of the process $\\{X_{t};t\in{\mathbb{Z}}\\}$, cf. Rosenblatt (1985).
Notice that if $f=f_{\theta_{0}}$, that is, if the spectral density of the
process belongs to the parametric family ${\mathcal{F}}_{\theta}$ and the
model is correctly specified, then $V_{1}=2W$ and the covariance matrix of the
limiting Gaussian distribution (2.8) is given by $W^{-1}(2I_{m}+V_{2}W^{-1})$,
where $I_{m}$ denotes the $m\times m$ unit matrix.

The matrix $V_{2}$ in general does not entirely disappear even if the
underlying process is linear, that is, if $X_{t}$ is generated as
$X_{t}=\sum_{j=-\infty}^{\infty}\psi_{j}\varepsilon_{t-j}$, where
$\sum_{j=-\infty}^{\infty}|\psi_{j}|<\infty$ and the $\varepsilon_{t}$’s are
zero mean, i.i.d. innovations with variance $\sigma^{2}>0$. However, in this
case, the matrix $V_{2}$ simplifies considerably. To elaborate, recall that
under the assumption of linearity and of validity of Kolmogorov’s formulae,
(see Blockwell and Davis (1991), Ch. 5.8), $\theta=(\sigma^{2},\tau)$ where
$\tau\in{\mathbb{R}}^{m-1}$ is free of the innovation variance $\sigma^{2}$.
The dependence on the fourth order moment structure of Whittle estimators
disappears then if one is solely interested in the distribution of the
estimators of the part $\tau$ of the parameter vector. This is due to the fact
that for linear processes, the spectral density $f_{\theta}(\cdot)$ factorizes
as $f_{\theta}(\cdot)=h_{\tau}(\cdot)\sigma^{2}/(2\pi)$, where the function
$h_{\tau}(\cdot)$ depends on $\tau$ only. Then, and since for the same class
of processes,
$f_{4}(\lambda_{1},\lambda_{2},-\lambda_{2})=(2\pi)^{-1}\eta_{4}f(\lambda_{1})f(\lambda_{2})$,
with $\eta_{4}=E(\varepsilon^{4}_{1}/\sigma^{4}-3)$, the rescaled fourth order
cumulants (kurtosis) of the i.i.d. innovations, we get that,

$\displaystyle 2\pi\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}$ $\displaystyle
g_{j,\tau_{0}}(\lambda_{1})g_{k,\tau_{0}}(\lambda_{2})f_{4}(\lambda_{1},\lambda_{2},-\lambda_{2})d\lambda_{1}d\lambda_{2}$
$\displaystyle=\eta_{4}\int_{-\pi}^{\pi}g_{j,\tau_{0}}(\lambda_{1})f(\lambda_{1})d\lambda_{1}\int_{-\pi}^{\pi}g_{k,\tau_{0}}(\lambda_{2})f(\lambda_{2})\lambda_{2}=0.$

The last equality follows since the score function implies
$\int_{-\pi}^{\pi}g_{s,\tau_{0}}(\lambda)f(\lambda)d\lambda=0$ for every
$s=1,2,\ldots,m-1$, where $g_{s,\tau_{0}}(\cdot)$ denotes the partial
derivative of $-h^{-1}_{\tau}(\cdot)/(2\pi)$ with respect to the $s$-th
variable of the $m-1$ dimensional vector $\tau$, evaluated at $\tau=\tau_{0}$.
Thus the $(m-1)\times(m-1)$ submatrix of $V_{2}$ which corresponds to the
elements of the vector $\tau$ only, consists of zeros. Notice that this
simplification does not hold true for the components of the matrix $V_{2}$
which are affected by the estimator of $\sigma^{2}$; see Dahlhaus and Janas
(1996) for more details.

Now, a close look at the derivations leading to the covariance formulae
$W^{-1}(V_{1}+V_{2})W^{-1}$ of the limiting Gaussian distribution (2.8),
reveals that the term $V_{2}$ is solely due to the weak and asymptotically
vanishing covariance of the periodogram ordinates across frequencies. Recall
the basic expression of the covariance of the periodogram for nonzero Fourier
frequencies $|\lambda_{j,n}|\neq|\lambda_{k,n}|$,

(2.9) ${\rm
cov}(I_{n}(\lambda_{j,n}),I_{n}(\lambda_{k,n}))=\frac{1}{n}f_{4}(\lambda_{j,n},\lambda_{k,n},-\lambda_{k,n})(1+o(1))+O(n^{-2}).$

Summing up these covariances over all frequencies
$|\lambda_{j,n}|\neq|\lambda_{k,n}|$ in the set ${\mathcal{G}}(n)$, leads to a
non-vanishing contribution to the limiting distribution of Whittle estimators
as expressed by the matrix $V_{2}$. As already mentioned, frequency domain
bootstrap procedures which generate independent periodogram ordinates, like
for instance the multiplicative perodogram bootstrap, can not imitate the weak
dependence structure (2.9). Therefore, and besides the special case $f_{4}=0$,
e.g. Gaussian time series, such procedures do not appropriately capture the
term $V_{2}$ and they fail in consistently estimating the distribution of
$\sqrt{n}(\widehat{\theta}_{n}-\theta_{0})$.

##  3\. The basic bootstrap procedure

Our goal is to develop a consistent bootstrap estimator of the distribution of
$L_{n}=\sqrt{n}(\widehat{\theta}_{n}-\theta_{0})$ without imposing restrictive
structural assumptions on the process $\\{X_{t},t\in{\mathbb{Z}}\\}$ and
allowing at the same time for the case of model misspecification. Toward this
goal, we propose a hybrid, frequency domain bootstrap procedure which builds
upon the multiplicative periodogram bootstrap and appropriately extends it, in
order to overcome its limitations for the class of Whittle estimators. The
procedure consists of two main parts. The first part is based on the
multiplicative bootstrap approach as proposed by Franke and Härdle (1992) and
Dahlhaus and Janas (1996). It is used to estimate all features of the
distribution of $L_{n}$ including the parts of the covariance matrix which do
not depend on the fourth order characteristics of the process
$\\{X_{t};t\in{\mathbb{Z}}\\}$ and in particular the matrices $W$ and $V_{1}$.
This is done in Step 2 and Step 3 of the following algorithm. However, and as
already mentioned, since the multiplicative periodogram bootstrap generates
independent pseudo periodogram ordinates, it can not be used to imitate the
fourth order characteristics of the process which affect the distribution of
$L_{n}$ and more specifically the matrix $V_{2}$. The second part of our
bootstrap procedure corrects for this shortcoming. This is achieved by
generating pseudo periodograms of subsamples of length $b$, $b<n$, using
randomly selected sets of appropriately defined frequency domain residuals.
The advantage of these pseudo periodograms of subsamples is that they retain
the weak dependence structure of the periodogram across the Fourier
frequencies corresponding to the subsamples. They can, therefore, be used to
consistently estimate the missing part $V_{2}$. This is done in Step 4 and
Step 5 of the bootstrap algorithm presented bellow. Putting these two parts
together in an appropriate way, finally, leads in Step 6, to a consistent
estimator of the distribution of the random sequence
$\sqrt{n}(\widehat{\theta}_{n}-\theta_{0})$ of interest.

The following algorithm implements the above ideas and is the basic hybrid
frequency domain bootstrap procedure proposed in this paper.

  1. Step 1:

Calculate Whittle’s estimator $\widehat{\theta}_{n}$.

  2. Step 2:

Let $\widehat{f}$ be a nonparametric estimator of $f$. For $j=1,2,\ldots,N$,
generate

$I^{\ast}_{n}(\lambda_{j,n})=\widehat{f}(\lambda_{j,n})\cdot U^{\ast}_{j},$

where $U_{j}^{\ast}$ are i.i.d. standard exponential distributed random
variables. Set $I^{\ast}_{n}(\lambda_{j,n})=I^{\ast}_{n}(-\lambda_{j,n})$ for
$j=-1,-2,\ldots,-N$.

  3. Step 3:

Let

$\widehat{\theta}^{\ast}_{n}=\arg\min_{\theta\in\Theta}\,D_{n}(\theta,I^{\ast}_{n})\
\mbox{and}\
\widehat{\theta}_{0}=\arg\min_{\theta\in\Theta}\,D_{n}(\theta,\widehat{f}).$

Define

$V_{1,n}^{\ast}=Var^{\ast}(M_{n}^{\ast})\ \ \mbox{and}\ \
W_{n}^{\ast}=\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}D_{n}(\theta,I^{*}_{n})\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)},$

where

$M^{\ast}_{n}=\frac{2\pi}{\sqrt{n}}\sum_{j\in{\mathcal{G}}(n)}g_{\widehat{\theta}_{0}}(\lambda_{j,n})\big{(}I_{n}^{\ast}(\lambda_{j,n})-\widehat{f}(\lambda_{j,n})\big{)},$

and $g_{\widehat{\theta}_{0}}(\lambda)$ the $m$-dimensional vector

$g_{\widehat{\theta}_{0}}(\lambda)=\Big{(}g_{j,\widehat{\theta}_{0}}(\lambda)=-\frac{1}{2\pi}\frac{\partial}{\partial\theta_{j}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\widehat{\theta}_{0}},\
j=1,2,\ldots,m\Big{)}^{\top}.$

Calculate the pseudo random vector $Z_{n}^{\ast}$ defined by

$Z^{\ast}_{n}=\big{(}V^{\ast}_{1,n}\big{)}^{-1/2}W^{\ast}_{n}\sqrt{n}(\widehat{\theta}^{*}_{n}-\widehat{\theta}_{0}).$

  4. Step 4:

Select an integer $b<n$ and generate $k=[n/b]$ pseudo periodograms pseudo
periodograms $I_{b}^{(\ell)}(\lambda_{j,b})$, $\ell=1,2,\ldots,k$, where

$I_{b}^{(\ell)}(\lambda_{j,b})=\widehat{f}(\lambda_{j,b})\cdot
U^{(\ell)}_{b}(\lambda_{j,b}).$

Here,
$U^{(\ell)}_{b}(\lambda_{j,b})=I_{b}^{(\ell)}(\lambda_{j,b})/\widetilde{f}(\lambda_{j,b})$,
where

$\hskip
42.67912ptI_{b}^{(\ell)}(\lambda_{j,b})=\frac{1}{2\pi}\big{|}\sum_{s=1}^{b}X_{i_{\ell}+s-1}e^{-is\lambda_{j,b}}\big{|}^{2}$

and $i_{\ell}$, $\ell=1,2,\ldots,k$, are i.i.d. random variables uniformly
distributed on the set $\\{1,2,\ldots,n-b+1\\}$. Furthermore,

$\widetilde{f}(\lambda_{j,b})=\frac{1}{n-b+1}\sum_{t=1}^{n-b+1}I_{b}^{(t)}(\lambda_{j,b}),$

with
$I_{b}^{(t)}(\lambda)=(2\pi)^{-1}\big{|}\sum_{s=1}^{b}X_{t+s-1}e^{-it\lambda}\big{|}^{2}$.
Calculate the pseudo random variables $M_{n}^{+}$ as

$M^{+}_{n}=\sqrt{kb}\frac{1}{k}\sum_{l=1}^{k}\frac{2\pi}{b}\sum_{j\in{\mathcal{G}}(b)}g_{\widehat{\theta}_{0}}(\lambda_{j,b})\big{(}I_{b}^{(\ell)}(\lambda_{j,b})-\widehat{f}(\lambda_{j,b})\big{)}.$

  5. Step 5:

Calculate the $m\times m$ matrix $V_{2,n}^{+}$ as

$V_{2,n}^{+}=\Sigma^{+}_{n}-C_{n}^{+},$

where $\Sigma^{+}_{n}={\rm Var}^{\ast}(M_{n}^{+})$,
$C_{n}^{+}=\big{(}c_{n}^{(r,s)}\big{)}_{r,s=1,2,\ldots,m}$ and the elements
$c_{n}^{(r,s)}$ are given by

$\displaystyle c_{n}^{(r,s)}$
$\displaystyle=\frac{8\pi^{2}}{b}\sum_{j\in{\mathcal{G}}(b)}g_{r,\widehat{\theta}_{0}}(\lambda_{j,b})g_{s,\widehat{\theta}_{0}}(\lambda_{j,b})\widehat{f}(\lambda_{j,b})^{2}$
$\displaystyle\ \ \ \ \
\times\Big{(}\frac{1}{n-b+1}\sum_{t=1}^{n-b+1}\frac{I_{b}^{(t)}(\lambda_{j,b})^{2}}{\widetilde{f}_{b}(\lambda_{j,b})^{2}}-1\Big{)}.$

  6. Step 6:

Appoximate the distribution of
$L_{n}=\sqrt{n}(\widehat{\theta}_{n}-\theta_{0})$ by the distribution of

$L_{n}^{\ast}=\big{(}W_{n}^{\ast}\big{)}^{-1}\big{(}V_{1,n}^{\ast}+V_{2,n}^{+}\big{)}^{1/2}\cdot
Z^{\ast}_{n}.$

Several aspects of the above bootstrap algorithm are clarified by the
following series of remarks and comments.

######  Remark 3.1.

  1. (i)

Observe that in order to appropriately capture the effect of model
misspecification, i.e., the case $f\notin{\mathcal{F}}_{\theta}$, a
nonparametric estimator $\widehat{f}$ is used in Step 2 to generate the pseudo
periodogram ordinates $I^{\ast}_{n}$ and not the estimated parametric spectral
density $f_{\widehat{\theta}_{n}}$.

  2. (ii)

The estimator $\widehat{\theta}_{0}$ in Step 3 is defined in a way which
imitates the properties of $\theta_{0}$. This estimator also delivers the
appropriate centering of the bootstrap estimator
$\widehat{\theta}_{n}^{\ast}$, that is
$\sqrt{n}(\widehat{\theta}_{n}^{\ast}-\widehat{\theta}_{0})$ is used as a
bootstrap analogue of $\sqrt{n}(\widehat{\theta}_{n}-\theta_{0})$.

  3. (iii)

As the proof of Theorem 4.1 shows,

$\sqrt{n}(\widehat{\theta}^{\ast}_{n}-\widehat{\theta}_{0})\stackrel{{\scriptstyle
D}}{{\rightarrow}}{\mathcal{N}}(0,W^{-1}V_{1}W^{-1}),$

in probability. Furthermore, Lemma 7.2 of Section 7, shows that
$W_{n}^{\ast}\stackrel{{\scriptstyle P}}{{\rightarrow}}W$ and
$V_{1,n}^{\ast}\stackrel{{\scriptstyle P}}{{\rightarrow}}V_{1}$. These facts
imply that the limiting distribution of the standardized pseudo random
variable
$Z^{\ast}_{n}=\big{(}V^{\ast}_{1,n}\big{)}^{-1/2}W^{\ast}_{n}\sqrt{n}(\widehat{\theta}^{*}_{n}-\widehat{\theta}_{0})$
appearing at the end of Step 3, has covariance matrix the $m\times m$ identity
matrix.

######  Remark 3.2.

  1. (i)

In Step 4 the resampled periodograms of subsamples of length $b$, i.e.,
$I_{b}^{(\ell)}(\lambda_{j,b})$, are obtained by using the same spectral
density estimator $\widehat{f}$ as in Step 2 but evaluated at the Fourier
frequencies $\lambda_{j,b}$ corresponding to the length $b$ of the subsamples.
Furthermore, in order to generate $I_{b}^{(\ell)}(\lambda_{j,b})$,
$j=1,2,\ldots,B$, the estimated spectral density $\widehat{f}(\lambda_{j,b})$
is multiplied with the entire set of frequency domain residuals, denoted by
$U_{b}^{(\ell)}(\lambda_{j,b})$, $j=1,2,\ldots,B$. Because these residuals are
obtained by rescaling the entire set of periodogram ordinates of a subsample,
they retain for the Fourier frequencies $\lambda_{j,b}$ of the subsample, the
weak dependence structure of the periodogram. Notice that the particular
rescaling of these residuals applied, ensures that
$E^{\ast}(U^{(\ell)}_{b}(\lambda_{j,b}))=1$, that is
$E^{\ast}(I^{(\ell)}_{b}(\lambda_{j,b}))=\widehat{f}(\lambda_{j,b})$.

  2. (ii)

As Lemma 7.2 of Section 7 shows, $V_{2,n}^{+}\stackrel{{\scriptstyle
P}}{{\rightarrow}}V_{2}$. This implies that
$V_{1,n}^{\ast}+V_{2,n}^{+}\stackrel{{\scriptstyle
P}}{{\rightarrow}}V_{1}+V_{2}$, in probability. That is, the pseudo random
variables $M^{+}_{n}$ based on the convolved periodograms generated in Step 4
of the above algorithm, are solely used to estimate the part $V_{2}$ of the
covariance matrix of the distribution of $L_{n}$. As already mentioned, this
part can not be captured by the distribution of
$\sqrt{n}(\widehat{\theta}^{*}_{n}-\widehat{\theta}_{0})$ generated in Step 3
due to the independence of the pseudo periodograms
$I^{\ast}_{n}(\lambda_{j,n})$ across the Fourier frequencies $\lambda_{j,n}$.

  3. (iii)

The statistic $\widetilde{f}$ appearing in Step 4 is itself a nonparametric
estimator of the spectral density $f$ the properties of which have been
investigated in the literature; see Dahlhaus (1985) and the references
therein. However, we use in this step the estimator $\widehat{f}$ for
generating the pseudo periodograms $I_{b}^{(\ell)}$ in order to ensure that
the same spectral density estimator is used here as in Step 2 of the bootstrap
procedure.

######  Remark 3.3.

To understand the motivation behind the displayed equation in Step 6, notice
that as the proof of Theorem 4.1 shows and by Lemma 7.2 of Section 7, we have
that $Z^{\ast}_{n}\stackrel{{\scriptstyle
D}}{{\rightarrow}}{\mathcal{N}}(0,I_{m})$, in probability. Furthermore and by
the same lemma, it holds true that

$(W_{n}^{\ast})^{-1}(V_{1,n}^{\ast}+V_{2,n}^{+})^{1/2}\stackrel{{\scriptstyle
P}}{{\rightarrow}}W^{-1}(V_{1}+V_{2})^{1/2}.$

These results imply that $L_{n}^{\ast}$ defined in Step 6 satisfies
$L_{n}^{\ast}\stackrel{{\scriptstyle
D}}{{\rightarrow}}{\mathcal{N}}(0,W^{-1}(V_{1}+V_{2})W^{-1}),$ in probability,
which coincides with the limiting distribution of the random sequence
$L_{n}=\sqrt{n}(\widehat{\theta}_{n}-\theta_{0})$.

######  Remark 3.4.

The $m\times m$ matrix $W_{n}^{\ast}$ containing the second order partial
derivatives may be difficult to calculate in some situations. In this case an
additional step can be included in the above procedure the aim of which will
be to directly estimate $W_{n}^{\ast}$. To elaborate, let
$\Sigma^{\ast}_{n}={\rm
Var}^{\ast}(\sqrt{n}\big{(}\widehat{\theta}_{n}^{\ast}-\widehat{\theta}_{0})\big{)}$
and recall the definition of $V_{1,n}^{\ast}$ in Step 3. The matrix
$W^{\ast}_{n}$ can then be estimated by

(3.1) $\widehat{W}_{n}^{\ast}=V_{1,n}^{\ast}\big{(}\Sigma_{n}^{\ast}\cdot
V_{1,n}^{\ast}\big{)}^{-1/2}.$

By the property $\Sigma^{\ast}_{n}\stackrel{{\scriptstyle
P}}{{\rightarrow}}W^{-1}V_{1}W^{-1}$, we have by Lemma 7.2 of Section 7, that

$V^{\ast}_{1,n}\big{(}\Sigma_{n}^{\ast}\cdot
V_{1,n}^{\ast}\big{)}^{-1/2}\stackrel{{\scriptstyle
P}}{{\rightarrow}}V_{1}\big{(}W^{-1}V_{1}W^{-1}V_{1}\big{)}^{-1/2}=V_{1}\big{(}W^{-1}V_{1}\big{)}^{-1}=W.$

That is, $\widehat{W}^{\ast}_{n}$ given in (3.1) consistently estimates $W$
and can, therefore, be used to replace $W_{n}^{\ast}$ in the bootstrap
algoritm.

##  4\. Bootstrap Validity

In this section we establish the asymptotic validity of the bootstrap
procedure proposed. Toward this goal we need to impose some conditions on the
dependence structure of the process $\\{X_{t},t\in{\mathbb{Z}}\\}$ as well as
on the smoothness properties of the functions and of the spectral densities
involved. These conditions are summarized in the following assumptions.

Assumption 1:

  1. (i)

The process $\\{X_{t},t\in{\mathbb{Z}}\\}$ has mean zero, is eight-order
stationary, i.e., the joint cumulants up to eighth-order, ${\rm
cum}(X_{t},X_{t+h_{1}},\ldots,X_{t+h_{7}})$, do not depend on $t$ for any
$h_{1},h_{2},\ldots,h_{7}\in{\mathbb{Z}}$. Furthermore,
$\sum_{h\in{\mathbb{Z}}}|h||{\rm cum}(X_{0},X_{h})|<\infty$,
$\inf_{\lambda\in[0,\pi]}f(\lambda)>0$,

$\sum_{h_{1},h_{2},h_{3}\in{\mathbb{Z}}}(|h_{1}|+|h_{2}|+|h_{3}|)|{\rm
cum}(X_{0},X_{h_{1}},X_{h_{2}},X_{h_{3}})|<\infty,$

and

$\sum_{h_{1},\ldots,h_{7}\in{\mathbb{Z}}}|{\rm
cum}(X_{0},X_{h_{1}},\ldots,X_{h_{7}})|<\infty.$

  2. (ii)

The sequence of Whittle estimators
$\\{\sqrt{n}(\widehat{\theta}_{n}-\theta_{0}),n\in{\mathbb{N}}\\}$, satisfies
(2.8).

The above assumptions on the dependence structure of the process
$\\{X_{t},t\in{\mathbb{Z}}\\}$ are rather weak and cover a wide range of
processes considered in the literature; see among others, Rosenblatt (1985),
Doukhan and León (1989), Wu and Shao (2004) for summability properties of
cumulants for processes satisfying different weak dependent conditions. The
requirement of eighth-order stationarity seems inavoidable taking into account
the fact that our derivations include calculations of the variance of time
averaged products of periodograms of subsamples at different frequencies.
Notice that beyond the summability requirements of Assumption 1 and in order
to be as flexible as possible, we do not impose any further conditions on the
moment or on the dependence structure of the underlying process. Instead, we
directly require that the asymptotic normality of the sequence
$\\{\sqrt{n}(\widehat{\theta}_{n}-\theta_{0}),n\in{\mathbb{N}}\\}$, as stated
in Assumption 1(ii), holds true. This also covers a wide range of processes
satisfying a variety of weak dependence conditions; see the discussion before
equation (2.8) in Section 2.

Assumption 2:

  1. (i)

${\mathcal{F}}_{\theta}=\\{f_{\theta},\theta\in\Theta,\inf_{\theta\in\Theta}\inf_{\lambda\in[-\pi,\pi]}f_{\theta}(\lambda)\geq\delta>0\\},$

where $\Theta$ is a compact subset of ${\mathbb{R}}^{m}$ and $\theta_{0}$
defined in (2.5) is unique and belongs to the interior of $\Theta$.

  2. (ii)

$D(\theta,f)$ is twice differentiable in $\theta\in\Theta$ under the integral
sign.

  3. (iii)

$f_{\theta}(\lambda)$ is continuous at any
$(\lambda,\theta)\in[-\pi,\pi]\times\Theta$.

  4. (iv)

The first and second order partial derivatives of
$f^{-1}_{\theta}(\cdot)=1/f_{\theta}(\cdot)$ with respect to $\theta$ are
continuous at any $(\lambda,\theta)\in[-\pi,\pi]\times\Theta$.

  5. (v)

The $m\times m$ matrix
$W=\Big{(}\frac{\displaystyle\partial^{2}}{\displaystyle\partial\theta_{j}\partial\theta_{k}}D(\theta_{0},f)\Big{)}_{j,k=1,\ldots,m}$
is non singular.

Assumption 2 specifies the conditions imposed on the family
${\mathcal{F}}_{\theta}$ of parametric spectral densities considered.
Assumption 2(i) requires that the spectral densities
$f_{\theta}\in{\mathcal{F}}_{\theta}$ are bounded from bellow away from zero
for all frequencies $\lambda\in[-\pi,\pi]$ and for all $\theta\in\Theta$. This
part of the assumption as well as the smoothness properties of $f_{\theta}$
imposed in part (ii) to part (v), are standard and common in the literature;
see Taniguchi (1987) and Dahlhaus and Wefelmeyer (1996).

Our next assumption specifies the required consistency properties of the
nonparametric spectral density estimator $\widehat{f}$ used to generate the
pseudo periodograms $I^{\ast}_{n}$ and $I_{b}^{(\ell)}$. It is a standard
requirement of uniform consistency.

Assumption 3: The nonparametric spectral density estimator $\widehat{f}$
satisfies

$\sup_{\lambda\in[-\pi,\pi]}\big{|}\widehat{f}(\lambda)-f(\lambda)\big{|}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0.$

Our last assumption concerns the rate at which the subsampling size $b$,
involved in the generation of the convolved periodograms of subsamples, is
allowed to increase to infinity with the sample size $n$ in order to ensure
consistency of the estimator $\widehat{V}_{2,n}^{+}$ used.

Assumption 4: $b\rightarrow\infty$ as $n\rightarrow\infty$ such that
$b^{3}/n\rightarrow 0$.

We now state the main result of this paper which establishes consistency of
the bootstrap proposal $L_{n}^{\ast}$ defined in Step 6 of the basic bootstrap
algorithm and used to estimate the distribution of $L_{n}$.

######  Theorem 4.1.

Let Assumptions 1 to 4 be satisfied. Then, as $n\rightarrow\infty$,

$\sup_{x\in{\mathbb{R}}^{m}}\Big{|}P(L_{n}^{\ast}\leq
x\big{|}X_{1},X_{2},\ldots,X_{n})-P(L_{n}\leq x)\Big{|}\rightarrow 0,$

in probability, where
$P(L_{n}^{\ast}\leq\cdot\big{|}X_{1},X_{2},\ldots,X_{n})$ denotes the
distribution function of the random variable $L^{\ast}_{n}$ given the time
series $X_{1},X_{2},\ldots,X_{n}$.

##  5\. Incorporating Whittle Likelihood Modifications

It has been observed that despite their nice properties, Whittle estimators
may behave, in certain, small samples situations, inferior compared to the
exact, time domain, maximum likelihood estimators. More specifically, Whittle
estimators may be biased in small samples due to errors inherited in Whittle’s
frequency domain approximation of the time domain Gaussian maximum likelihood
or in cases where the spectral density of the underlying process contains
(strong) peaks or the periodogram suffers from the well known blurring or
aliasing effects. These drawbacks motivated many researchers to investigate
modifications of the basic Whittle likelihood in order to improve the finite
sample performance of the estimators. Dahlhaus (1988) proposed and
investigated the use of tapered periodograms, while Velasco and Robinson
(2000) combined tapering with differencing the time series before obtaining
Whittle’s estimators. Sykulski et al. (2019) introduced a de-biased Whittle
likelihood to reduce leakage and blurring effects and more recently, Subba Rao
and Yang (2020) proposed the boundary corrected Whittle likelihood and also
combined this with tapering. In this section we will propose modifications,
respectively extensions, of the basic bootstrap algorithm presented in Section
3, which appropriately take into account such modifications of Whittle
estimators.

###  5.1. Tapered Periodograms

Applying a data taper to the time series observed, leads to the replacement of
the periodogram $I_{n}(\lambda_{j,n})$ used in Whittle’s likelihood
approximation (2.1) by a tapered periodogram, denoted by
$I_{n,T}(\lambda_{j,n})$. Let $h_{t,n}=h(t/n)$ be a data taper, that is
$h:{\mathbb{R}}\rightarrow[0,1]$ is a function of bounded variation satisfying
$h(x)=0$ for $x\notin(0,1]$. For any frequency
$\lambda_{j,n}\in{\mathcal{F}}_{n}$, the tapered periodogram is then defined
as

$I_{n,T}(\lambda_{j,n})=\frac{1}{2\pi H_{2,n}}J_{n,T}(\lambda_{j,n})\cdot
J_{n,T}(-\lambda_{j,n}),$

where
$J_{n,T}(\lambda_{j,n})=\sum_{t=1}^{n}h_{t,n}X_{t}\exp\\{-i\lambda_{j,n}t\\}$
is the finite Fourier transform of the tapered time series. Here and for any
$k\in{\mathbb{N}}$, $H_{k,n}=\sum_{t=1}^{n}h_{t,n}^{k}$.

To incorporate data tapering in the bootstrap procedure presented in Section
3, the pseudo periodograms $I_{n}^{\ast}$ and $I_{b}^{(\ell)}$ have to be
replaced by tapered versions, denoted by $I_{n,T}^{\ast}$ and
$I_{b,T}^{(\ell)}$, respectively.

The tapered version $I_{b,T}^{(\ell)}$ of the periodogram $I_{b}^{(\ell)}$,
$\ell=1,2,\ldots,k$, based on subsamples of length $b$ and used in Step 4 and
Step 5, can be generated as

$I_{b,T}^{(\ell)}(\lambda_{j,b})=\widehat{f}(\lambda_{j,b})\cdot
U_{b,T}^{(\ell)}(\lambda_{j,b}),$

where
$U_{b,T}^{(t)}(\lambda_{j,b})=I_{b,T}^{(t)}(\lambda_{j,b})/\widetilde{f}_{T}(\lambda_{j,b})$
for $t=1,2,\ldots,n-b+1$. In the last expression,
$\widetilde{f}_{T}(\lambda_{j,b})=\sum_{t=1}^{n-b+1}I^{(t)}_{b,T}(\lambda_{j,b})/(n-b+1)$
and $I_{b,T}^{(t)}(\lambda_{j,b})=J_{b,T}^{(t)}(\lambda_{j,b})\cdot
J_{b,T}^{(t)}(-\lambda_{j,b})/(2\pi H_{2,b})$ with
$J_{b,T}^{(t)}(\lambda_{j,b})$ the finite Fourier transform of the tapered
subsample $X_{t},X_{t+1},\ldots,X_{t+b-1}$, that is,

$J_{b,T}^{(t)}(\lambda_{j,b})=\sum_{s=1}^{b}h_{s,b}X_{t+s-1}e^{-i\lambda_{j,b}s}.$

The generation of the tapered pseudo periodogram $I_{n,T}^{\ast}$ used to
replace $I^{\ast}_{n}$ in Step 2 and Step 3 is more involved. For this
purpose, the following procedure consisting of Step 2a) to Step 2c) can be
used to generate the tapered pseudo periodograms $I_{n,T}^{\ast}$ and to
replace Step 2 of the basic bootstrap algorithm presented in Section 3 .

  1. Step 2a:

Generate
$\varepsilon_{1}^{\ast},\varepsilon_{2}^{\ast},\ldots,\varepsilon_{n}^{\ast}$
i.i.d., ${\mathcal{N}}(0,1)$ distributed random variables and calculate for
$\lambda_{s,n}\in{\mathcal{F}}_{n}$, the normalized finite Fourier transform,

$Z^{\ast}_{s,n}=\frac{1}{\sqrt{2\pi
n}}\sum_{t=1}^{n}\varepsilon_{t}^{\ast}e^{-it\lambda_{s,n}}.$

  2. Step 2b:

For $t=1,2,\ldots,n$, calculate the pseudo random variables

$X_{t}^{\ast}=\sqrt{\frac{2\pi}{n}}\sum_{\lambda_{s,n}\in{\mathcal{F}}_{n}}\widehat{f}^{1/2}(\lambda_{s,n})Z^{\ast}_{s,n}e^{it\lambda_{s,n}}.$

  3. Step 2c:

For $\lambda_{j,n}\in{\mathcal{F}}_{n}$, calculate the finite Fourier
transform of the tapered pseudo time series
$X_{1}^{\ast},X^{\ast}_{2},\ldots,X_{n}^{\ast}$, that is,

$J_{n,T}^{\ast}(\lambda_{j,n})=\sum_{t=1}^{n}h_{t,n}X_{t}^{\ast}e^{-it\lambda_{j,n}}.$

The tapered pseudo periodogram $I^{\ast}_{n,T}(\lambda_{j,n})$ is then defined
as

$I^{\ast}_{n,T}(\lambda_{j,n})=\frac{1}{2\pi
H_{2,n}}J_{n,T}^{\ast}(\lambda_{j,n})\cdot\overline{J}_{n,T}^{\ast}(\lambda_{j,n}).$

######  Remark 5.1.

Notice that if we set $h_{t,n}\equiv 1$ for $t=1,2,\ldots,n$, in Step 2c) of
the above algorithm, which corresponds to the case of no taper, then we get

$I^{\ast}_{n,T}(\lambda_{j,n})=\widehat{f}(\lambda_{j,n})|Z^{\ast}_{j,n}|^{2}.$

The random variables $|Z_{j,n}^{\ast}|^{2}$, $j=1,2,\ldots,N$, are
independent, which implies that the periodogram ordinates
$I^{\ast}_{n,T}(\lambda_{j,n})$ are independent across the frequencies
$\lambda_{j,n}$, $j=1,2,\ldots,N$. Furthermore and since the
$|Z_{j,n}^{\ast}|^{2}$ also have a standard exponential distribution for every
$j$, we get that if $h_{t,n}\equiv 1$ for all $t=1,2,\ldots,n$, then the
pseudo periodograms $I^{\ast}_{n,T}(\lambda_{j,n})$ generated following Step
2a) to Step 2c) have exactly the same properties as the (non tapered) pseudo
periodograms $I_{n}^{\ast}(\lambda_{j,n})$ generated in Step 2 of the
bootstrap algorithm presented in Section 3.

Now, the tapered pseudo periodograms $I_{b,T}^{(\ell)}$ and $I_{n,T}^{\ast}$
can be used in the bootstrap algorithm to approximate the distribution of
$\sqrt{n}(\widehat{\theta}_{n,T}-\theta_{0})$, where
$\widehat{\theta}_{n,T}=\arg\min_{\theta\in\Theta}D_{n}(\theta,I_{n,T})$ and
$D_{n}(\cdot)$ the function given in (2.1). To elaborate,
$\widehat{\theta}_{n}$ appearing in Step 1 and elsewhere in this algorithm is
replaced by the tapered estimator $\widehat{\theta}_{n,T}$. In Step 2 the
bootstrap periodogram $I^{\ast}_{n}$ is replaced by $I^{\ast}_{n,T}$ and
consequently $\widehat{\theta}_{n}^{\ast}$ in Step 3 by

$\widehat{\theta}_{n,T}^{\ast}=\arg\min_{\theta\in\Theta}D_{n}(\theta,I^{\ast}_{n,T}).$

Thus the bootstrap approximation of
$\sqrt{n}(\widehat{\theta}_{n,T}-\theta_{0})$ in the same step is given by
$\widetilde{L}_{n}^{\ast}=\sqrt{n}(\widehat{\theta}^{\ast}_{n,T}-\widehat{\theta}_{0})$.
Furthermore, the tapered pseudo periodogram $I^{\ast}_{n,T}$ is used in the
expression of the vector $M_{n}^{\ast}$ in Step 3, while the matrix
$W^{\ast}_{n}$ in the same step is calculated using
$D_{n}(\theta,I^{\ast}_{n,T})$, that is $W_{n}^{\ast}$ is replaced by

$W^{\ast}_{n,T}=\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}D_{n}(\theta,I^{*}_{n,T})\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)}.$

Finally, in Step 4, the periodogram of the random subsamples $I_{b}^{(\ell)}$
in the expression for $M_{n}^{+}$ as well as the periodogram $I_{b}^{(t)}$ in
the same step and in Step 5, are replaced by their tapered versions
$I_{b,T}^{(\ell)}$ and $I_{b,T}^{(t)}$, respectively.

###  5.2. Debiased Whittle Likelihood

Debiasing the Whittle likelihood has been proposed by Sykulski et al. (2019).
The basic idea is to replace the parametric spectral density $f_{\theta}$
appearing in $D_{n}(\theta,I_{n})$ by a smoothed version which equals the
expectation of the periodogram $I_{n}(\lambda)$ under the assumption that
$f=f_{\theta}$. More specifically, the objective function considered by this
modification of Whittle’s approximation of the quasi Gaussian likelihood is
given by

(5.1)
$D_{n}^{(db)}(\theta,I_{n})=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{\\{}\log\overline{f}_{\theta}(\lambda_{j,n})+\frac{\displaystyle
I_{n}(\lambda_{j,n})}{\displaystyle\overline{f}_{\theta}(\lambda_{j,n})}\Big{\\}},$

where

$\overline{f}_{\theta}(\lambda_{j,n})=\int_{-\pi}^{\pi}K_{n}(w-\lambda_{j,n})f_{\theta}(\omega)d\omega$

and $K_{n}(\cdot)$ is the Fejer-kernel,

$K_{n}(x)={\bf 1}_{\\{x=0\\}}n/2\pi+{\bf 1}_{\\{x\neq
0\\}}\sin^{2}(nx/2)\big{/}(2\pi n\sin^{2}(x/2)\big{)}.$

Notice that if $f=f_{\theta}$, then ${\rm
E}(I_{n}(\lambda_{j,n}))=\overline{f}_{\theta}(\lambda_{j,n})$; see for
instance, Rosenblatt (1963). That is, in this case, the periodogram $I_{n}$
apearing in $D^{(db)}(\theta,I_{n})$ is an unbiased estimator of
$\overline{f}_{\theta}$, which justifies the name given to this modification.

In the following steps we summarize the modifications needed in order to adapt
the basic bootstrap algorithm presented in Section 3 to imitate the random
properties of the estimator $\widetilde{\theta}_{n}$.

  1. Step I: Calculate $\widetilde{\theta}_{n}=\arg\min_{\theta\in\Theta}D_{n}^{(db)}(\theta,I_{n})$.

  2. Step II: The same as Step 2 of the bootstrap algorithm in Section 3.

  3. Step III: Calculate $\widetilde{\theta}^{\ast}_{n}=\arg\min_{\theta\in\Theta}\,D^{(db)}_{n}(\theta,I^{\ast}_{n})$, $\widetilde{\theta}_{0}=\arg\min_{\theta\in\Theta}\,D^{(db)}_{n}(\theta,\widehat{f})$ and $\widetilde{V}_{1,n}^{\ast}=Var^{\ast}(\widetilde{M}_{n}^{\ast})$, where

$\widetilde{M}^{\ast}_{n}=\frac{2\pi}{\sqrt{n}}\sum_{j\in{\mathcal{G}}(n)}\overline{g}_{\widetilde{\theta}_{0}}(\lambda_{j,n})\big{(}I_{n}^{\ast}(\lambda_{j,n})-\widehat{f}(\lambda_{j,n})\big{)}$

and $\overline{g}_{\widetilde{\theta}_{0}}(\lambda)$ is the $m$-dimensional
vector

$\overline{g}_{\widehat{\theta}_{0}}(\lambda)=\Big{(}\overline{g}_{j,\widehat{\theta}_{0}}(\lambda)=-\frac{1}{2\pi}\frac{\partial}{\partial\theta_{j}}\overline{f}^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\widetilde{\theta}_{0}},\
j=1,2,\ldots,m\Big{)}^{\top}.$

Calculate

$\widetilde{W}_{n}^{\ast}=\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}D^{(db)}_{n}(\theta,I^{*}_{n})\Big{|}_{\theta=\widetilde{\theta}_{0}}\Big{)}$

and the pseudo random vector $\widetilde{Z}_{n}^{\ast}$ defined by

$\widetilde{Z}^{\ast}_{n}=\big{(}\widetilde{V}^{\ast}_{1,n}\big{)}^{-1/2}\widetilde{W}^{\ast}_{n}\sqrt{n}(\widetilde{\theta}^{*}_{n}-\widetilde{\theta}_{0}).$

  4. Step IV: The same as Step 4 of the bootstrap algorithm in Section 3 but by replacing $M_{n}^{+}$ by

$\widetilde{M}^{+}_{n}=\sqrt{kb}\frac{1}{k}\sum_{l=1}^{k}\frac{2\pi}{b}\sum_{j\in{\mathcal{G}}(b)}\overline{g}_{\widetilde{\theta}_{0}}(\lambda_{j,b})\big{(}I_{b}^{(\ell)}(\lambda_{j,b})-\widehat{f}(\lambda_{j,b})\big{)}.$

  5. Step V: Calculate the $m\times m$ matrix $\widetilde{V}_{2,n}^{+}=\widetilde{\Sigma}^{+}_{n}-\widetilde{C}_{n}^{+}$, where $\widetilde{\Sigma}^{+}_{n}={\rm Var}^{\ast}(\widetilde{M}_{n}^{+})$, $\widetilde{C}_{n}^{+}=\big{(}\widetilde{c}_{n}^{(r,s)}\big{)}_{r,s=1,2,\ldots,m}$ and the elements $\widetilde{c}_{n}^{(r,s)}$ given by

$\displaystyle\widetilde{c}_{n}^{(r,s)}$
$\displaystyle=\frac{8\pi^{2}}{b}\sum_{j\in{\mathcal{G}}(b)}\overline{g}_{r,\widetilde{\theta}_{0}}(\lambda_{j,b})\overline{g}_{s,\widetilde{\theta}_{0}}(\lambda_{j,b})\widehat{f}(\lambda_{j,b})^{2}$
$\displaystyle\ \ \ \ \
\times\Big{(}\frac{1}{n-b+1}\sum_{t=1}^{n-b+1}\frac{I_{b}^{(t)}(\lambda_{j,b})^{2}}{\widetilde{f}_{b}(\lambda_{j,b})^{2}}-1\Big{)}.$

  6. Step VI: Appoximate the distribution of $L_{n}=\sqrt{n}(\widetilde{\theta}_{n}-\theta_{0})$ by the distribution of

$\widetilde{L}_{n}^{\ast}=\big{(}\widetilde{W}_{n}^{\ast}\big{)}^{-1}\big{(}\widetilde{V}_{1,n}^{\ast}+\widetilde{V}_{2,n}^{+}\big{)}^{1/2}\cdot\widetilde{Z}^{\ast}_{n}.$

###  5.3. Boundary Corrected Whittle Likelihood

In order to reduce the bias of Whittle estimators caused by boundary effects,
Subba Rao and Yang (2021) proposed the so called, boundary corrected Whittle
likelihood. The objective function to be minimized in this case is given by

(5.2)
$D_{n}^{(bc)}(\theta,\widetilde{I}_{n})=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{\\{}\log
f_{\theta}(\lambda_{j,n})+\frac{\displaystyle\widetilde{I}_{n}(\lambda_{j,n})}{\displaystyle
f_{\theta}(\lambda_{j,n})}\Big{\\}},$

where $\widetilde{I}_{n}(\lambda_{j,n})=(2\pi
n)^{-1}\widetilde{J}_{n}(\lambda_{j,n})\overline{J}_{n}(\lambda_{j,n})$. Here,
$J_{n}(\lambda_{j,n})=\sum_{t=1}^{n}X_{t}e^{-i\lambda_{j,n}t}$, is the finite
Fourier transform of the time series $X_{1},X_{2},\ldots,X_{n}$ while
$\widetilde{J}_{n}(\lambda_{j,n})$ is the boundary extended finite Fourier
transform given by
$\widetilde{J}_{n}(\lambda_{j,n})=J_{n}(\lambda_{j,n})+\widehat{J}_{n}(\lambda_{j,n})$.
The boundary extended finite Fourier transform
$\widetilde{J}_{n}(\lambda_{j,n})$ is obtained by calculating the finite
Fourier transform of the out of sample extended time series

(5.3)
$\ldots,\widehat{X}_{-1},\widehat{X}_{0},X_{1},X_{2},\ldots,X_{n},\widehat{X}_{n+1},\widehat{X}_{n+2},\ldots,$

where the pseudo observations $\widehat{X}_{t}$ are the best linear predictors
of the corresponding (not observed) values $X_{t}$ based on an AR(p) model. To
elaborate, let $(\widehat{\phi}_{s,p},s=1,2,\ldots,p)^{\top}$ be the vector of
Yule-Walker estimators obtained by fitting an AR(p) model to the time series
$X_{1}X_{2},\dots,X_{n}$. Then
$\widehat{X}_{n+s}=\sum_{j=1}^{p}\widehat{\phi}_{j,p}\widetilde{X}_{n+s-j}$
with $\widetilde{X}_{t}=X_{t}$ if $t\in\\{1,2,\ldots,n\\}$ and
$\widetilde{X}_{t}=\widehat{X}_{t}$ if $t>n$. An analogue expression yields
for $\widehat{X}_{-s}$, $s\leq 0$. This leads to the expression
$\widetilde{J}_{n}(\lambda_{j,n})=J_{n}(\lambda_{j,n})+\widehat{J}_{n}(\lambda_{j,n})$,
where the “extension” part $\widehat{J}_{n}(\lambda_{j,n})$ can be written as

$\displaystyle\widehat{J}_{n}(\lambda_{j,n})$
$\displaystyle=\frac{1}{\widehat{\phi}_{p}(\lambda_{j,n})}\sum_{\ell=1}^{p}X_{\ell}\sum_{s=0}^{p-\ell}\widehat{\phi}_{\ell+s,p}e^{-is\lambda_{j,n}}$
$\displaystyle\ \ \ \ \
+e^{in\lambda_{j,n}}\frac{1}{\overline{\widehat{\phi}_{p}(\lambda_{j,n})}}\sum_{\ell=1}^{p}X_{n+1-\ell}\sum_{s=0}^{p-\ell}\widehat{\phi}_{\ell+s,p}e^{i(s+1)\lambda_{j,n}},$

with
$\widehat{\phi}_{p}(\lambda)=1-\sum_{s=1}^{p}\widehat{\phi}_{s,p}e^{-is\lambda}$;
see Subba Rao and Yang (2021).

In the following we only describe how to modify the first part of the basic
bootstrap algorithm, that is Step 1 to Step 3, in order to get replicates of
the boundary corrected Whittle estimators
$\breve{\theta}_{n}=\arg\min_{\theta\in\Theta}\,D^{(bc)}_{n}(\theta,\widetilde{I}_{n})$.
The modifications needed for the second part of the basic bootstrap procedure
which uses convolved periodograms of subsamples (Step 4 to Step 5), easily
follow from those presented for the first part.

  1. Step 1${}^{{}^{\prime}}$:

Calculate Whittle’s estimator
$\breve{\theta}_{n}=\arg\min_{\theta\in\Theta}\,D^{(bc)}_{n}(\theta,\widetilde{I}_{n})$.

  2. Step 2${}^{{}^{\prime}}$:

For $t=1,2,\ldots,n$, calculate the pseudo random variables

$X_{t}^{\ast}=\sqrt{\frac{2\pi}{n}}\sum_{\lambda_{s,n}\in{\mathcal{F}}_{n}}\widehat{f}^{1/2}(\lambda_{s,n})Z^{\ast}_{s,n}e^{it\lambda_{s,n}},$

where

$Z^{\ast}_{s,n}=\frac{1}{\sqrt{2\pi
n}}\sum_{t=1}^{n}\varepsilon_{t}^{\ast}e^{-it\lambda_{s,n}}$

and
$\varepsilon_{1}^{\ast},\varepsilon_{2}^{\ast},\ldots,\varepsilon_{n}^{\ast}$
are i.i.d., ${\mathcal{N}}(0,1)$ distributed random variables. Calculate

$\widetilde{J}_{n}^{\ast}(\lambda_{j,n})=J^{\ast}_{n}(\lambda_{j,n})+\widehat{J}_{n}^{\ast}(\lambda_{j,n}),$

where
$J^{\ast}_{n}(\lambda_{j,n})=\sum_{t=1}^{n}X^{\ast}_{t}\exp\\{-i\lambda_{j,n}t\\}$
and

$\displaystyle\widehat{J}^{\ast}_{n}(\lambda_{j,n})$
$\displaystyle=\frac{1}{\widehat{\phi}^{\ast}_{p}(\lambda_{j,n})}\sum_{\ell=1}^{p}X^{\ast}_{\ell}\sum_{s=0}^{p-\ell}\widehat{\phi}^{\ast}_{\ell+s,p}e^{-is\lambda_{j,n}}$
$\displaystyle\ \ \ \ \
+e^{in\lambda_{j,n}}\frac{1}{\overline{\phi^{\ast}_{p}(\lambda_{j,})}}\sum_{\ell=1}^{p}X^{\ast}_{n+1-\ell}\sum_{s=0}^{p-\ell}\widehat{\phi}^{\ast}_{\ell+s,p}e^{i(s+1)\lambda_{j,n}}.$

In the above expression,
$\widehat{\phi}^{\ast}_{p}(\lambda)=1-\sum_{s=1}^{p}\widehat{\phi}^{\ast}_{s,p}e^{-is\lambda}$
and $(\widehat{\phi}^{\ast}_{s,p},s=1,2,\ldots,p)^{\top}$ is the vector of
Yule-Walker estimators obtained by fitting an AR(p) model to the pseudo time
series $X_{1}^{\ast},X_{2}^{\ast},\dots,X_{n}^{\ast}$. Define,

$\widetilde{I}_{n}^{\ast}(\lambda_{j,n})=\frac{1}{2\pi
n}\widetilde{J}_{n}^{\ast}(\lambda_{j,n})\overline{J_{n}^{\ast}}(\lambda_{j,n}).$

  3. Step 3${}^{{}^{\prime}}$:

Let

$\breve{\theta}^{\ast}_{n}=\arg\min_{\theta\in\Theta}\,D_{n}(\theta,\widetilde{I}^{\ast}_{n}).$

Define

$\widetilde{V}_{1,n}^{\ast}=Var^{\ast}(\widetilde{M}_{n}^{\ast})\ \
\mbox{and}\ \
\widetilde{W}_{n}^{\ast}=\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}D_{n}(\theta,\widetilde{I}^{*}_{n})\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)},$

where

$\widetilde{M}^{\ast}_{n}=\frac{2\pi}{\sqrt{n}}\sum_{j\in{\mathcal{G}}(n)}g_{\widehat{\theta}_{0}}(\lambda_{j,n})\big{(}\widetilde{I}_{n}^{\ast}(\lambda_{j,n})-\widehat{f}(\lambda_{j,n})\big{)},$

and $g_{\widehat{\theta}_{0}}(\lambda)$ the $m$-dimensional vector

$g_{\widehat{\theta}_{0}}(\lambda)=\Big{(}g_{j,\widehat{\theta}_{0}}(\lambda)=-\frac{1}{2\pi}\frac{\partial}{\partial\theta_{j}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\widehat{\theta}_{0}},\
j=1,2,\ldots,m\Big{)}^{\top}.$

Calculate the pseudo random vector $\widetilde{Z}_{n}^{\ast}$ defined by

$\widetilde{Z}^{\ast}_{n}=\big{(}\widetilde{V}^{\ast}_{1,n}\big{)}^{-1/2}\widetilde{W}^{\ast}_{n}\sqrt{n}(\breve{\theta}^{*}_{n}-\widehat{\theta}_{0}).$

##  6\. Practical Implementation and Numerical Results

###  6.1. Choice of bootstrap parameters

To implement our procedure we need to choose two parameters, the bandwidth
involved in obtaining the spectral density estimator $\widehat{f}$ and the
blocksize $b$, used in the convolved part of our procedure. Regarding
$\widehat{f}$ we use a kernel type estimator obtained by smoothing the
periodogram by means of the Bartlett-Priestley kernel. The bandwidth used for
this estimator is obtained using a frequency domain cross validation
procedure; see Beltrão and Bloomfield (1987). For the choice of the
subsampling parameter $b$, observe first that, as we will see in the next
section, our simulation results seem not to be very sensitive with respect to
the choice of this parameter, provided $b$ not chosen too large with respect
to $n$. Based on this empirical observation we use the following practical
rule to select this parameter: $b=4\cdot n^{0.25}$. This rule delivers a
subsampling size which is not too large and at the same time also satisfies
the conditions of Assumption 4. However, in order to see the sensitivity of
the bootstrap approximations obtained with respect to the choice of the
subsampling parameter, we present in Section 6.2 results for a large range of
values of $b$. The aforementioned rule for choosing $b$ has been applied for
obtaining the numerical results presented in Section 6.3 for analyzing the
real-life data example.

###  6.2. Simulations

Let ${\mathcal{F}}_{\theta}$ be the family containing the spectral densities
of the first order autoregressive processes given by

$f_{\theta}(\lambda)=\frac{\sigma^{2}}{2\pi}(1+a^{2}-2a\cos(\lambda))^{-1},$

where $\theta=(\sigma^{2},a)\in(0,+\infty)\times(-1,1)$. To demonstrate the
advances of the bootstrap procedure proposed in this paper, we concentrate in
the following on the (standard) Whittle estimator of $a$ given by

$\widehat{a}_{n}=\sum_{j\in{\mathcal{G}}(n)}I_{n}(\lambda_{j,n})\cos(\lambda_{j,n})\big{/}\sum_{j\in{\mathcal{G}}(n)}I_{n}(\lambda_{j,n}).$

Observe that
$a_{0}=\int_{-\pi}^{\pi}\cos(\lambda)f_{\theta_{0}}(\lambda)d\lambda/\int_{-\pi}^{\pi}f_{\theta_{0}}(\lambda)d\lambda=\rho(1)$,
where $\rho(1)$ denotes the first order autocorrelation of the “best fitting”
AR(1) process; see Section 2. As we have seen, the distribution of
$\sqrt{n}(\widehat{a}_{n}-a_{0})$ depends on the dependence properties of the
underlying process $\\{X_{t};t\in{\mathbb{Z}}\\}$. In order to demonstrate the
finite sample behavior and the capabilities of the bootstrap procedure
proposed to approximate the distribution of $\sqrt{n}(\widehat{a}_{n}-a_{0})$
for a variety of situations, we consider time series
$X_{1},X_{2},\ldots,X_{n}$ stemming from the following three processes:

  1. Model I: $X_{t}=0.8X_{t-1}+\varepsilon_{t}$, and i.i.d. innovations $\varepsilon_{t}\sim{\mathcal{N}}(0,1)$.

  2. Model II: $X_{t}=0.75X_{t-1}+0.6X_{t-1}\cdot\varepsilon_{t-1}+\varepsilon_{t}$, and i.i.d. innovations $\varepsilon_{t}\sim{\mathcal{L}aplace}(0,0.1)$.

  3. Model III: $X_{t}=\left\\{\begin{array}[]{lll}-0.3X_{t-1}+\varepsilon_{t}&&\mbox{if $X_{t-1}\leq 0$}\\\ &&\\\ 0.8X_{t-1}+\varepsilon_{t}&&\mbox{if $X_{t-1}>0$,}\end{array}\right.$
and i.i.d. innovations $\varepsilon_{t}\sim{\mathcal{L}aplace}(0,0.1)$.

Model I is a Gaussian AR(1) model. Model II and Model III are nonlinear and
have been considered in Fan and Yao (2005) and Auestad and Tjøstheim (1990),
respectively. They have been modified so that they are driven by Laplace
(double exponential) distributed innovations with mean zero and parameter
$0.1$. Notice that fitting a linear AR(1) model to time series stemming from
the nonlinear Models II and III, resamples a situation of model
misspecification.

Two sample sizes, $n=50$ and $n=1000$ are considered in order to investigate
the small sample behavior of the bootstrap procedure proposed as well as its
consistency behavior when the length of the time series becomes large. In
order to see the effects of the convolved subsampling step implemented in Step
4 through generating the pseudo random variables $M^{+}_{n}$, we also present
results for the bootstrap approximation of $\sqrt{n}(\widehat{a}_{n}-a_{0})$
using the multiplicative periodogram bootstrap only, that is, the pseudo
random variable
$L_{n,MB}^{\ast}=\sqrt{n}(\widehat{a}^{\ast}_{n}-\widehat{a}_{0})$, generated
in Step 3 of the bootstrap algorithm of Section 3. $B=1000$ bootstrap
replications have been used in each run and the $d_{1}$-distances, between the
exact distribution and the two bootstrap approximations, that is the
multiplicative periodogram bootstrap $L_{n,MB}^{\ast}$ and the hybrid
periodogram bootstrap $L_{n}^{\ast}$ as generated in Step 6, have been
calculated. Notice that for $F$ and $G$ distribution functions,
$d_{1}=\int_{0}^{1}|F^{-1}(u)-G^{-1}(u)|du$. To estimate the exact
distribution of $\sqrt{n}(\widehat{a}_{n}-a_{0})$, $R=10,000$ replications
have been used. Figure 1 shows averages of the $d_{1}$ distances calculated
over $500$ repetitions for each of the three different time series models and
for each of the two sample sizes considered.

Figure 1. Average $d_{1}$-distances between the exact and the bootstrap
distribution of the Whittle estimator $\sqrt{n}(\widehat{a}_{n}-a_{0})$ for
various block sizes $b$. Left column n=50, right column, n=1,000. First row
Model I, second row Model II and third row Model III. The crosses denote the
$d_{1}$-distance of the multiplicative periodogram bootstrap estimation
$L^{\ast}_{n,MB}$ and the circles of the hybrid periodogram bootstrap
estimation $L^{\ast}_{n}$. The dashed lines with the plus symbol in the first
row refer to the average $d_{1}$-distance of the asymptotic Gaussian
approximation.

As it is seen from Figure 1, for the case of the Gaussian AR(1) model and for
both sample sizes considered, the bootstrap approximations $L_{n}^{\ast}$ and
$L_{n,MB}^{\ast}$ behave very similar and the corresponding $d_{1}$ distances
are very close to each other. Recall that in this case the distribution of
$\sqrt{n}(\widehat{a}_{n}-a_{0})$ only depends on the second order
characteristics of the underlying AR(1) process and therefore, the
multiplicative bootstrap estimation $L_{n,MB}^{\ast}$ also provides a
consistent estimation of the distribution of $L_{n}$. As it is also seen, for
both sample sizes considered and for the case of the Gaussian AR(1) model,
both frequency domain bootstrap procedures outperform the asymptotic Gaussian
approximation. This is due to the skewness of the distribution of
$\sqrt{n}(\widehat{a}-a_{0})$ which does not vanish even for $n=1000$
observations. In the case of the nonlinear models considered, that is for
Model II and Model III, the behavior of the multiplicative and of the hybrid
bootstrap is very different. Recall that in these cases, the multiplicative
bootstrap fails to appropriately capture the fourth order characteristics of
the underlying nonlinear processes that affect the distribution of
$\sqrt{n}(\widehat{a}_{n}-a_{0})$. This leads to a larger $d_{1}$ distance of
the bootstrap estimator $L^{\ast}_{n,MB}$ compared to the hybrid bootstrap
estimator $L^{\ast}_{n}$. The hybrid bootstrap captures these characteristics
and performs much better leading to an overall smaller $d_{1}$-distance. For
the case $n=50$, this is true for all block sizes $b$ and for all models
considered. Only for Model III and for the sample size of $n=50$, the behavior
of the hybrid periodogram bootstrap gets closer to that of the multiplicative
bootstrap when the block size $b$ becomes too large. Finally, for the sample
size of $n=1,000$ observations, the advantages of the hybrid bootstrap
procedure are clearly seen in the corresponding exhibits of Figure 1 for all
models.

###  6.3. Periodicity of Annual Sunspot Data

We consider the yearly mean total sunspot numbers from 1700 to 2020 available
at www.sidc.be/silso/datafiles. A plot of the corresponding time series
consisting of $n=321$ observations, is shown in Figure 2(a). Our aim is to
estimate the main periodicity of this time series and to infer properties of
the estimator used by applying the frequency domain bootstrap procedure
proposed in this paper. To make things precise, suppose that the stochastic
process generating the observed yearly mean sunspot data possesses a spectral
density $f$ and assume that a unique frequency $\lambda_{\max}\in(0,\pi)$
exists such that $\lambda_{\max}=\arg\min_{\lambda}f(\lambda)$. We are
interested in the main periodicity of the yearly sunspot numbers defined as
the parameter $P_{X}=2\pi/\lambda_{\max}$.

One approach to estimate this parameter is to use the class of linear $AR(p)$
process to get an estimate of the frequency $\lambda_{\max}$. To elaborate,
suppose that an AR(p) model is fitted to the time series of sunspot numbers
and that $\widehat{\lambda}_{\max,AR}$ is the (unique) frequency in $(0,\pi)$
defined by
$\widehat{\lambda}_{\max,AR}=\arg\min_{\lambda}\widehat{f}_{AR}(\lambda)$,
where $\widehat{f}_{AR}$ is the spectral density of the estimated AR(p) model.
The estimator of $P_{X}$ obtained following this approach is then defined as
$\widehat{P}_{X}=2\pi/\widehat{\lambda}_{\max,AR}$. Observe that consistency
of the estimator $\widehat{P}_{X}$ only requires that
$\widehat{\lambda}_{\max,AR}\stackrel{{\scriptstyle
P}}{{\rightarrow}}\lambda_{\max}$, as $n\rightarrow\infty$. This can be
achieved if the spectral density, say $f_{L}$, to which $\widehat{f}_{AR}$
uniformly converges, that is,
$\sup_{\lambda\in[0,\pi]}|\widehat{f}_{AR}(\lambda)-f_{L}(\lambda)|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$, satisfies $\lambda_{\max,L}=\lambda_{\max}$, where
$\lambda_{\max,L}=\arg\min_{\lambda}f_{L}(\lambda)$. Hence consistency of
$\widehat{\lambda}_{\max,AR}$ only requires that the limiting spectral density
$f_{L}$ to which $\widehat{f}_{AR}$ uniformly converges in probability, has
its largest peak at the same frequency as the spectral density $f$ of
interest. One way to achieve this, is to allow for the order $p$ of the AR
model fitted, to increase to infinity at some appropriate rate as $n$
increases to infinity. Under certain conditions it can then be shown that
$f_{L}=f$, i.e.,
$\sup_{\lambda\in[0,\pi]}|\widehat{f}_{AR}(\lambda)-f(\lambda)|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$ and that the corresponding estimator
$\widehat{P}_{X}=2\pi/\widehat{\lambda}_{\max,AR}$ achieves the rate
$\widehat{P}_{X}=P_{X}+O_{P}(p^{3/2}/n^{1/2})$; see Newton and Pagano (1983).

However, an alternative way to consistently estimate $P_{X}$ is the following.
Suppose that there exists a finite order AR(p) model possessing a spectral
density $f_{AR}$ such that $\lambda_{\max,AR}=\lambda_{\max}$, where
$\lambda_{\max,AR}=\arg\min_{\lambda}f_{AR}(\lambda)$ and $f_{AR}$ denotes the
spectral density of the AR(p) process. Then
$\sup_{\lambda\in[0,\pi]}|\widehat{f}_{AR}(\lambda)-f_{AR}(\lambda)|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$ implies
$\widehat{P}_{X}=2\pi/\widehat{\lambda}_{\max,AR}\stackrel{{\scriptstyle
P}}{{\rightarrow}}P_{X}$ and this estimator converges at the parametric rate
$\widehat{P}_{X}=P_{X}+O_{P}(1/n^{1/2})$. Notice that consistency of the
described approach, does not rely on the assumption that the AR(p) model
correctly describes the entire stochastic structure of the process generating
the sunspot data. Not even the entire autocovariance structure of the sunspot
time series has to appropriately be captured by the AR(p) model. What is
solely required is that the spectral density $f_{AR}$ of the AR(p) process has
its main peak at the same frequency as the spectral density $f$ of the
stochastic process generating the sunspot time series. The AR(p) model is then
solely used as a vehicle to construct an estimator of the main periodicity
$P_{X}$. Moreover, in conjunction with the frequency domain bootstrap, this
approach also allows for the investigation of the sampling properties of the
estimator $\widehat{P}_{X}$ and for quantifying the uncertainty associated
with estimating the parameter $P_{X}$ of interest.

For the yearly sunspot time series, selecting an AR(p) model using Akaike’s
Information Criterion (AIC), leads to an AR(9) model. However, and as already
mentioned, since we are not interested in parametrizing the entire
autocovariance structure of the yearly sunspot numbers but solely in
consistently estimating the frequency $\lambda_{\max}$, an AR(2) model,
$X_{t}=a_{1}X_{t-1}+a_{2}X_{t-2}+\varepsilon_{t}$, also can be used for this
purpose. Figure 2(b) demonstrates this by showing the periodogram of the
sunspot time series together with the spectral densities of the fitted AR(2)
and AR(9) models. Notice that the periodogram of this time series takes its
maximum value at the Fourier frequency $\lambda_{j,n}=0.090625$, which
corresponds to a periodicity of $11.034$ years. Fitting an AR(2) model also
leads to the estimate $\widehat{\lambda}_{\max,AR}=0.090625$, which
corresponds to the same estimate of the main periodicity
$\widehat{P}_{X}=11.034$. Fitting the AR(9) model leads to the estimates
$0.09375$ for the frequency $\lambda_{\max}$ and $10.667$ years for the main
periodicity $P_{X}$. We, therefore, proceed by using the more parsimonious
AR(2) model for our analysis. In particular, we apply the frequency domain
bootstrap procedure proposed in this paper to generate replicates of the the
estimated parameters
$\widehat{\theta}_{n}=(\widehat{\sigma}^{2},\widehat{a}_{1},\widehat{a}_{2})^{\top}$
of the AR(2) model. Clearly and since $f=f_{AR}$ is not a reasonable
assumption in our context, we are in the setting of model misspecification.
Using the bootstrap replicates of the estimated parameters, we get replicates
of the estimated spectral density of the AR(2) model, say
$\widehat{f}_{AR}^{\ast}$. Bootstrap replicates of
$\widehat{\lambda}_{\max,AR}$ can then be obtained as
$\widehat{\lambda}_{\max,AR}^{\ast}=\arg\min_{\lambda}\widehat{f}^{\ast}_{AR}(\lambda)$
which lead to bootstrap replicates
$\widehat{P}^{\ast}_{X}=2\pi/\widehat{\lambda}_{\max,AR}^{\ast}$ of the
estimator of the main periodicity
$\widehat{P}_{X}=2\pi/\widehat{\lambda}_{\max,AR}$ used. By repeating these
steps a large number, say $B$, of times, bootstrap estimators of the
distribution of $\widehat{\lambda}_{\max,AR}$ and, conseqeuntly, of
$\widehat{P}_{X}$, are obtained. Figure 2(c) presents a histogram of $B=1,000$
bootstrap replicates of $\widehat{\lambda}_{\max,AR}^{\ast}$ and Figure 2(d)
of the corresponding estimates
$\widehat{P}^{\ast}_{X}=2\pi/\widehat{\lambda}_{\max,AR}^{\ast}$ obtained by
using a grid of $500$ equidistant frequencies in the interval $(0,\pi)$. The
corresponding $95\%$ confidence interval for the main periodicity of the
sunspot time series, based on bootstrap percentages, is then given by
$[9.90,12.98]$.

Figure 2. Reading clockwise from top to bottom: (a) Time series of yearly mean
sunspot numbers. (b) Periodogram of the time series with estimated spectral
densities of the AR(2) model (solid line) and of the AR(9) model (dashed
line), log-scale. (c) Histogram of $B=1000$ replications of
$\widehat{\lambda}^{\ast}_{\max,AR}$ and (d) histogram of the corresponding
replications of the estimated main periodicity$\widehat{P}^{\ast}_{X}$, both
using the AR(2) model. The estimated values $\widehat{\lambda}_{\max,AR}$ and
$\widehat{P}_{X}$ are indicated in (c) and (d) by vertical dashed lines.

##  7\. Auxiliary Results and Proofs

To simplify notation we write $D_{n}^{(j)}(\mathring{\theta},I_{n}^{\ast})$,
$\mathring{\theta}\in\Theta$, for the $j$th derivative of
$D_{n}(\theta,I^{\ast}_{n})$ with respect to $\theta$ evaluated at
$\theta=\mathring{\theta}$. We also write $\|\cdot\|$ for the Euclidean norm
in ${\mathbb{R}}^{m}$ and $\|A\|_{F}$ for the Frobenius norm of a matrix
$A\in{\mathbb{R}}^{m\times m}$. Furthermore, at different places we will use
the expansion

(7.1)
$D^{(1)}_{n}(\theta,I^{\ast}_{n})=D^{(1)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})+D^{(2)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})(\theta-\widehat{\theta}_{0})+R_{n}(\theta)(\theta-\widehat{\theta}_{0}),$

where
$R_{n}(\theta)=D^{(2)}_{n}(\theta^{+},I^{\ast}_{n})-D^{(2)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})$
for some $\theta^{+}\in\Theta$ such that
$\|\theta^{+}-\widehat{\theta}_{0}\|\leq\|\theta-\widehat{\theta}_{0}\|$.

We first establish the following two useful lemmas.

######  Lemma 7.1.

If Assumption 2 and Assumption 3 are satisfied, then the following assertions
hold true.

  1. (i)

$\widehat{\theta}_{0}\stackrel{{\scriptstyle P}}{{\rightarrow}}\theta_{0}$.

  2. (ii)

$\big{\|}\widehat{\theta}_{n}^{\ast}-\widehat{\theta}_{0}\big{\|}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$, in probability.

###### Proof.

Consider (i). Since $\widehat{\theta}_{0}\in\Theta$ and
$\Theta\subset{\mathbb{R}}^{m}$ is compact,
$\\{\widehat{\theta}_{0},n\in{\mathbb{N}}\\}$ is a bounded sequence. By the
continuity of $D(\cdot,f)$ as a function on $\theta\in\Theta$, it suffices to
show that

$D(\widehat{\theta}_{0},f)\stackrel{{\scriptstyle
P}}{{\rightarrow}}D(\theta_{0},f).$

It yields

$\displaystyle|D(\widehat{\theta}_{0},f)-D(\theta_{0},f)|$
$\displaystyle\leq|D(\widehat{\theta}_{0},f)-D_{n}(\widehat{\theta}_{0},f)|+|D_{n}(\widehat{\theta}_{0},f)-D_{n}(\widehat{\theta}_{0},\widehat{f})|$
$\displaystyle\ \ \ \
+|D_{n}(\widehat{\theta}_{0},\widehat{f})-D(\widetilde{\theta}_{0},\widehat{f})|+|D(\widetilde{\theta}_{0},\widehat{f})-D(\theta_{0},f)|$
$\displaystyle=\sum_{j=1}^{4}D_{j,n},$

where $\widetilde{\theta}_{0}=\mbox{argmin}_{\theta}D(\theta,\widehat{f})$ and
with an obvious notation for $D_{j,n}$, $j=1,\ldots,4$. We show that
$D_{j,n}\stackrel{{\scriptstyle P}}{{\rightarrow}}0$, as $n\rightarrow\infty$,
for $j=1,\ldots,4$. We have

$\displaystyle D_{1,n}$
$\displaystyle\leq\sup_{\theta\in\Theta}\big{|}\frac{1}{2\pi}\int_{-\pi}^{\pi}\log
f_{\theta}(\lambda)d\lambda-\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\log
f_{\theta}(\lambda_{j,n})\big{|}$ $\displaystyle\ \
+\sup_{\theta\in\Theta}\big{|}\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{f(\lambda)}{f_{\theta}(\lambda)}d\lambda-\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{f(\lambda_{j,n})}{f_{\theta}(\lambda_{j,n})}\big{|}=O(n^{-1}),$

where the last equality follows by the differentiability of
$1/f_{\theta}(\lambda)$ with respect to $\theta$, the boundedness properties
of $f_{\theta}\in{\mathcal{F}}_{\theta}$ and the fact that
$\sup_{\theta\in\Theta}\big{|}\partial/\partial\theta
f^{-1}_{\theta}(\lambda)\big{|}$ is bounded uniformly in $\lambda$. Using
$1/f_{\theta}(\lambda)\leq 1/\delta$ and Assumption 3, we get

$\displaystyle D_{2,n}$
$\displaystyle\leq\sup_{\lambda\in[-\pi,\pi]}\big{|}\widehat{f}(\lambda)-f(\lambda)\big{|}\sup_{\theta\in\Theta}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{1}{f_{\theta}(\lambda_{j,n})}$
$\displaystyle=O(1)\sup_{\lambda\in[-\pi,\pi]}\big{|}\widehat{f}(\lambda)-f(\lambda)\big{|}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0.$

To establish $D_{3,n}\stackrel{{\scriptstyle P}}{{\rightarrow}}0$, it suffices
to show that
$\sup_{\theta\in\Theta}|D_{n}(\theta,\widehat{f})-D(\theta,\widehat{f})|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$. For this we have

$\displaystyle\sup_{\theta\in\Theta}|D_{n}(\theta,\widehat{f})-D(\theta,\widehat{f})|$
$\displaystyle\leq\sup_{\theta\in\Theta}\big{|}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\log
f_{\theta}(\lambda_{j,n})-\frac{1}{2\pi}\int_{-\pi}^{\pi}\log
f_{\theta}(\lambda)d\lambda\big{|}$ $\displaystyle\ \
+\sup_{\theta\in\Theta}\big{|}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{\widehat{f}(\lambda_{j,n})}{f_{\theta}(\lambda_{j,n})}-\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{\widehat{f}(\lambda)}{f_{\theta}(\lambda)}d\lambda\big{|}$
$\displaystyle\leq
O(n^{-1})+\sup_{\theta\in\Theta}\big{|}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{1}{f_{\theta}(\lambda)}(\widehat{f}(\lambda_{j,n})-f(\lambda_{j,n}))\big{|}$
$\displaystyle\ \
+\sup_{\theta\in\Theta}\big{|}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{f(\lambda_{j,n})}{f_{\theta}(\lambda_{j,n})}-\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{f(\lambda)}{f_{\theta}(\lambda)}d\lambda\big{|}$
$\displaystyle\ \
+\sup_{\theta\in\Theta}\big{|}\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{1}{f_{\theta}(\lambda)}(f(\lambda)-\widehat{f}(\lambda))d\lambda\big{|}$
$\displaystyle=O(n^{-1})+O(1)\sup_{\lambda\in[-\pi,\pi]}\big{|}\widehat{f}(\lambda)-f(\lambda)\big{|},$

where the last equality follows because the second and the last term of the
last bound above is
$O(1)\sup_{\lambda\in[-\pi,\pi]}\big{|}\widehat{f}(\lambda)-f(\lambda)\big{|}$
and the third term is $O(n^{-1})$. Finally, $D_{4,n}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$ follows from
$\sup_{\theta\in\Theta}|D(\theta,\widehat{f})-D(\theta,f)|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$, which holds true since

$\displaystyle\sup_{\theta\in\Theta}|D(\theta,\widehat{f})-D(\theta,f)|$
$\displaystyle\leq\sup_{\theta\in\Theta}\big{|}\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{1}{f_{\theta}(\lambda)}(f(\lambda)-\widehat{f}(\lambda))d\lambda\big{|}$
$\displaystyle\leq
O(1)\sup_{\lambda\in[-\pi,\pi]}\big{|}\widehat{f}(\lambda)-f(\lambda)\big{|}.$

Consider (ii). Recall that by Assumption 3, the matrix $W$ is nonsingular. By
Lemma 4.2 of Lahiri (2003), $D^{(2)}_{n}(\widehat{\theta}_{0},I_{n}^{\ast})$
is nonsingular, if for $\delta>0$,

$\|D^{(2)}_{n}(\widehat{\theta}_{0},I_{n}^{\ast})-W\|_{F}\leq\delta/\|W^{-1}\|_{F},$

on a set with probability arbitrarily close to one for $n$ large enough. This
holds true since
$\|D^{(2)}_{n}(\widehat{\theta}_{0},I_{n}^{\ast})-W\|_{F}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$, in probability, see Lemma 7.2(i). Furthermore, by the
same Lemma 4.2, we have, on the same set, that, with probability arbitrarily
close to one,

(7.2)
$\|(D^{(2)}_{n}(\widehat{\theta}_{0},I_{n}^{\ast}))^{-1}\|_{F}\leq\|W^{-1}\|_{F}/(1-\delta)=2\|W^{-1}\|_{F},$

for $\delta=1/2$. Recall equation (7.1) and define on the set on which (7.2)
holds true the function

(7.3)
$w(\widehat{\theta}_{0}-\theta)=\big{(}D_{n}^{(2)}(\widehat{\theta}_{0},I_{n}^{\ast})\big{)}^{-1}\big{[}D_{n}^{(1)}(\widehat{\theta}_{0},I^{\ast}_{n})+\widetilde{R}_{n}(\theta)\big{]},\
\ \ \theta\in{\mathcal{B}}(\widehat{\theta}_{0},\delta),$

where
${\mathcal{B}}(\widehat{\theta}_{0},\delta)=\\{x\in{\mathbb{R}}^{m}|\|x\|\leq\delta\\}$
and

$\widetilde{R}_{n}(\theta)=\big{(}D_{n}^{(2)}(\theta^{+}_{n},I_{n}^{\ast})-D_{n}^{(2)}(\widehat{\theta}_{0},I_{n}^{\ast})\big{)}(\theta-\widehat{\theta}_{0}),$

with
$\|\widehat{\theta}_{0}-\theta^{+}_{n}\|\leq\|\widehat{\theta}_{0}-\theta\|$.
Using the expansion

$D_{n}^{(1)}(\theta,I^{\ast}_{n})=D_{n}^{(1)}(\widehat{\theta}_{0},I^{\ast}_{n})+D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})(\theta-\widehat{\theta}_{0})+\widetilde{R}_{n}(\theta),$

we can also express $w(\cdot)$ as

$w(\widehat{\theta}_{0}-\theta)=\big{(}D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})\big{)}^{-1}\big{[}D_{n}^{(1)}(\theta,I^{\ast}_{n})-D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})(\theta-\widehat{\theta}_{0})\big{]}.$

From this and because, by assumption, $\widehat{\theta}_{n}^{\ast}$ is the
unique solution of $D_{n}^{(1)}(\theta,I^{\ast}_{n})=0$, we get that for
$\theta=\widehat{\theta}_{n}^{\ast}$ it holds true that
$w(\widehat{\theta}_{0}-\widehat{\theta}_{n}^{\ast})=\widehat{\theta}_{0}-\widehat{\theta}_{n}^{\ast}$
and this is the unique solution of
$w(\widehat{\theta}_{0}-\theta)=\widehat{\theta}_{0}-\theta$ .

We next show that for $n$ large enough and with probability arbitrarily close
to one, a constant $C>0$ exists such that
$\|w(\widehat{\theta}_{0}-\theta)\|\leq C\log(n)/\sqrt{n}$ if
$\|\widehat{\theta}_{0}-\theta\|\leq C\log(n)/\sqrt{n}$. Toward this goal, we
get using (7.3) and the fact that
$D_{n}^{(1)}(\widehat{\theta}_{0},\widehat{f})=0$, the bound

(7.4) $\displaystyle\|w(\widehat{\theta}_{0}-\theta)\|\leq$
$\displaystyle\|\big{(}D_{n}^{(2)}(\widehat{\theta}_{0},I_{n}^{\ast})\big{)}^{-1}\|_{F}\Big{(}\|D_{n}^{(1)}(\widehat{\theta}_{0},I^{\ast}_{n})-D_{n}^{(1)}(\widehat{\theta}_{0},\widehat{f})\|+\|\widetilde{R}_{n}(\theta)\big{)}\|\Big{)}.$

Recall the bound
$\|(D^{(2)}_{n}(\widehat{\theta}_{0},I_{n}^{\ast}))^{-1}\|_{F}\leq
2\|W^{-1}\|_{F}$. Furthermore, for the first term in parentheses on the right
hand side of (7.4), we have

$\|D_{n}^{(1)}(\widehat{\theta}_{0},I^{\ast}_{n})-D_{n}^{(1)}(\widehat{\theta}_{0},\widehat{f})\|=\|Y^{\ast}_{n}\|,$

where

$Y^{\ast}_{n}=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{\partial}{\partial\theta}\frac{1}{f_{\theta}(\lambda_{j,n})}\Big{|}_{\theta=\widehat{\theta}_{0}}\widehat{f}(\lambda_{j,n})(U_{j}^{\ast}-1).$

It yields,

$\displaystyle n{\rm E}^{\ast}\|Y^{\ast}_{n}\|^{2}=$
$\displaystyle\frac{2}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{(}\frac{\partial}{\partial\theta}\frac{1}{f_{\theta}(\lambda_{j,n})}\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)}^{\top}\Big{(}\frac{\partial}{\partial\theta}\frac{1}{f_{\theta}(\lambda_{j,n})}\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)}\widehat{f}(\lambda_{j,n})^{2}$
$\displaystyle\stackrel{{\scriptstyle
P}}{{\rightarrow}}\frac{1}{\pi}\int_{-\pi}^{\pi}\Big{(}\frac{\partial}{\partial\theta}\frac{1}{f_{\theta}(\lambda)}\Big{|}_{\theta=\theta_{0}}\Big{)}^{\top}\Big{(}\frac{\partial}{\partial\theta}\frac{1}{f_{\theta}(\lambda)}\Big{|}_{\theta=\theta_{0}}\Big{)}f(\lambda)^{2},$

by Assumption 3, the continuity of the derivative and Lemma 7.1(i). This
implies by Markov’s inequality, that,

$\displaystyle{\rm
P}^{\ast}\big{(}\|D_{n}^{(1)}(\widehat{\theta}_{0},I^{\ast}_{n})-D_{n}^{(1)}(\widehat{\theta}_{0},\widehat{f})\|\geq
C\log(n)/\sqrt{n}\big{)}$ $\displaystyle\leq\frac{n{\rm
E}\|Y_{n}^{\ast}\|^{2}}{C^{2}\log^{2}(n)}=O_{P}(1/\log^{2}(n)).$

For the term $\widetilde{R}_{n}(\theta)$ in (7.4) we use the bound

$\|\widetilde{R}_{n}(\theta)\|\leq\big{(}\|M_{1,n}\|_{F}+\|M_{2,n}\|_{F}\big{)}\|\theta-\widehat{\theta}_{0}\|,$

where

$M_{1,n}=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\log
f_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\theta_{n}^{+}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\log
f_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)}$

and

$M_{2,n}=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\frac{1}{f_{\theta}(\lambda_{j,n})}\Big{|}_{\theta=\theta_{n}^{+}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\frac{1}{f_{\theta}(\lambda_{j,n})}\Big{|}_{\theta=\widehat{\theta}_{0}}\Big{)}I^{\ast}_{n}(\lambda_{j,n}).$

By the Lipschitz continuity of the second order derivatives, following from
Assumption 3, and since
$\|\theta_{n}^{+}-\widehat{\theta}_{0}\|\leq\|\theta-\widehat{\theta}_{0}\|$,
we get $\|\widetilde{R}_{n}(\theta)\|\leq C\|\widehat{\theta}_{0}-\theta\|$.
Hence if $\|\widehat{\theta}_{0}-\theta\|\leq C\log(n)/\sqrt{n}$ then, for $n$
large enough and with probability arbitrarily close to one, we have that,
$\|w(\widehat{\theta}_{0}-\theta)\|\leq C\log(n)/\sqrt{n}$. Consider next the
function $g:{\mathcal{B}}(0,1)\rightarrow{\mathcal{B}}(0,1)$ defined as

$g(x)=\frac{\sqrt{n}}{C\log(n)}w\big{(}\frac{C\log(n)}{\sqrt{n}}\cdot
x\big{)}.$

Notice that $g$ is continuous and that because for $x\in{\mathcal{B}}(0,1)$,
$\|C\log(n)x/\sqrt{n}\|\leq C\log(n)/\sqrt{n}$, we have,

$\|g(x)\|=\frac{\sqrt{n}}{C\log(n)}\|w(\frac{C\log(n)}{\sqrt{n}}x)\|\leq\frac{\sqrt{n}}{C\log(n)}\cdot\frac{C\log(n)}{\sqrt{n}}=1.$

By Bronwer’s fixed point Theorem, see Lahiri (2003), Proposition 4.1, there
exists $x_{0}\in{\mathcal{B}}(0,1)$ such that $g(x_{0})=x_{0}$, that is,

$w\big{(}\frac{C\log(n)}{\sqrt{n}}\cdot
x_{0}\big{)}=\frac{{C\log(n)}}{\sqrt{n}}x_{0}.$

Since $\widehat{\theta}_{0}-\widehat{\theta}^{\ast}_{n}$ is the unique
solution of $w(\widehat{\theta}_{0}-\theta)=\widehat{\theta}_{0}-\theta$, we
have that
$\widehat{\theta}_{0}-\widehat{\theta}^{\ast}_{n}=x_{0}C\log(n)/\sqrt{n}$,
that is,

$\|\widehat{\theta}_{0}-\widehat{\theta}^{\ast}_{n}\|\leq\frac{C\log(n)}{\sqrt{n}}\|x_{0}\|\leq\frac{C\log(n)}{\sqrt{n}}.$

Hence for $n$ large enough and with probability arbitrarily close to one, we
have

$\|\widehat{\theta}_{0}-\widehat{\theta}^{\ast}_{n}\|=O_{P}\Big{(}\frac{\log(n)}{\sqrt{n}}\Big{)},$

which converges to zero, as $n\rightarrow\infty$. ∎

######  Lemma 7.2.

Suppose that Assumption 1 to Assumption 4 are satisfied. Then, as
$n\rightarrow\infty$,

  1. (i)

$W^{\ast}_{n}\stackrel{{\scriptstyle P}}{{\rightarrow}}W$,

  2. (ii)

$V_{1,n}^{\ast}\stackrel{{\scriptstyle P}}{{\rightarrow}}V_{1}$,

  3. (iii)

$V^{+}_{2,n}\stackrel{{\scriptstyle P}}{{\rightarrow}}V_{2}$.

###### Proof.

Notice first that

(7.5)
$\sup_{\lambda\in[\pi,\pi]}\big{\|}g_{\widehat{\theta}_{0}}(\lambda)-g_{\theta_{0}}(\lambda)\big{\|}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0.$

and

(7.6)
$\sup_{\lambda\in[\pi,\pi]}\Big{\|}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\widehat{\theta}_{0}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\theta_{0}}\Big{\|}_{F}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0.$

To see why the above assertions hold true, observe first that by Lemma 7.1(i)
we have for $n$ large enough, $P(\widehat{\theta}_{0}\in
B(\theta_{0},\epsilon))\geq 1-\epsilon$, where
$B(\theta_{0},\epsilon)=\\{\theta:\|\theta-\theta_{0}\|\leq\epsilon\\}\subset\Theta$.
(7.5) and (7.6) follow then because the functions $g_{\theta}(\lambda)$ and
$\partial^{2}\big{/}(\partial\theta\partial\theta^{\top})f^{-1}_{\theta}(\lambda)$
are uniformly continuous on the compact set $[-\pi,\pi]\times
B(\theta_{0},\epsilon)$.

Consider (i). We have

$\displaystyle W^{\ast}_{n}$
$\displaystyle=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\Big{(}\log
f_{\theta}(\lambda_{j,n})+I^{\ast}_{n}(\lambda_{j,n})f^{-1}_{\theta}(\lambda)\Big{)}\Big{|}_{\theta=\theta_{0}}$
$\displaystyle\ \
+\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\log
f_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\widehat{\theta}_{0}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\log
f_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\theta_{0}}\Big{)}$ $\displaystyle\ \
+\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\widehat{\theta}_{0}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\theta_{0}}\Big{)}I_{n}^{\ast}(\lambda_{j,n})$
$\displaystyle=\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}\log
f_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\theta_{0}}+\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\theta_{0}}I^{\ast}_{n}(\lambda_{j,n})+o_{P}(1)$
$\displaystyle\stackrel{{\scriptstyle
P}}{{\rightarrow}}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}D(\theta_{0},f).$

The last equality follows because by (7.6),

$\displaystyle\Big{\|}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}$
$\displaystyle\Big{(}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\widehat{\theta}_{0}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda_{j,n})\Big{|}_{\theta=\theta_{0}}\Big{)}I_{n}^{\ast}(\lambda_{j,n})\Big{\|}_{F}$
$\displaystyle\leq\sup_{\lambda\in[-\pi,\pi]}\Big{\|}\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda)\big{|}_{\theta=\widehat{\theta}_{0}}-\frac{\partial^{2}}{\partial\theta\partial\theta^{\top}}f^{-1}_{\theta}(\lambda)\big{|}_{\theta=\theta_{0}}\Big{\|}_{F}\times\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}I^{\ast}_{n}(\lambda_{j,n})$
$\displaystyle=o_{P}(1)\cdot O_{P}(1).$

Consider (ii). Recall that $V_{1,n}^{\ast}={\rm Var}^{\ast}(M^{\ast}_{n})$ and
that ${\rm Cov}(I^{\ast}_{n}(\lambda_{j,n}),I^{\ast}_{n}(\lambda_{k,n}))={\bf
1}_{\\{|j|=|k|\\}}\widehat{f}(\lambda_{j,n})^{2}$. Hence

$\displaystyle{\rm Var}^{\ast}(M_{n}^{\ast})$
$\displaystyle=\frac{4\pi^{2}}{n}\sum_{j=1}^{N}\big{(}g_{\widehat{\theta}_{0}}(-\lambda_{j,n})g^{\top}_{\widehat{\theta}_{0}}(-\lambda_{j,n})+g_{\widehat{\theta}_{0}}(\lambda_{j,n})g^{\top}_{\widehat{\theta}_{0}}(-\lambda_{j,n})$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \
+g_{\widehat{\theta}_{0}}(-\lambda_{j,n})g^{\top}_{\widehat{\theta}_{0}}(\lambda_{j,n})+g_{\widehat{\theta}_{0}}(\lambda_{j,n})g^{\top}_{\widehat{\theta}_{0}}(\lambda_{j,n})\big{)}\widehat{f}(\lambda_{j,n})^{2}$
$\displaystyle=\frac{8\pi^{2}}{n}\sum_{j\in{\mathcal{G}}(n)}g_{\theta_{0}}(\lambda_{j,n})g^{\top}_{\theta_{0}}(\lambda_{j,n})f(\lambda_{j,n})^{2}+o_{P}(1)$
$\displaystyle\stackrel{{\scriptstyle P}}{{\rightarrow}}V_{1},$

where the last equality follows using the symmetry of $g_{\theta}(\lambda)$
with respect to $\lambda$, assertion (7.5) and Assumption 3.

Consider (iii). We have

$\displaystyle\Sigma_{n}^{+}$ $\displaystyle=\frac{b}{k}\sum_{\ell=1}^{k}{\rm
Var}^{*}\Big{(}\frac{2\pi}{b}\sum_{j\in{\mathcal{G}}(b)}g_{\widehat{\theta}_{0}}(\lambda_{j,b})I_{b}^{(\ell)}(\lambda_{j,b})\Big{)}$
$\displaystyle=\frac{4\pi^{2}}{b}\sum_{j_{1}\in{\mathcal{G}}(b)}\sum_{j_{2}\in{\mathcal{G}}(b)}g_{\widehat{\theta}_{0}}(\lambda_{j_{1},b})g^{\top}_{\widehat{\theta}_{0}}(\lambda_{j_{2},b}){\rm
Cov}^{\ast}\big{(}I_{b}^{(i_{1})}(\lambda_{j_{1},b}),I_{b}^{(i_{1})}(\lambda_{j_{2},b})\big{)}$
$\displaystyle=\frac{4\pi^{2}}{b}\sum_{j_{1}\in{\mathcal{G}}(b)}\sum_{j_{2}\in{\mathcal{G}}(b)}g_{\widehat{\theta}_{0}}(\lambda_{j_{1},b})g^{\top}_{\widehat{\theta}_{0}}(\lambda_{j_{2},b})$
$\displaystyle\ \ \ \ \ \ \
\times\frac{1}{n-b+1}\sum_{t=1}^{n-b+1}\big{\\{}I_{b}^{(t)}(\lambda_{j_{1},b}),I_{b}^{(t)}(\lambda_{j_{2},b}-\widetilde{f}(\lambda_{j_{1},b})\widetilde{f}(\lambda_{j_{1},b})\big{\\}}$
$\displaystyle=\frac{4\pi^{2}}{b}\sum_{j_{1}\in{\mathcal{G}}(b)}\sum_{j_{2}\in{\mathcal{G}}(b)}g_{\widehat{\theta}_{0}}(\lambda_{j_{1},b})g^{\top}_{\widehat{\theta}_{0}}(\lambda_{j_{2},b}){\rm
Cov}\big{(}I_{b}^{(1)}(\lambda_{j_{1},b}),I_{b}^{(1)}(\lambda_{j_{2},b})\big{)}+o_{P}(1)$

where the $o_{P}(1)$ term follows using Lemma 4.1 of Meyer et al. (2020) and
Assumption 4. Thus using the covariance properties of the periodogram
$I_{b}^{(1)}(\lambda_{j,b})$ of the subsample $X_{1},X_{2},\ldots,X_{b}$ and
(7.5), we get that, as $b\rightarrow\infty$,

$\frac{4\pi^{2}}{b}\sum_{j_{1}\in{\mathcal{G}}(b)}\sum_{j_{2}\in{\mathcal{G}}(b)}g_{\widehat{\theta}_{0}}(\lambda_{j_{1},b})g^{\top}_{\widehat{\theta}_{0}}(\lambda_{j_{2},b}){\rm
Cov}\big{(}I_{b}^{(1)}(\lambda_{j_{1},b}),I_{b}^{(1)}(\lambda_{j_{2},b})\big{)}\stackrel{{\scriptstyle
P}}{{\rightarrow}}V_{1}+V_{2}.$

By the same lemma we also have

$(n-b+1)^{-1}\sum_{t=1}^{n-b+1}I_{b}^{(t)}(\lambda_{j,b})^{2}/\widetilde{f}_{b}(\lambda_{j,b})^{2}\stackrel{{\scriptstyle
P}}{{\rightarrow}}2,$

which implies using (7.5) again and Assumption 3, that
$C_{n}^{+}\stackrel{{\scriptstyle P}}{{\rightarrow}}V_{1}$, in probability.
The assertion follows then since $V^{+}_{2,n}=\Sigma^{+}_{n}-C_{n}^{+}$. ∎

Proof of Theorem 4.1: In view of the definition of $L_{n}^{\ast}$ and Lemma
7.2, to establish the assertion of the theorem it suffices to show that,

(7.7)
$\sqrt{n}(\widehat{\theta}_{n}^{\ast}-\widehat{\theta}_{0})\stackrel{{\scriptstyle
D}}{{\rightarrow}}{\mathcal{N}}(0,W^{-1}V_{1}W^{-1}).$

Since $D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})=W^{\ast}_{n}$ we get by
Lemma 7.2(i), that
$\|D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})-W\|_{F}\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$ and that, for $n$ large enough,
$D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})$ is nonsingular; see also the
proof of Lemma 7.1(ii). Recall that
$D^{(1)}_{n}(\widehat{\theta}_{n}^{\ast},I^{\ast}_{n})=0$ and that
$D^{(1)}_{n}(\widehat{\theta}_{0},\widehat{f})=0$. We get using the expansion
(7.1), that

$\displaystyle\big{[}I_{m}+\big{(}D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})\big{)}^{-1}$
$\displaystyle
R_{n}(\widehat{\theta}^{\ast}_{n})\big{]}\sqrt{n}(\widehat{\theta}_{n}^{\ast}-\widehat{\theta}_{0})$
(7.8)
$\displaystyle=-\big{(}D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})\big{)}^{-1}\sqrt{n}\big{(}D^{(1)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})-D^{(1)}_{n}(\widehat{\theta}_{0},\widehat{f})\big{)}.$

Since
$\big{(}D_{n}^{(2)}(\widehat{\theta}_{0},I^{\ast}_{n})\big{)}^{-1}\stackrel{{\scriptstyle
P}}{{\rightarrow}}W^{-1}$, in probability, assertion (7.7) follows from
expression (7), if we show that

(7.9)
$-\sqrt{n}\big{(}D^{(1)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})-D^{(1)}_{n}(\widehat{\theta}_{0},\widehat{f})\big{)}\stackrel{{\scriptstyle
D}}{{\rightarrow}}{\mathcal{N}}(0,V_{1}),$

and

(7.10) $R_{n}(\widehat{\theta}_{n}^{\ast})=o_{P}(1),$

in probability.

For (7.9) we have,

$\displaystyle-\sqrt{n}\big{(}D^{(1)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})-D^{(1)}_{n}(\widehat{\theta}_{0},\widehat{f})\big{)}$
$\displaystyle=-\frac{1}{\sqrt{n}}\sum_{j\in{\mathcal{G}}(n)}\frac{\partial}{\partial\theta}\frac{1}{f_{\theta}(\lambda_{j,n})}\Big{|}_{\theta=\widehat{\theta}_{0}}(I^{\ast}_{n}(\lambda_{j,n})-\widehat{f}(\lambda_{j,n}))$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=1}^{N}W_{j,n}(U^{\ast}_{j}-1),$

where
$W_{j,n}=(g_{\widehat{\theta}_{0}}(\lambda_{j,n})+g_{\widehat{\theta}_{0}}(-\lambda_{j,n}))\widehat{f}(\lambda_{j,n})$
and the $U^{\ast}_{j}$’s are i.i.d. (7.9) follows then because ${\rm
E}^{\ast}(W_{j,n}(U^{\ast}_{j}-1))=0$, ${\rm
Var}^{\ast}(W_{j,n}(U^{\ast}_{j}-1))\stackrel{{\scriptstyle
P}}{{\rightarrow}}V_{1}$ by the same arguments as those used in the proof that
${\rm Var}^{\ast}(M_{n}^{\ast})\stackrel{{\scriptstyle
P}}{{\rightarrow}}V_{1}$ in Lemma 7.2(ii), and because

$\displaystyle\sum_{j=1}^{N}{\rm
E}^{\ast}\|\frac{1}{\sqrt{n}}W_{j,n}(U^{\ast}_{j}-1)\|^{3}$
$\displaystyle={\rm
E}^{\ast}(|U^{\ast}_{j}-1|^{3})\frac{1}{n^{3/2}}\sum_{j=1}^{N}\|W_{j,n}\|^{3}=O_{P}(n^{-1/2}),$

verifies Liapunov’s condition.

Consider (7.10) and observe that
$R_{n}(\widehat{\theta}^{\ast}_{n})=D^{(2)}_{n}(\theta^{+}_{n},I^{\ast}_{n})-D^{(2)}_{n}(\widehat{\theta}_{0},I^{\ast}_{n})$
for some $\theta^{+}_{n}\in\Theta$ such that
$\|\theta^{+}_{n}-\widehat{\theta}_{0}\|\leq\|\widehat{\theta}^{\ast}_{n}-\widehat{\theta}_{0}\|$.
Denote by $d_{j,k}(\mathring{\theta})$ the $(j,k)$th element of the matrix
$D^{(2)}_{n}(\mathring{\theta},I^{\ast}_{n})$, that is,
$d_{j,k}(\mathring{\theta})=\partial^{2}\big{/}(\partial\theta_{j}\partial\theta_{k})D_{n}(\theta,I^{\ast})|_{\theta=\mathring{\theta}}$.
We then have for the $(j,k)$th element $r_{j,k}(\widehat{\theta}_{n}^{\ast})$
of the matrix $R_{n}(\widehat{\theta}^{\ast}_{n})$,

$\displaystyle|r_{j,k}(\widehat{\theta}_{n}^{\ast})|$
$\displaystyle=|d_{j,k}(\theta^{+}_{n})-d_{j,k}(\widehat{\theta}_{0})|$
$\displaystyle\leq\sup_{\lambda\in[-\pi,\pi]}\big{|}\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}\log
f_{\theta}(\lambda)\Big{|}_{\theta=\theta^{+}_{n}}-\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}\log
f_{\theta}(\lambda)\Big{|}_{\theta=\widehat{\theta}_{0}}\big{|}$
$\displaystyle\ \
+\sup_{\lambda\in[-\pi,\pi]}\big{|}\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\theta^{+}_{0}}-\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{k}}f^{-1}_{\theta}(\lambda)\Big{|}_{\theta=\widehat{\theta}_{0}}\big{|}\frac{1}{n}\sum_{j\in{\mathcal{G}}(n)}I_{n}^{\ast}(\lambda_{j,n})$
$\displaystyle=o_{P}(1),$

as $n\rightarrow\infty$, since
$n^{-1}\sum_{j\in{\mathcal{G}}(n)}I_{n}^{\ast}(\lambda_{j,n})=O_{P}(1)$, the
second order partial derivative functions
$\partial^{2}\big{/}(\partial\theta_{j}\partial\theta_{k})\log
f_{\theta}(\lambda)$ and
$\partial^{2}\big{/}(\partial\theta_{j}\partial\theta_{k})f^{-1}_{\theta}(\lambda)$
are uniformly continuous of the compact set $[-\pi,\pi]\times\Theta$ and
$\|\theta^{+}_{n}-\widehat{\theta}_{0}\|\leq\|\widehat{\theta}^{\ast}_{n}-\widehat{\theta}_{0}\|\stackrel{{\scriptstyle
P}}{{\rightarrow}}0$, in probability, by Lemma 7.1(ii).  o

## References

  * [1] AuestadTjostheim1990AuesTjos190 Auestad, B. and Tjostheim, D. (1990). Identification of Nonlinear Time Series: First Order Characterization and Order Determination. Biometrika 77, 669–687.
  * [2] Beltrão and Bloomfield1987BeltBloo1987 Beltrão, K. I. and Bloomfield, P. (1987). Determining the Bandwidth of a Kernel Spectrum Estimate. Journal of Time Series Analysis 8, 21–38.
  * [3] Brockwell and Davis1991BrockwellDavis1991 Brockwell, P. and Davis, R.A. (1991). Time Series: Theory and Methods (2nd edition). Springer, New York.
  * [4] Dahlhaus1985Dahlhaus1985 Dahlhaus, R. (1985). Asymptotic Normality of Spectral Estimates. Journal of Multivariate Analysis 16, 412–431.
  * [5] Dahlhaus1988Dahl1988 Dahlhaus, R. (1988). Small Sample Effects in Time Series Analysis: A New Asymptotic Theory and a New Estimate. Annals of Statistics 16, 808–841.
  * [6] Dahlhaus and Janas1996DahlJanas1996 Dahlhaus, R. and Janas, D. (1996). A Frequency Domain Bootstrap for Ratio Statistics in Time Series Analysis. The Annals of Statistics 24, 1934–1963.
  * [7] DahlhausWeelmeyer1996DahlWefe1996 Dahlhaus, R. and Wefelmeyer, W. (1996). Asymptotically Optimal Estimation in Misspecified Time Series Models. Annals of Statistics 24, 952–974.
  * [8] Doukhan and Leon1989DoukLeon1989 Doukhan, P. and León, J. R. (1989). Cumulants for Stationary Mixing Random Sequences and Applications to Empirical Spectral Density. Probability and Mathematical Statistics 10, 11–26.
  * [9] FanYao2005FanYao2005 Fan, Y. and Yao, Q. (2005). Nonlinear Time Series. Springer-Verlag, New York.
  * [10] Franke and Härdle1992FrankeHaerdle1992 Franke, J. and Härdle, W. (1992). On Bootstrapping Kernel Spectral Estimates. The Annals of Statistics 20, 121–145.
  * [11] Hannan1973Hann1973 Hannan, E. (1973). The Asymptotic Theory of Linear Time Series Models. Journal of Applied Probability, 10, 130–145.
  * [12] Hosoya1979Hosoya1979 Hosoya, Y. (1979). High-order Efficiency in the Estimation of Linear Processes. Annals of Statistics 7, 516–530.
  * [13] Hurvich and Zeger1987HurvichZeger1987 Hurvich, C.M. and Zeger, S.L. (1987). Frequency Domain Bootstrap Methods for Time Series. Preprint, Department of Statistics and Operations Research, New York University.
  * [14] KimNordman2013KimNord2013 Kim, Y. M. and Nordman, D. (2013). A Frequency Domain Bootstrap for Whittle Estimation Under Long-Range Dependence. Journal of Multivariate Analysis 115, 405–420.
  * [15] Lahiri2003Lahi 2003 Lahiri, S. N. (2003). Resampling Methods for Dependent Data. Springer-Verlag, New York.
  * [16] MeyerEtAl2020MeyerEtAl2020 Meyer, M., E. Paparoditis and Kreiss, J-P. (2020). Extending the Validity of Frequency Domain Bootstrap Methods to General Stationary Processes. Annals of Statistics, 48, 2404–2427.
  * [17] Newton and Pagano1983NewPag1983 Newton, H. J. and Pagano, M. (1983). A Method for Determing Periods in Time Series. Journal of the American Statistical Association 78, 152–157.
  * [18] Paparoditis and Politis1999PapPol1999 Paparoditis, E. and Politis, D.N. (1999). The Local Bootstrap for Periodogram Statistics. Journal of Time Series Analysis 20, 193–222.
  * [19] Rosenblatt1985Rose1985 Rosenblatt, M. (1985). Stationary Sequences and Random Fields. Birkhäuser, Boston.
  * [20] Shao, X.2010Shao2010 Shao, X. (2010). Nonstationarity-extended Whittle estimators. Econometric Theory, 26. 1060–1087.
  * [21] SubbaRaoYang2020SubbYang2020 Subba Rao, S. and Yang, J. (2020). Reconciling the Gaussian and Whittle Likelihood with an Application to Estimation in the Frequency Domain. Annals of Statistics, to appear.
  * [22] SykulskiEtAl2019Suku2019 Sykulski, A. M., S. C., Ohlede, A. P. Guillaumin, J. M. Lilly and Early, J. J. (2019). The Debiased Whittle Likelihood. Biometrika 106, 251–266.
  * [23] Taniguchi1987Taniguchi1987 Taniguchi, M. (1987). Minimum Contrast Estimation for Spectral Densities of Stationary Processes. Journal of the Royal Statistical Society, Series B 49, 315–325.
  * [24] Tewes, J., Nordmann D. J. and Politis, D. N.2019TeNoPo2019 Tewes, J., Nordman, D. J. and Politis, D. N. (2019). Convolved Subsampling Estimation with Applications to Block Bootstrap. Annals of Statistics, 47, 468-496.
  * [25] VelascoRobinson2000VelRobi2000 Velasco, C. and Robinson, P. (2000). Whittle Pseudo-Maximum Likelihood Estimation for Nonstationary Time Series. Journal of the American Statistical Association 95, 1229–1243.
  * [26] Walker1964Walk1964 Walker, A. M. (1964). Asymptotic Properties of Least-Squares Estimates of Parameters of the Spectrum of a Stationary Non-Deterministic Time-Series. Journal of the Australian Mathematical Society, 4, 363–384.
  * [27] Whittlel1951Whittle51 Whittle, P. (1951). Hypothesis Testing in Time Series Analysis. Thesis, Uppsala University.
  * [28] Whittlel1953Whittle53 Whittle, p. (1953). The Analysis of Multiple Stationary Time Series. Journal of the Royal Statistical Society, Series B, 15, 125–139.
  * [29] WuShao2004WuShao2004 Wu, B. W. and Shao, X. (2004). Limit Theorems and Iterated Random Functions. Journal of Applied Probability, 41, 425–436.
  * [30]