*
* @EqnToACF(options)  equation  acf
*
* Creates a series with the theoretical autocovariance/autocorrelation
* function for the ARMA model described by an equation. For the algorithm used,
* see Brockwell&Davis, Introduction to Time Series and Forecasting, 2nd ed, section
* 3.2. This uses "Third Method" from page 91.
*
* Parameters:
*
*   equation=input equation. This can be produced by the DEFINE option on BOXJENK or
*    LINREG, or it could be input directly using EQUATION.
*
*   acf=output series of covariances (assuming unit innovation variance) or correlations
*
* Options:
*
*   number=number of lags to compute. By default, 1 more than the maximum AR or MA lag.
*   correlations/[nocorrelations]  By default, the procedure returns autocovariances. To
*     get autocorrelations (which will have a 1.0 in the first slot, since that's actually
*     lag zero), use the correlations option.
*
*  Example:
*    This produces the first 10 autocovariances from the ARMA(1,1) process,
*     y(t)=.8y(t-1)+u(t)+.3u(t-1). Note that a unit variance for u is assumed.
*
*    equation(noconst) arma11 y
*    associate arma11
*    # .8 .3
*    @eqntoacf(number=10) arma11 armaacf
*
*  Revision Schedule
*    5/03 Written by Tom Doan, Estima
*
procedure EqnToACF eq acf
type equation eq
type series *acf
*
option integer number
option switch correlations 0
*
local rect[int] eqtable
local vect eqcoeffs
local integer i j l p q r
local series psi
local vector theta phi psitheta base
local rect cx
local integer n
local real v
*
compute eqtable=%eqntable(eq),eqcoeffs=%eqncoeffs(eq)
compute q=p=0
do i=1,%eqnsize(eq)
   if eqtable(1,i)==%eqndepvar(eq).and.eqtable(2,i)>p
      compute p=eqtable(2,i)
   if eqtable(1,i)==%mvgavge.and.eqtable(2,i)>q
      compute q=eqtable(2,i)
end do i
*
compute r=%imax(p+1,q+1)
if %defined(number).and.number>r
   compute n=number
else
   compute n=r

impulse(noprint,shocks=||1.0||) 1 r
# eq psi
*
dim theta(r) phi(r)
ewise theta(i)=(i==1)
ewise phi(i)=(i==1)
*
do i=1,%eqnsize(eq)
   if eqtable(1,i)==%eqndepvar(eq)
      compute phi(eqtable(2,i)+1)=-eqcoeffs(i)
   else
   if eqtable(1,i)==%mvgavge
      compute theta(eqtable(2,i)+1)=eqcoeffs(i)
end do i
*
dim psitheta(r) cx(r,r)
do i=1,r
   ewise psitheta(i)=%dot(%xsubvec(theta,i,r),%xsubvec(psi,1,r+1-i))
end do i
*
compute cx=%const(0.0)
do i=0,r-1
   do j=0,r-1
      compute l=%iabs(i-j)
      compute cx(i+1,l+1)=cx(i+1,l+1)+phi(1+j)
   end do j
end do i
compute base=%solve(cx,psitheta)
set acf 1 r   = base(t)
set acf r+1 n = v=0.0,%do(i,2,r,v=v-phi(i)*acf{i-1}),v
*
if correlations
   set(scratch) acf 1 n = acf(t)/acf(1)
end eqntoacf



