SECTION 18 . TIME SERIES
The name of the function SUBROUTINE which computes the function value for conditional likelihood maximization for an arma(p,q) model is ARMAPQ; hence the optimization algorithms (such as GRADX, see Section 2.) should be called with the function subroutine name ARMAPQ. Since function subroutines to be used with optimization algorithms are called by GQOPT as CALL FUNC(X,NP,FU0,*) where X is the variable vector with respect to which optimization is carried out, NP is the number of variables, and FU0 is the value of the function, other parameters of interest for an arma(p,q) estimation problems must be communicated to the function subroutine via labelled COMMON (see below). The model, the parameters of which are to be estimated, is in general
YT = X(1)+X(2)*YT-1+...+X(1+IP)*YT-IP+UT+X(2+IP)*UT-1+ +...+X(1+IP+IQ)*UT-IQ
Hence, the paramaters to be estimated are the constant, IP autoregressive parameters and IQ moving average parameters. The various paramaters to be communicated to ARMAPQ have to be set in the MAIN program, which must include the following labelled COMMON blocks:
COMMON/BARMA1/IP,IQ,IBJ,IT,IROOTS
COMMON/BARMA2/Y(IT)
COMMON/BARMA3/YY(IT)
COMMON/BARMA4/EPSIL(IT)
where
IP = the number of autoregressive lags
IQ = the order of the moving average
IBJ = 0 if the user wishes to use the actual values the Y's for
the first IP vaues, and
IBJ = 1 if the first IP values are to be replaced by their
expected values (Box-Jenkins)
IT = the total number of observations
IROOTS = 0 if the program should check at every function
evaluation whether the roots of the moving average polynomial are
inside the unit circle,
IROOTS = 1 if no such checking is desired
Y = an array (REAL*8) containing the data
YY = a work array (REAL*8) where
EPSIL = a work array (REAL*8)
To maximize the conditional likelihood, issue the usual call to OPT with ARMAPQ as the name of the function subroutine (ARMAPQ, as well as the METHOD of optimization, must be declared in an EXTERNAL statement):
CALL OPT(P,NPE,F,METHOD,ITERL,MAX,IER,ACC,ARMAPQ,ALABEL)
SECTION 18.2 AUTOCOVARIANCE FUNCTION
To compute the autocovariance function or autocorrelation function of a set of time series observations, use
CALL AUTOCOV(X,COV,N,M,IER,IAUT)
where
X = the (REAL*8) array of data dimensioned X(0:N)
COV = the (REAL*8) array of autocovariances or autocorrelations
computed
N = the number of data points
M = the maximum lag desired (.LT. N - 1)
IER = error return (= -3 input error, = -10 variance is zero)
IAUT = 0 compute autocovariances, = 1 compute autocorrelations
The results are returned in the array COV. COV(0) will contain the variance of the series. In addition, AUTOCOV generates a MATHEMATICA® program in the file AUTOCOV.PR which displays the autocovariance (autocorrelation) function on the computer screen. This file is written to UNIT 16 and the open statement is
OPEN (UNIT=16,FILE='AUTOCOV.PR',STATUS= 'UNKNOWN')
and should not be used for other purposes. Start up MATHEMATICA, to which MATHEMATICA would normally respond by typing on the screen the prompt
In[1]:=
The user should now direct MATHEMATICA to take its input from the file AUTOCOV.PR, which is accomplished by replying to the prompt as follows:
In[1]:=<<autocov.pr
The autocovariance (autorrelation) function will then be displayed.
************************************************************
MATHEMATICA may be obtained from Wolfram Research, Inc.
100 Trade Center Drive, Champaign, IL 61820-7237
************************************************************
A test problem for AUTOCOV is XTIMSER.FOR which uses data in the file LYNX.DAT from M. J. Campbell and A. M. Walker, "A Survey of Statistical Work on the Mackenzie River Series of Annual Canadian Lynx Trappings for the Years 1821-1934 and a New Analysis," J. of the Royal Satistical Society, Series A (1977), 140, Part 4, pp. 411-431.
SECTION 18.3 FOURIER TRANSFORMS
See, for example, Peter Bloomfield, Fourier Analysis of Time Series: An Introduction, Wiley, 1976 and James S. Walker, Fast Fourier Transforms, CRC Press, 1991.
In order to compute the fast Fourier Transform of a (possibly) complex series, use
CALL CXFFT(M,IR,F,G,SN,MM1,MM2,IER)
where
M = the lowest power of 2 at least as great as the number of
data points
IR = the exponent of 2 such that M = 2**IR
MM1 = M - 1
MM2 = M/4 - 1
F = a (REAL*8) array, dimensioned F(0:MM1), containing the real
parts of the series to be fast fourier transformed
G = a (REAL*8) array, dimensioned G(0:MM1), containing the
imaginary parts of the series
SN = a (REAL*8) work array dimensioned SN(0:MM2)
IER = error return = -3 if IR <= 3
The real and imaginary parts of the Fourier Transform are returned in F amd G respectively. Note that to transform a real series, the elements of G should be set equal to zero. Note also that if the actual number of data points, say N, is less than M, the last M - N elements of F and G should be set equal to zero. The inverse fast Fourier Transform is computed by
CALL CXINVFFT(M,IR,F,G,SN,MM1,MM2)
where the meaning of the arguments is the same as for CXFFT, except that F and G contain on input the real and imaginary components of a series the inverse Fourier Transform of which is desired and contain on output the Inverse Fourier transform.
The compute the periodogram of a time series, use
CALL PERIOD00(X,G,PER,MM1,MM2,N,IR)
X = a (REAL*8) array, dimensioned X(0:MM1), containing the
real parts of the series
G = a (REAL*8) array, dimensioned G(0:MM1), containing the
imaginary parts of the series
PER = a (REAL*8) array in which the periodogram is returned
N = the number of data points
IR = the exponent of 2 such that M = 2**IR is the lowest power of
2 that is at least as great as N
MM1 = M - 1
MM2 = M/4 - 1
Before calling PERIOD00, it is useful to
CALL DATACHK(X,G,M,MM1,MM2,IR,N,IER)
where
M = the smallest value of 2**IR (IR integer) such that M >=
N which takes X, G, M, and N as its inputs and computes the
appropriate values of IR, MM1, and MM2.
IER is an error return and returns -3 if N > 2**20 (thus
limiting the length of a series to slightly over one million
datapoints). It also checks whether the value of M input into the
routine is correct.
In addition, PERIOD00 generates a MATHEMATICA program in the file PERGRAM.PR which displays the periodogram on the user's computer screen. This file is written to UNIT 14 and the open statement is
OPEN (UNIT=14,FILE='PERGRAM.PR',STATUS= 'UNKNOWN')
and should not be used for other purposes. Start up MATHEMATICA, to which MATHEMATICA would normally respond by typing on the screen the prompt
In[1]:=
The user should now direct MATHEMATICA to take its input from the file PERGRAM.PR, which is accomplished by replying to the prompt as follows:
In[1]:=pergram.pr
A test problem for PERIOD00 is XTIMSER.FOR which uses data in the file LYNX.DAT from M. J. Campbell and A. M. Walker, "A Survey of Statistical Work on the Mackenzie River Series of Annual Canadian Lynx Trappings for the Years 1821-1934 and a New Analysis," J. of the Royal Satistical Society, Series A (1977), 140, Part 4, pp. 411-431.
SECTION 18.5 SPECTRAL ANALYSIS
See M. B. Priestley, Spectral Analysis and Time Series, Academic Press, 1981. To compute the spectrum obtained by the Fourier Transform of the periodogram, use
CALL SPEC00(N,M,MM1,MM2,IR,WIN,X,G,PER,SN,IER)
where
N = number of data points
M = 2**IR, where IR is the smallest integer such that N <= M
MM1 = M - 1
MM2 = M/4 - 1 IR = smallest integer such that 2**IR >= N
WIN = a CHARACTER*1 variable set to 'p' for a Parzen window, 'h'
for a Hanning window, 'w' for a Welch window, '4' for a Welch-type
window with exponent of 4, '6' for a Welch-type window with
exponent of 6
X = a REAL*8 array dimensioned X(0:MM1) containing the data
G = a REAL*8 work array dimensioned G(0:MM1)
PER = a REAL*8 array dimensioned PER(0:MM1) which will contain
the output, i.e., the ordinates of the spectrum
SN = a REAL*8 work array dimensioned SN(0:MM2)
IER = an error code = -3 if IR< 3
In addition, SPEC00 generates a MATHEMATICA program in the file SPECTRUM.PR which displays the spectrum on the user's computer screen. This file is written to UNIT 13 and the open statement is
OPEN (UNIT=13,FILE='SPECTRUM.PR',STATUS= 'UNKNOWN')
and should not be used for other purposes. Start up MATHEMATICA, to which MATHEMATICA would normally respond by typing on the screen the prompt
In[1]:=
The user should now direct MATHEMATICA to take its input from the file SPECTRUM.PR, which is accomplished by replying to the prompt as follows:
In[1]:=<<spectrum.pr
To compute the spectrum obtained by the Cosine Transform of the autocovariance function, use
CALL SPEC01(F,G,COV,TCOV,m,m2,win,ier)
where
F = a REAL*8 array of data dimensioned F(0:M)
G = a REAL*8 array dimensioned G(0:M) in which the
COV = a REAL*8 array dimensioned COV(0:M2) which will contain the
autocovariances
TCOV = a REAL*8 array dimensioned TCOV(0:M2) which will contain
the windowed autocovariances
M = the number of observations
M2 = the number of lags beyond which window assigns zero weight
to autocovariances. If M2 = 0, this is computed automatically,
otherwise the M2 value in the call is used
WIN = CHARACTER*1 variable, if 'p' use Parzen window, if 'h' use
Hanning window, if 'b' use Bartlett window
IER = error return; 0 if normal, -10 if series is a constant
A test problem for SPEC01 is YTIMSER.FOR which uses data in the file LYNX1.DAT from M. J. Campbell and A. M. Walker, "A Survey of Statistical Work on the Mackenzie River Series of Annual Canadian Lynx Trappings for the Years 1821-1934 and a New Analysis," J. of the Royal Satistical Society, Series A (1977), 140, Part 4, pp. 411-431. Note that a peak in the spectrum at a point k on the horizontal axis reflects a frequency of 2*PI*k/m.
Return to
|Sect. 19.1|Sect.
19.2|Beginning|