[R] discrepancy between periodogram implementations ? per and spec.pgram
Lieven Desmet
lieven.desmet at wis.kuleuven.be
Fri Dec 14 15:31:59 CET 2007
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
More information about the R-help
mailing list