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@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.
