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

[Lieven Desmet]

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