/* ------------------------------------------------------------------------ */
/*                                univar.g                                  */
/* runs Dx on lags of Dx.  finds impulse response function                  */
/* and beveridge nelson decomposition                                       */

/*
   Usage:
     {b1,dxir1,xir1,stat} = univar(x,lags,horiz,prnt);
   Inputs:
     x:         series to do.  levels, e.g. log GNP
     lags:      how many lagged differences to use
     horiz:     how far to do impulse response
     prnt:      1: print out results, 0: just return b and ir fn.
   Returns:
     b1:        regression coefficients.  [const Dx(t-1) Dx(t-2)...];
     dxir1:     impulse response of differences.
     xir1:      impulse response of levels.

    note lag polynomial is normalized to unit shock variance, not c(o) = 1;

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

clearg __output;

proc(4) = univar(x,lags,horiz,names,prnt);

local vanm,m,b1,stb,vc,stderr,sigma,cx,rsq,resid,dwstat,
       b1,dxir1,xir1,T,Dx,rhv,j,lhv,compan,
       amat,cofst,stat;

    T  = rows(x);
    Dx = x[2:T]-x[1:T-1];
    T  = rows(Dx);
    rhv = Dx[lags:T-1];                             /* lags+1:T  = (t)       */
    j = 2;                                          /*  1:T-lags  = (t-lags) */
    do while j <= lags;
        rhv = rhv~Dx[lags+1-j:T-j];
        j = j+1;
    endo;
    lhv = Dx[lags+1:T];
    if prnt; __output = 1; endif;
    if not prnt; __output = 0; endif;
    { vanm,m,b1,stb,vc,stderr,sigma,cx,rsq,resid,dwstat } = OLS("",lhv,rhv);

    /* -------------------------------------------------- */
    /*         simulate to impulse response               */
    /* simulate companion system dy(t) = A dy(t-1)
                                           dy(t-2)
                                            ....          */
    /* -------------------------------------------------- */

   dxir1 = zeros(horiz,1);
   xir1 = zeros(horiz,1);
   compan = zeros(lags,1);
   dxir1[1] = 1*sigma;                   /* feed a 1 sd shock */
   xir1[1] = dxir1[1];
   j = 2;
   do while j <= horiz;
        compan  = dxir1[j-1]|compan;
        compan  = compan[1:rows(compan)-1];
        dxir1[j] = b1[2:lags+1]'compan;
        xir1[j] = xir1[j-1] + dxir1[j];
        j = j+1;
   endo;

   /* ------------------------ */
   /* print impulse response   */
   /* ------------------------ */

   /*
   if prnt;
        "Impulse response function";
        format /RD 10,4;
        seqa(1,1,horiz)~dxir1~xir1 ;
   endif;
   */

   /* ------------------------ */
   /* find b-n trend.          */
   /* ------------------------ */
   /* ------------------------------------------ */
   /* 1) set up big companion  of no-mean system */
   /* dy(t)    =  [  --- b --- ] dy(t-1)
      dy(t-1)     [  1 0 0 0   ] dy(t-2)
       ...        [  0 1 0 0   ]   ....          */
   /* ------------------------------------------ */

    amat = (b1[2:rows(b1)]');
    cofst = 1;
    if lags > 1;
        amat = amat|(eye(lags-1)~zeros(lags-1,1));
        cofst = (1~zeros(1,lags-1));
    endif;
    cofst = cofst*amat*inv(eye(lags)-amat);
    stat  = rhv*(cofst');


retp(b1,dxir1,xir1,stat);
endp;

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


/* test call    */
/*
     output file=trash.out reset;
     let names = "y";
     prnt = 1;
     {b1,dxir1,xir1,stat} = univar(y,lags,horiz,names,prnt);
     output off;
*/