rolf@natmlab.OZ (Rolf Turner) (01/26/86)
A rather strange phenomenon has developed around a home-grown
S function, named 'spaut' (spectrum -- autoregressive). The code for this
function is given below. It worked fine under an older version of S,
but when I recompiled it under a newer version (tape dated approx. Dec.
1984) it gave totally wrong answers. After fiddling with it for a while I
discovered (more or less by accident) that if the declaration for `at'
and `sincos' was changed from
STATIC( real at(2), sincos(2) )
to
STRUCTURE(at/REAL,2/,sincos/REAL,2/)
it worked again. (Either one works fine if compiled using the old S.)
Can anyone out there explain why this phenomenon happens, and how to fix
things so that it doesn't happen? That STATIC declaration ***should***
work! One wonders about the validity of results produced by other functions
containing STATIC declarations.
The newer S makes use of a version 4.2 Fortran 77 compiler, with SOME of
the bugs fixed, whereas the older one uses a version 4.1 compiler.
It has been suggested that the problem has something to do with the
'optimization' which the 4.2 compiler performs.
To test things, try doing
v1<-c(-2.7607,3.8106,-2.6535,0.9238)
v2<-spaut(arp=v1,var=1.)
plot(v2,type="l")
The correct graph is shown on page 609 of "Spectral Analysis and Time Series",
by M. B. Priestley, Academic Press, 1981.
##############################################################################
Code for spaut:
FUNCTION spaut(
arstr /STR,OPTIONAL/
arp /REAL,OPTIONAL/
var /REAL,OPTIONAL/
nfreq /INT,1,100/
)
if(!MISSING(arstr)) {
if(!MISSING(arp)|!MISSING(var))FATAL(Specify arstr OR arp and var.)
STRUCTURE(
arp/REAL,FROM(arstr)/
pvar/REAL,FROM(arstr)/
)
STRUCTURE(
var/REAL,1/
)
}
lp = LENGTH(arp)
if(MISSING(var)) var[1] = pvar[lp]
zero=0.0
pi=3.14159
STRUCTURE(spec/REAL,nfreq/,freq/REAL,nfreq/)
STATIC( real at(2), sincos(2) ) # This is the line that makes
# things go funny.
do j=1,nfreq{
omega= pi*float(j)/float(nfreq)
freq[j] = omega
do n=1,2{ at[n] = zero }
do k=1,lp{
ang=omega*float(k)
sincos[1]=cos(ang)
sincos[2]=sin(ang)
at[1] = sincos[1] * arp[k] + at[1]
at[2] = sincos[2] * arp[k] + at[2]
}
spec[j]=(at[1]+1)**2+at[2]**2
spec[j]=var/(2.*pi*spec[j])
}
RETURN(x=freq,y=spec)
END