Confidence interval of Spectrum

[ivory4]

Is there a procedure to calculate the confidence interval for estimated spectrum?

[TomDoan]

The following estimates the spectrum with confidence bands for the series of sunspot numbers. This is based upon the standard smoothed periodogram procedure.

bd_tsap354.rpf

(1.07 KiB) Downloaded 1174 times

with data file

sunspots.dat

(393 Bytes) Downloaded 1150 times

If what you want is a confidence interval for the spectral density estimated indirectly through an ARMA model, that is discussed on page 334 of Hannan's Multiple Time Series with the comment: "This is a relatively complicated expression reflecting, it appears, the use of a method of estimation not basically designed for spectral estimation."

[ivory4]

Thanks. It seems that this one is not included in 7.3's textbook example folder.

why it seems much smoother in @armaspectrum calculated spectrum than @spectrum.

I am asking because the Manual says that @armaspectrum calculates transfer function which I thought is a function of Cosine? Like AR(1) case 1/(1+phi^2+2phi*cos(w))? That's why it is smoother?

[ivory4]

How the transfer function is caculated in Rats?

[TomDoan]

Thanks. It seems that this one is not included in 7.3's textbook example folder.

No. The more advanced of the two Brockwell & Davis books (from which this came) doesn't really have many worked examples, and the less advanced doesn't deal with the asymptotic distribution of the spectral density.

why it seems much smoother in @armaspectrum calculated spectrum than @spectrum.

I am asking because the Manual says that @armaspectrum calculates transfer function which I thought is a function of Cosine? Like AR(1) case 1/(1+phi^2+2phi*cos(w))? That's why it is smoother?

From a low order ARMA process, the entire spectrum is a function of just a very small number of parameters. By constrast, the periodogram with N data points has N (theoretically) independently estimated parameters. (That's exactly true if you use N ordinates; if you pad, it's almost true). Smoothing reduces the independence (trading some bias for variance reduction), but doesn't change the fact that you're still basing your calculation on quite a few parameters.

cos(w) is quite a smooth function over [0,2pi]. It's only when you get into cos(kw) for k large that you get a very erratic function. And those will be components of the power spectrum computed using Fourier transforms.

[ivory4]

I read in the manual that transfer function is calculated using FT of the lag polynomial [1-A(L)] and then divided by B(L), given yt=A(L)yt+B(L)ut.
In a low dimention case, like AR(2) or ARMA(2,2) I would like to calculate the spectrum in an analytical way using cosine function and ARMA coefficients because if armaspectrum is not working for a slightly modified ARMA model

[ivory4]

or is there a way to plot analytical spectrum using Phi(exp(-iw))^2 explicitly?

[TomDoan]

There's no analytical spectrum for an ARIMA model with differences. (The power at 0 frequency is "infinite" and it isn't integrable). If you try to apply ARMASpectrum to a non-stationary equation, it won't work because it will try to divide by 0 at 0 frequency. You need to difference first, then estimate the underlying ARMA model---it will have a well-defined spectrum.

[ivory4]

Code: Select all

    open data lgdp.txt
    calendar(q) 1947
    data(format=free,org=columns) 1947:01 1998:02 lgdp
    diff lgdp / dlgdp
    compute obs=206
    * ARMA(2,2) Record the coefficients
    boxjenk(ar=2,diffs=1,ma=2,maxl,constant,define=lgdpeq) dlgdp
    compute phi1=%beta(2),phi2=%beta(3),theta1=%beta(4),theta2=%beta(5),se=%seesq

    *Analytical form


    compute cst1 = 1+theta1^2-2*theta2+theta2^2
    compute cst2 = 1+phi1^2+2*phi2+phi2^2
    set time / = t

    set cosine / = cos((time-1)/obs*2*%pi)

    set nume / = cst1+2*(theta1+theta1*theta2)*cosine+4*theta2*cosine*cosine
    set deno / =  cst2+2*(phi1*phi2-phi1)*cosine           -4*phi2^2*cosine*cosine
    set spec / = se/(2*%pi)*nume/deno


    set frequencies / = (time-1)/104
    scatter(style=lines,header="Spectrum of Differenced Yt",hlabel="Fractions of pi") 1
       # frequencies spect 1 104

I changed it and to calculate the analytical one explicitly in this way, is it working? Especially the way to define cosine(w). I want to calculate the analytical spectrum in this way because for some modified model, @ARMAspectrum may not be available.

[TomDoan]

First off, I wouldn't recommend trying to expand out 2nd degree polynomials in z transforms---the chances for an error are rather high. You'll probably find it easier to use CSET with the %ZLAG function to create your own transfer functions, and then let CMULT multiply everything out in the frequency domain. The FRACTINT.PRG example uses %ZLAG rather extensively.

If you want to check your code, take a case where the armaspectrum procedure applies, and see if you get the same answer.

[ivory4]

If what you want is a confidence interval for the spectral density estimated indirectly through an ARMA model, that is discussed on page 334 of Hannan's Multiple Time Series with the comment: "This is a relatively complicated expression reflecting, it appears, the use of a method of estimation not basically designed for spectral estimation."

Analytical spectrum implied by a ARMA model is hard to obtain in an explicit way as you mentioned. How about I generate data (Boot Strapping errors) from the ARMA model and estimate the ARMA coefficients and Std of error iteratively and use the estimated parameters(ARMA coefficients and Std) to draw spectrum.

For each point of frequency, use a grid search to find out the XX% confidence interval?

[TomDoan]

That would work. An example of bootstrapping an ARMA model is given at:

http://www.estima.com/forum/viewtopic.php?f=5&t=752

Note that for what you're doing, you don't need to add back the mean in generating bootstrapped data.