/* ------------------------------------------------------------------------ */
/*                               finduniv.g                                 */
/* finds univariate impulse-response function implied by a var.             */
/* does it by factoring spectral density                                    */
/*
    Usage:
        {lagpoly,z} = finduniv(dxir,varnum,y);
    Inputs:
        dxir: multivariate impluse-response to orthogonalized shocks, as
              produced by var.g
        varnum: which variable to calculate univariate i-r: usually 2
    Outputs:
        lagpoly: lag polynomial = impulse response itself (all real)
        z: beveridge-nelson trend to series y.

   How the program works.
   start with an orhtogonalized wold representation
        Dx(t) = D(L)v(t) E(vv') = I.
   the program receives dxir which are
        dxir = d11(L)~d21(L)~d12(L)~d22(L).
   It picks out
        dofl = d11(L)~d12(L).
   The spectral density of Dx1 is d11d11* + d12d12*. thus we have to find
        a(z)a(z^-1) = d11(z)d11(z^-1) + d12(z)d12(z^-1)
   it does it by finding the roots of the right hand side.  specifically,
        d(z)d(z^-1) = (d0 + d1z + d2z^2 + .. dkz^k)(d0 + d1z^-1 + .. )
                    = (d0^2 + d1^2 ..+dk^2) +
                      (d0d1 + d1d2 ..) (z + z^-1) + ...
                      (d0dk)(z^k + z^-k)
   thus, the first step is construct the matrix poly, of coefficients.
        poly  =  [ d0^2 + d1^2 ..+dk^2 ]   + (same for d12(L))
                 [ d0d1 + d1d2 ..      ]
                 [  ..                 ]
                 [ d0dk                ]

   Next, convert to the polynomial of z^k*d(L).  this is

                 [ d0dk        ]
                 [ ..          ]
                 [ d0d1 + ..   ]
                 [ d0^2 + ...  ]
                 [ d0d1 + ..   ]
                 [ ..          ]
                 [ d0dk        ]
   Find roots of this polynomial, r.  half are more than one, and these are
   the roots of a(L).  It then creates a(L) from its roots using polymkjc.

    Beveridge nelson decomposition.  One possibility is to invert to an ar
    representation.  Instead, the program sets pre-sample errors and growth
    rates to zero.  Then, recursively,

    e(t) = Dy(t) - [ c1 e(t-1) + c2 e(t-2) + ... ]
    z(t) = y(t) + [ c1+c2+..] e(t) + [c2+c3+..]  e(t-1) + ...
*/
/* ------------------------------------------------------------------------ */

/* proc to make polynomials from roots.  Gauss's polymake adpted to complex
   roots.  the only difference is r is Nx2 and the multiplication is complex.
   someday make this smarter to exploit the fact that you konw roots come
   in pairs and hence the answer must be real                              */

proc polymkjc(r);
    local n,c,j;
    n = rows(r);
    c = zeros(n+1,2);
    j = 1;
    c[1,.] = 1~0;
    do until j > n;
        c[2:j+1,.] = c[2:j+1,.] -
                    ( ( r[j,1]*c[1:j,1] - r[j,2]*c[1:j,2] ) ~
                      ( r[j,2]*c[1:j,1] + r[j,1]*c[1:j,2] ) );
        j = j+1;
    endo;
    retp(c);
endp;

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


proc(2) = finduniv(dxir,varnum,y);

local nseries, horizl, dofl, poly, lagpoly, lagpoly2, i, bigd, r, r2, test,
       cstar, e, dy, t, z;

    nseries = cols(dxir)^(1/2);
    horizl = rows(dxir);

    /* part 1: find lagpoly, lag polynomial of univariate representaiton */

    dofl = dxir[.,seqa(varnum,nseries,nseries)];   /* pick out varnum's resp */
    poly = zeros(horizl,1);
    i = 1;
    do while i <= cols(dofl);
        bigd = ones(horizl,1)*(dofl[.,i]');
        bigd = (shiftr(bigd,-seqa(0,1,horizl),0));
        poly = poly + bigd*dofl[.,i];
        i = i+1;
    endo;
    poly = rev(poly)|poly[2:rows(poly)];
    r = polyroot(poly);
    r2 = selif(r,(r[.,1]^2+r[.,2]^2) .> 1 ); /* choose roots outside unit o */
    if rows(r2) /= (rows(r)/2); "Finduniv.g: lost roots or unit roots"; endif;
    lagpoly = polymkjc(r2);
    test = maxc(abs(lagpoly[.,2])) ;
    if test > 1E-5;
        "Finduniv.g: complex coeffs in lag polynomial";
        "Largest imaginary part is " test "and is ignored";
    endif;
    lagpoly2 = rev(lagpoly[.,1]);         /* real part from 0 to k not k to 0 */
    lagpoly2 = lagpoly2./lagpoly2[1];       /* normalize to 1st element = 1 */


    /* part 2: find beveridge-nelson trend */

    /* timing:  1 2 .... horizl horizl+1 ... horizl+rows(y)-1 */
    /*                    y[1]                   last y      */
    /*          0 0               1stdy          last dy     */

    cstar = rev(cumsumc(rev(lagpoly2)));
    cstar = cstar[2:rows(cstar)]|0;          /* coeffs of b-n stat component */

    e = zeros(horizl+rows(y)-1,1);           /* set up errors incl. presample */
    dy = e ;
    z = e;
    dy[horizl+1:rows(dy)] = y[2:rows(y)]-y[1:rows(y)-1];
    dy[horizl+1:rows(dy)] = dy[horizl+1:rows(dy)]-meanc(dy[horizl+1:rows(dy)]);
    z[horizl] = y[1];
    t = horizl+1;
    do while t <= horizl+rows(y)-1;
        e[t] = dy[t] -
               sumc(lagpoly2[2:rows(lagpoly2)].*rev(e[t-rows(lagpoly2)+1:t-1]));
        z[t] = y[t-horizl] + sumc(cstar.*rev(e[t-rows(cstar)+1:t]));
        t = t+1;
   endo;
   z = z[horizl:rows(z)];

   retp(lagpoly2,z);

endp;

/* -------------------------------------------------------------------------- */
/* test call */

/*

{lagpoly,z} = finduniv(dxir,varnum,y);
_pltype = 6;
xy(seqa(1,1,10),lagpoly[1:10]);
xy(seqa(47.25,1,rows(y)),y~z);

*/