*
* @VARSpectrum  model sigma spectrum
*
* Returns the estimated multivariate spectrum from a VAR given by the
* combination of model and sigma. It returns in its third argument a series of
* complex matrices.
*
* Parameters:
*   model  = input model (generally created by SYSTEM and estimated using
*     ESTIMATE)
*   sigma  = covariance matrix of the VAR
*   spectrum [output] = SERIES[CMATRIX] with the estimated spectral
*     density. 0 frequency will be in entry 1; in general, frequency
*     2xpix(t-1)/ORDINATES will be in entry t.
*
* Option:
*   ORDINATES=# of frequency ordinates on which calculations are based [256].
*     This will return values only out to frequency pi, thus to entries ORDINATES/2+1.
*   FACTOR=(output) SERIES[CMATRIX] which has a frequency by frequency (complex matrix)
*     factorization of the spectral density matrix.
*
* Revision Schedule:
*   09/2004 Written by Tom Doan, Estima
*   12/2004 Added factor option
*
procedure VARSpectrum model sigma spectrum
type model model
type symm sigma
type series[cmatrix] *spectrum
*
option integer ords 256
option series[cmatrix] *factor
*
local integer i j k n last
local vec[integer] mdepvar
local series[cmatrix] atilde
local equation eqn
local rect[int] etable
local vect ecoeff
local rect comat
local cmatrix eye
*
compute n      =%modelsize(model)
compute mdepvar=%modeldepvars(model)
compute comat  =%zeros(n,n)
if %rows(sigma)<>n {
   disp "VARSpectrum Error. dim(sigma)=" %rows(sigma) "Doesn't Match Model Size" n
   return
}
*
dim eye(n,n)
ewise eye(i,j)=(i==j)
compute last=ords/2+1
gset atilde 1 last = eye
*
do i=1,n
   compute eqn   =%modeleqn(model,i)
   compute etable=%eqntable(eqn)
   compute ecoeff=%eqncoeffs(eqn)
   do j=1,%eqnsize(eqn)
      do k=1,n
         if etable(1,j)==mdepvar(k) {
            compute comat(i,k)=ecoeff(j)
            gset atilde 1 last = atilde(t)-%unit(-2*%pi*(t-1)*etable(2,j)/ords)*comat
            compute comat(i,k)=0.0
         }
      end do k
   end do j
end do i
gset spectrum 1 last = %cxinv(atilde)*sigma*%cxadj(%cxinv(atilde))
if %defined(factor)
   gset factor 1 last = %cxinv(atilde)*%decomp(sigma)
end

