*************************************************************
*
* PROCEDURE ACF invec corrs
*
*  Parameters:
*     invec     Vector of AR and MA parameters
*     corrs     The name of the series to which the calculated
*                  autocorrelations/autocovsariances are saved
*
*   Options
*     COV/[NOCOV]         Produce autocovariances or autocorrelations
*     MA      [0]         The number of MA lags
*     NCORRS [12]         The number of autocorrelations
*                           to calculate
*     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.
*     [PRINT]/NOPRINT     Print the calculated autocorrelations
*     [GRAPH]/NOGRAPH     Graph them
*     STYLE = [BAR]/LINES The Style for the graphs
*
*
* Given a vector containing the a(i) and b(i) from the ARMA(p,q) model
*  for y(t) =
*   a1*y(t-1) + ... + ap*y(t-p) + e(t) + b1*e(t-1) + ..+ bq*e(t-q)
* this procedure calculates the theoretical autocorrelations or
* autocovariances and will print and/or graph them.
*
* NOTE:  THE AR PARAMETERS MUST PRECEDE THE MA PARAMETERS
*
*     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.
*
*
*
*   References:
*           McLeod, A. [1975], Applied Statistics, 24: 255-256
*           McLeod, A. [1977], Applied Statistics, 26: 194
*
*   Norman Morin, nmorin@frb.gov
*           v2, October 1997
*           v1, April 1997
*
*
*************************************************************


PROCEDURE ACF invec corrs

TYPE VECTOR invec
TYPE SERIES *corrs

OPTION INTEGER MA 0
OPTION INTEGER NCORRS 12
OPTION INTEGER BDROP 0
OPTION INTEGER EDROP 0
OPTION SWITCH  COV 0
OPTION SWITCH  PRINT 1
OPTION SWITCH  GRAPH 1
OPTION CHOICE  STYLE 1 BAR LINES

LOCAL INTEGER p q ncorr ar r counter i j k ii jj
LOCAL SERIES ars mas c b covs corrs
LOCAL RECT Amat Zm Ainv Covmat Corrmat
LOCAL VECT coeffs

SET corrs 1 ncorrs+1 = 0.
SET covs 1 ncorrs+1 = 0.

DIM coeffs(%rows(invec)-BDROP-EDROP)
compute counter = (%rows(invec)) ; dis 'count = ' counter
dis 'rows,cols of coeffs = ' %rows(coeffs) %cols(coeffs)
DO i = 1,%rows(invec)-BDROP-EDROP
  dis 'i = ' i 'bdrop, edrop = ' bdrop edrop
  COMPUTE coeffs(i) = invec(i+BDROP)
END DO

COMPUTE p = %rows(coeffs)-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) = coeffs(i)
  END DO
  }

IF q.gt.0 {
  DO i = 1,q
    COMPUTE MAS(i+1) = -1.*(coeffs(i+p))
  END DO
 }

SET c 1 r+1 = 0.
COMPUTE c(1) = 1
SET b 1 r+1 = 0.
COMPUTE b(1) = 0.


DO k = 2,r+1
   COMPUTE c(k) = -MAS(k)
   DO i = 1,%IMIN(p,k-1)
      COMPUTE c(k) = c(k) + ars(i+1)*c(k-i)
   END DO
END DO

DO k = 1,q+1
   DO i = k,q+1
      COMPUTE b(k) = b(k) - mas(i)*c(i+1-k)
   END DO
END DO


IF p.eq.0 {
  SET covs 1 q+1 = b
  SET corrs 1 q+1 = covs/covs(1)
 }

ELSE {
  DIM Amat(r+1,r+1) Zm(r+1,1)
  COMPUTE Amat = %CONST(0.)
  DO ii = 1,r+1
     COMPUTE Amat(ii,1) = ars(ii)
  END DO
  DO ii = 1,r+1
    DO jj = 2,r+1
      IF ((ii-jj).ge.0).and.((ii-jj).le.p)  {
         COMPUTE Amat(ii,jj) = ars(ii-jj+1)
         }
      IF ((ii+jj-2).ge.0).and.((ii+jj-2).le.p) {
         COMPUTE Amat(ii,jj) = ars(ii+jj-2+1)+Amat(ii,jj)
         }
    END DO
  END DO
  COMPUTE Zm = %CONST(0.)
  DO ii = 1,r+1
     COMPUTE  Zm(ii,1) = -1.*b(ii)
  END DO
  COMPUTE AINV = inv(AMAT)
  COMPUTE COVMAT = Ainv*Zm
  COMPUTE CORRMAT = COVMAT
  EWISE CORRMAT(i,j) = COVMAT(i,j)/COVMAT(1,1)
  SET corrs 1 r+1 = CORRMAT(t,1)
  SET covs 1 r+1 = COVMAT(t,1)
  }

DO k = r+2,ncorrs+1
  DO jj = 1,p
    COMPUTE corrs(k) = corrs(k) + ars(jj+1)*corrs(k-jj)
	 COMPUTE covs(k) = covs(k) + ars(jj+1)*covs(k-jj)
  END DO
END DO

IF COV.eq.1
  SET corrs 1 ncorrs+1 = covs

IF GRAPH.eq.1 {
   GRAPH(number=0,style=STYLE) 1 ; #corrs 1 ncorrs+1
 }

DISP
IF COV.eq.1
  DISP 'The autocovariances for an ARMA(' @-1 p @-1 ',' @-1 q @-1 ')'
ELSE
  DISP 'The autocorrelations for an ARMA(' @-1 p @-1 ',' @-1 q @-1 ')'
DISP
DISP '   y(t) = a1*y(t-1) + ... + ap*y(t-p) + e(t) + b1*e(t-1) + ..+ bq*e(t-q)'
DISP  @10  'a(i)         b(i)'
DO i = 1,%IMIN(p,q)
   DISP @7 ##.##### ars(i+1) @20 ##.##### -1*mas(i+1)
END DO
IF p.gt.q {
  DO i = q+1,p
     DISP @7 ##.##### ars(i+1)
  END DO
  }
IF q.gt.p {
  DO i = p+1,q
     DISP @20 ##.##### -1*mas(i+1)
  END DO
  }
DISP

IF PRINT.eq.1 {
  DISP
  IF COV.eq.1
    DISP '   Lag      Autocovariance'
  ELSE
	 DISP '   Lag      Autocorrelation'
  DO i = 1,ncorrs+1
    DISP  @3 ## i-1   @14 ##.#####  corrs(i)
  END DO
  DISP
 }

END




