[R] Implementation for selecting lag of a lag window spectral estimator using generalized cross validation (using deviance)

Ashim Kapoor @@h|mk@poor @end|ng |rom gm@||@com
Wed Jul 10 10:58:28 CEST 2024


Dear All,

I am looking for:

A software to select the lag length for a lag window spectral estimator.
Also, I have a small query in the reprex given below.

Background for the above, from the book by Percival and Walden:

1. We are given X_1,...,X_n which is one realization of a stochastic process.
2. We may compute the periodogram using FFT, for example by the
function spectrum in R.
3. The above is badly biased so we taper X_1,...,X_n to reduce the
bias in the periodogram.
4. Now that the bias in under control, we focus on reducing the
variance. So we take a window like for eg. the Parzen window, and
choose
a lag length m which controls the amount of smoothing across frequencies.
5. One way of choosing m is mentioned in :
https://web.archive.org/web/20080221221221id_/http://www.stat.uiuc.edu/~ombao/PAPERS.dir/gcvbmka.pdf

I am looking for an R package which implements 5.

Here is a reprex:

# 1.  Periodogram which may be biased
plot(spectrum(lh,taper= 0, method="pgram"),log="dB")

# 2. Using the default in built cosine taper
plot(spectrum(lh,taper = .3, method="pgram"),log="dB")

# 2. Again, using slepian taper
library(multitaper)
# I choose: n = length(lh), k =1, nw=2
mytaper = dpss(n=length(lh), k=1 , nw=2, returnEigenvalues=TRUE)
# Tapered series
lh * mytaper$v
# I may compute the spectrum with reduced bias as:
plot(spectrum(lh*mytaper$v,method="pgram"),log="dB")

# We now focus on the variance
# For a fixed m = 10, using a Parzen window.
library(gsignal)
parzenwin(10)

# The following 2 lines of code, where I try to do the same thing in 2
different ways, did not work for me:

kernapply( spectrum(lh*mytaper$v,method="pgram")$spec,parzenwin(10),circular=TRUE)
spectrum(lh*mytaper$v,kernel  = parzenwin(10),method="pgram")

# ?spec.pgram says
kernel: alternatively, a kernel smoother of class ‘"tskernel"’.

How can I see all available kernels of class tskernel ?

The important question here is how to choose m which implies a bias -
variance tradeoff. Ombao et al, have a generalized cross validation
method to choose m.
Please see point 5 above. Does that exist in R ?

Many thanks,
Ashim



More information about the R-help mailing list