/* ------------------------------------------------------------------------- */
/*                                 CALCSD.G                                  */
/* proc calcsd: caclulates spectral density from wold representation         */
/*
   Usage:
    {w,sd} = calcse(w,dxir);
   Inputs:
    w: where to calcluate (how fine a grid) in factions of pi.
       -1: default, w=seqa(0,0.01,100);
    dxir: impulse response to unit shocks of orthogonalized wold representation.
          organized as
                      response of
                    x1 x2 x3...  x1 x2 x3...
                      to shock in
                    x1 x1 x1     x2 x2 x2..
          if Wold is
                    x(t) = C(L)e(t)  E(ee') = sigma = r'r
          (r is choleski)
          orthogonlaized wold is
                    e(t) = r'v(t)   E(vv') = I
                    x(t) = C(L)r'v(t) = D(L)v(t).
          the impluse responses dxir should be the elements of D(L) not C(L).
                                                                              */
/* -------------------------------------------------------------------------- */


proc(2) = calcsd(w,dxir);
    local cosw,sinw,horiz,nseries,numw,j,sd,temp;

    if w == -1; w = seqa(0,.01,100); endif;
    w = w*pi;
    horiz = rows(dxir);
    nseries = cols(dxir)^(1/2);
    numw = rows(w);
    sd = zeros(numw,nseries);

/* -------------------------------------------------------------------------- */
/*It uses
    d(e-iw) d(eiw) = (d0 + d1cosw + d2cos2w +..)^2 +(d1sinw + d2sin2w +...)^2
  then from Dx = D(L)v, sd's are
    sx1 = d11^2+d12^2 and sx2 = d21^2 + d22^2
  program first forms temp'= [ d11^2 d12^2 ]
                             [ d21^2 d22^2 ]
  then sums to get spectral density at each frequency.  It does NOT calculate
  cross spectral density.  */
/* -------------------------------------------------------------------------- */

    j = 1;                                       /* find sd at each w         */
    do while j <= rows(w);
       cosw = cos(seqa(0,1,horiz)*w[j]);         /* 1 cosw cos2w...           */
       sinw = sin(seqa(0,1,horiz)*w[j]);

       temp = (sumc(dxir.*cosw))^2 + (sumc(dxir.*sinw))^2;
                                                  /* [d11^2 d21^2 d12^2 d22^2 */
       temp = reshape(temp,nseries,nseries);
       sd[j,.] = (sumc(temp))';
       j = j+1;
    endo;
    retp(w,sd);
endp;

/* --------------------------------------------------------------------------*/


/* test call:

library pgraph;
graphset;

{w,sd} = calcsd(-1,dxir);

_Pdate = "";
_Pltype = 6;
xy(w,sd);

*/