*
* @SSMSpectrum(A=A matrix,SW=SW matrix,other options) spectrum
*
* Returns the estimated multivariate spectrum from a state space model
* given by the A, F and SW options. It returns in its argument a series
* of complex matrices.
*
* Parameters:
*   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. In order to avoid problems
*     with non-stationary models, the 0 frequency isn't calculated.
*
* Options:
*   A =transition matrix for the state space model
*   F =mapping from fundamental shocks to shocks to the states [not used]
*   SW=covariance matrix of the fundamental shocks.
*
*   These are the same options as would be used in a DLM instruction.
*
*   ORDINATES=# of frequency ordinates on which calculations are based [512].
*    This will return values only out to frequency pi, thus to
*    entries ORDINATES/2+1.
*
*   COMPONENTS=VECT[INT] of components (numbered between 1,...,# of
*    states) which are actually wanted. For instance, COMPONENTS=||1|| will
*    return the estimate just for the first component.
*
* Revision Schedule:
*   09/2004 Written by Tom Doan, Estima
*   12/2004 Added FACTOR option
*   09/2010 Add F option
*
procedure SSMSpectrum spectrum
type series[cmatrix] *spectrum
*
option rect a
option symm sw
option rect f
*
option integer ords 512
option vec[integer] components
option series[cmatrix] *factor
*
local integer 	i j m n p t last
local cmatrix 	atilde ainv
local rect 		selector
local symm 		fswf
*
if .not.%defined(a).or..not.%defined(sw).or..not.%defined(spectrum) {
   disp "###@SSMSpectrum(A=A matrix,SW=SW matrix,other options) spectrum"
   return
}

compute n =%rows(a)
compute p =n
if %defined(f) {
   if %rows(f)<>n {
      disp "###@SSMSpectrum: dim(F)=" %rows(f) "Doesn't Match dim(A)=" n
      return
   }
   compute p=%cols(f)
}
if %rows(sw)<>p {
   disp "###@SSMSpectrum: dim(SW)=" %rows(sw) "Should Be" p
   return
}

if %defined(components) {
   compute m=%rows(components)
   dim selector(m,n)
   ewise selector(i,j)=(j==components(i))
}
*
compute last=ords/2+1
gset spectrum 1 last = %zeros(n,n)
if %defined(factor)
   gset factor 1 last = %zeros(n,n)
*
if %defined(f)
   compute fswf=f*sw*tr(f)
else
   compute fswf=sw
*
do t=2,last
   compute atilde=%identity(n)-%unit2(2-t,ords)*a
   compute ainv=%cxinv(atilde)
   if %defined(components)
      compute ainv=selector*ainv
   compute spectrum(t)=ainv*fswf*%cxadj(ainv)
   if %defined(factor)
      compute factor(t)=ainv*%decomp(fswf)
end do t
end

