**************************************************************
*PROCEDURE SPECTH invec specs freqs
*
* NOTE: Requires separate ACF.SRC procedure!
*
* Parameters:
* invec Vector of AR and MA coefficients
* specs The series to which the calculated
* spectrum is saved
* freqs The frequencies for which the spectrum
* is calculated is saved (equally
* spaced from 0 to PI, inclusive)
*
* That is, specs(k) is the spectrum calculated at
* freqency freqs(k)
*
* Note, the freq(k) are the freqencies. The series pind
* is also created which contains the freqencies divided by
* 2*PI, thus going from 0 to 1/2.
*
*
* Options
* NCOV 200 The number of covariances used to
* calculate the spectrum
* NFREQ 201 The number of freqnices over which
* the spectrum is calculated
* (using a round number + 1 makes the
* frequencies more interpretable)
* MA 0 The number of moving average parameter
* in the input vector
* LOGS/[NOLOGS] Plot the spectrum or log spectrum
* BDROP 0 The number of entries to drop from the
* EDROP 0 beginning and end of the input vector,
* as, for example, the input vector might
* be created by LINREG or BOXJENK and the
* first couple or last few entries might
* be CONSTANTs, TRENDS, DUMMY variables,
* or other regressors that one would want
* to exclude.
*
* SPECTH.SRC takes a vector of AR and MA coefficients and
* calculates and graphs the theoretical spectrum (or log-spectrum).
* It first uses ACF.SRC to calculate the autocovariances of the
* process given the AR and MA coefficients and uses them to
* calculate the spectrum. Since the sum should be over the
* autocovarinaces from 0 to infinity, and the procedure uses them
* from 0 to ncov, the spectrum will not be exact. The approximation
* is less close the smaller ncov is and the more persistent the
* ARMA process is. But using a couple hundred autocovariances
* usually makes the numbers close; if you're concerned 500 or so
* will make the spectrum close to many decimal places, but increases
* computing time.
*
*
* NOTES:
*
* THE AR PARAMETERS MUST PRECEDE THE MA PARAMETERS IN THE INPUT VECTOR
*
* AND
*
* THE SERIES MUST NOT HAVE ROOTS ON OR INSIDE THE UNIT CIRCLE OR THE
* PROCEDURE WILL YIELD NONSENSE, i.e, IT MUST BE STATIONARY.
*
*
* If you want to find the spectrum of a seasonal process, instead
* of mulitplying out the lag polynomials by hand, you can create
* the lag polynomials of the AR part and the SAR part separately
* and multiply them using POLYMULT.SRC (with the default NOZEROS
* and MINUS options). One can similarly do that for an MA and SMA
* part. Then concatenate the two products into one vector an input
* them into SPECTH.
*
* An example: suppose you want te spectrum of the Seasonal ARMA
* model (2,0,1)(1,0,1) for seasonal period of 12 and coeffients:
* AR(1), AR(2), -.3 .2
* SAR(1), .8
* MA(1), .5
*
****
* ALL 250
* DECL VECT a(2) b(1) c(1)
* INPUT(unit=keyboard) a b c
* -.3 .2
* .8
* .5
*
* SOURCE(NOECHO) POLYMULT.SRC
* @POLYMULT(span=12) a b ab
* DECL vect invec(15)
* DO i = 1,14
* COMPUTE invec(i) = ab(i)
* END DO
* COMPUTE invec(15) = c(1)
* SOURCE(noecho) ACF.SRC
* SOURCE(noecho) SPECTH.SRC
* @SPECTH(nfreq=201,ncov=250,ma=1) invec specs freqs
****
*
* References:
* Granger, C.W.J, and P. Newbold [1986]: Forecasting Economic
* Time Series (Second Edition). Academic Press, San Diego, CA
* Hamilton, J. [1994]: Time Series Analysis. Princeton University
* Press, Princeton, NJ
*
* Norman Morin, nmorin@frb.gov
* v1, October 1997
*
**************************************************************
****************************************************************************
PROCEDURE SPECTH invec specs freqs
TYPE VECT INVEC
TYPE SERIES *SPECS *FREQS
OPTION INTEGER ncov 20
OPTION INTEGER nfreq 201
OPTION INTEGER ma 0
OPTION SWITCH LOGS 0
OPTION INTEGER BDROP 0
OPTION INTEGER EDROP 0
LOCAL VECT invc
LOCAL SERIES acvs ars mas
LOCAL INTEGER i j k p q r
DIM invc(%rows(invec)-bdrop-edrop)
DO i = 1,%rows(invc)
COMPUTE invc(i) = invec(i+bdrop)
END DO
CLEAR ind pind
COMPUTE p = %rows(invc)-ma
COMPUTE q = ma
COMPUTE r = %IMAX(p,q)
SET ars 1 r+1 = 0.
COMPUTE ars(1) = -1.
SET mas 1 r+1 = 0.
COMPUTE mas(1) = -1.
IF p.gt.0 {
DO i = 1,p
COMPUTE ARS(i+1) = invc(i)
END DO
}
IF q.gt.0 {
DO i = 1,q
COMPUTE MAS(i+1) = -1.*(invc(i+p))
END DO
}
DISP ' '
DISP 'Plotting the spectrum for an ARMA(' @-1 p @-1 ',' @-1 q @-1 ')'
DISP ' over' NFREQ 'frequencies using' NCOV 'covariances.'
DISP
DISP ' y(t) = a1*y(t-1) + ... + ap*y(t-p) + e(t) + b1*e(t-1) + ..+ bq*e(t-q)'
DISP ' Lag' @13 'a(i) b(i)'
DO i = 1,%IMIN(p,q)
DISP #### i @10 ##.##### ars(i+1) @23 ##.##### -1*mas(i+1)
END DO
IF p.gt.q {
DO i = q+1,p
DISP #### i @10 ##.##### ars(i+1)
END DO
}
IF q.gt.p {
DO i = p+1,q
DISP #### i @23 ##.##### -1*mas(i+1)
END DO
}
DISP
@ACF(COV,ncorrs=ncov,noprint,nograph,ma=ma) invc acvs
SET spec 1 nfreq = 0.
SET ind 1 nfreq = (t-1)*%pi/(nfreq-1)
SET pind 1 nfreq = .5*ind/%pi
DO i = 1,nfreq
COMPUTE spec(i) = acvs(1)
DO j = 1,ncov
COMPUTE spec(i) = spec(i) + 2*acvs(j+1)*cos(ind(i)*j)
END DO
COMPUTE spec(i) = spec(i)/(2*%pi)
END DO
IF LOGS.eq.1
SET specs 1 nfreq = log(SPEC)
ELSE
SET specs 1 nfreq = spec
SET freqs 1 nfreq = ind
*
SCATTER(style=lines) 1 ; #pind specs
*
END