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

[Bert Gunter]

Have a look at the CRAN Time Series Task View:

https://cran.r-project.org/web/views/TimeSeries.html

Generally speaking, R-help is for help on R programming, not detailed
statistical questions, so it is less likely that your query would receive a
useful answer here.

Cheers,
Bert

On Wed, Jul 10, 2024 at 1:59 AM Ashim Kapoor <ashimkapoor using gmail.com> wrote:

> 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
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

	[[alternative HTML version deleted]]