Re: [R] discrepancy between periodogram implementations ? per and spec.pgram

From: Lieven Desmet <lieven.desmet_at_wis.kuleuven.be>
Date: Fri, 14 Dec 2007 15:31:59 +0100

Prof Brian Ripley wrote:

> There are several definitions of a periodgram. Note that
>
>> log(2*pi)
>
> [1] 1.837877
>
> See the comments in ?spectrum about scalings.
>
> I think the comments in ?per incorrectly ignore the scaling issues:
> per() does not take the base frequency into account and has an extra
> divisor of 2*pi. E.g.
>
>> x <- rnorm(64)
>> spec.pgram(x, taper=0, detrend=F)$spec/per(x)[-1]
>
> [1] 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185
> 6.283185
> [9] 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185
> 6.283185
> [17] 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185
> 6.283185
> [25] 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185 6.283185
> 6.283185
>
>
> On Wed, 12 Dec 2007, Lieven Desmet wrote:
>
>> hello,
>>
>> I have been using the per function in package longmemo to obtain a
>> simple raw periodogram.
>> I am considering to switch to the function spec.pgram since I want to be
>> able to do tapering.
>> To compare both I used spec.pgram with the options as suggested in the
>> documentation of per {longmemo} to make them correspond.
>>
>> Now I have found on a variety of examples that there is a shift between
>> the log of the periodogram with per and that with spec.pgram. This
>> vertical shift amounts to approx. 1.8 on the log scale (the graph of
>> spec.pgram being above the one from per).
>>
>> Is there some explanation for this ? Is the one from spec.pgram the
>> better one as suggested in the documentation of per {longmemo}? Finally
>> how are these related to an estimate of the spectral density obtained
>> from spec.arima ?
>
>
> What is spec.arima? If you meant spec.ar, that is on the same scale
> as spec.pgram for series with base frequency 1 (and for all series for
> R >= 2.7.0).
>
>
>> Many thanks for help and clarification.
>>
>> Lieven Desmet
>
>
Dear Prof. Ripley,

thanks very much for a quick and helpful response. In the last question I wanted to hint at
specARIMA which I am using to get the theoretical spectral density of an ARMA process.

This works very well in general, however, in a simple example

X_t=0.7*X_{t-1}+epsilon_t

I obtain a value 1.768253 for funscaled[1] ( the first Fourier frequency 0.003141593)

using

str(f<-specARIMA(eta=c(H=0.5,phi=c(0.7),psi=c()),p=1,q=0,m=2000)) funscaled<-numeric(length(f$freq))
funscaled<-f$spec*f$theta1

where the theoretical value should be 0.901878 with

b<-0.7
omega<-0.003141593
1/(2*pi)*(1-b^2)/(1+b^2-2*b*cos(omega))
[1] 0.9018088

using the formula (2.40) in Fan and Yao, Nonlinear Time Series ( Springer 2003 ), page 54-55

Is there also a simple explanation for this ? am I overlooking something ?

Thanks and best regards,

Lieven Desmet,

maths dept - KULeuven - Belgium

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm



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. Received on Fri 14 Dec 2007 - 14:39:49 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 14 Dec 2007 - 15:30:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.