Eric Thompson
2007-Nov-21 21:12 UTC
[R] Different freq returned by spec.ar() and spec.pgram()
Dear list, I've recently become interested in comparing the spectral estimates using the different methods ("pgram" and "ar") in the spectrum() function in the stats package. With many thanks to the authors of these complicated functions, I would like to point out what looks to me like a bit of an inconsistency -- but I would not be surprised if there is good reasoning that justifies it that I am just not seeing right now. If we use the lh data, the two methods return similar results:> spectrum(lh, col = "blue") > spec.ar(lh, add = TRUE)But using the ldeaths data:> spectrum(ldeaths, col = "blue") > spec.ar(ldeaths, add = TRUE)the resulting plots do not compare over the same frequency range. This results because spec.ar defines frequency as> freq <- seq.int(0, 0.5, length.out = n.freq)whereas spec.pgram uses> xfreq <- frequency(x) > N <- nrow(x) > Nspec <- floor(N/2) > freq <- seq.int(from = xfreq/N, by = xfreq/N, length.out = Nspec)And so the reason the spectral estimates of lh are similar is that frequency(lh) = 1, whereas frequency(ldeaths) = 12. The documentation seems more extensive for spec.pgram (and the pertinent section in MASS focuses on spec.pgram), and I realize that there is a warning in ?spec.ar that AR spectra can be misleading. But is there a reason that I am not aware of that the frequencies of the AR spectra are defined in this way? It seems to me that it would be desirable for frequency to be defined over the same range as in spec.pgram. All that would need to be added would be a line to scale the freq vector using the sampling frequency before it is returned. Eric Thompson Graduate Student Dept. of Civil & Env. Eng. Tufts University> sessionInfo()R version 2.6.0 (2007-10-03) i686-pc-linux-gnu locale: LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] datasets utils stats graphics grDevices methods base other attached packages: [1] MASS_7.2-37 loaded via a namespace (and not attached): [1] rcompgen_0.1-17
Prof Brian Ripley
2007-Nov-21 22:03 UTC
[R] Different freq returned by spec.ar() and spec.pgram()
Frequency has units (1/wavelength), and this is just a question of using different units. For monthly series, is this per month or per year? The functions were written for R at different times by different people: it would seem to make sense to alter spec.ar to use the per year interpretation. It would be unusual to use spec.ar on a series with known frequency (I have never seen it done), and it is arguable that the per-lag interpretation is more natural for that method. On Wed, 21 Nov 2007, Eric Thompson wrote:> Dear list, > > I've recently become interested in comparing the spectral estimates > using the different methods ("pgram" and "ar") in the spectrum() > function in the stats package. > > With many thanks to the authors of these complicated functions, I > would like to point out what looks to me like a bit of an > inconsistency -- but I would not be surprised if there is good > reasoning that justifies it that I am just not seeing right now. If we > use the lh data, the two methods return similar results: > >> spectrum(lh, col = "blue") >> spec.ar(lh, add = TRUE) > > But using the ldeaths data: > >> spectrum(ldeaths, col = "blue") >> spec.ar(ldeaths, add = TRUE) > > the resulting plots do not compare over the same frequency range. This > results because spec.ar defines frequency as > >> freq <- seq.int(0, 0.5, length.out = n.freq) > > whereas spec.pgram uses > >> xfreq <- frequency(x) >> N <- nrow(x) >> Nspec <- floor(N/2) >> freq <- seq.int(from = xfreq/N, by = xfreq/N, length.out = Nspec) > > And so the reason the spectral estimates of lh are similar is that > frequency(lh) = 1, whereas frequency(ldeaths) = 12. > > The documentation seems more extensive for spec.pgram (and the > pertinent section in MASS focuses on spec.pgram), and I realize that > there is a warning in ?spec.ar that AR spectra can be misleading. But > is there a reason that I am not aware of that the frequencies of the > AR spectra are defined in this way? It seems to me that it would be > desirable for frequency to be defined over the same range as in > spec.pgram. All that would need to be added would be a line to scale > the freq vector using the sampling frequency before it is returned. > > Eric Thompson > Graduate Student > Dept. of Civil & Env. Eng. > Tufts University > >> sessionInfo() > R version 2.6.0 (2007-10-03) > i686-pc-linux-gnu > > locale: > LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C > > attached base packages: > [1] datasets utils stats graphics grDevices methods base > > other attached packages: > [1] MASS_7.2-37 > > loaded via a namespace (and not attached): > [1] rcompgen_0.1-17 > > ______________________________________________ > R-help at r-project.org mailing list > 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. >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595