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

Y_{T }= X(1)+X(2)*Y_{T-1}+...+X(1+IP)*Y_{T-IP}+U_{T}+X(2+IP)*U_{T-1}+ +...+X(1+IP+IQ)*U_{T-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|