/* -------------------------------------------------------------------------- */
/*                                PIHTEST.G

  tests permanent income hypothesis based on univariate AR and VAR.

Usage:
    proc(0) = pihtest(dxirm,dxirum,dxiruu,r,x) ;

Arguments:
    dxirm = multivariate impulse response of differences.
    dxirum = univariate income i-r estimated from VAR
    dxiruu = univariate income i-r estimated from UAR
    r      = a vector of annual percent interest rates to try
              e.g. 5|3 for 5 and 3%
    x = data series, c~e~other.

dxirm is organized as
                     response of
               c  y  other  c  y  other   c      y   other
                        to
               c  c    c    y  y    y    other other other

dxirm and dxiruu are normalized to variance of error = 1;
dxirum is normalized to c(0) = 1, program finds var of error.
this is silly, but keeps the convention of the rest of the program.

Notes & warnings:
    Proc expects consumption , then income orthogonalization.  others will
    not work: it selects consumption as the first variable.  It should work
    with extra right hand variables.


THEORY
  if we write VAR in orthogonalized form,

     [ de(t) ] =  D(L) v(t)      E(vv') = I
     [ dc(t) ]

  the pih predicts

      dc(t) = D11(lambda)ve + D22(lambda)vc.

  and hence

      var(dc) = D11(l)^2 + D22(l)^2

  similarly, if

      de(t) = A(L) e(t)      var(e) = 1

  then

       dc(t) = A(l) e(t)

   so

        var(dc(t)) =  A(l)^2



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


proc(0) = pihtest(dxirm,dxirum,dxiruu,r,x) ;
local cgr, egr, lambda, deirm, sige, dxiruml,
      varpih, i, lamtoj, table, mask, fmt;

    cgr = x[2:rows(x),1]-x[1:rows(x)-1,1];
    egr = x[2:rows(x),2]-x[1:rows(x)-1,2];

    lambda = (1+r/100)^(-1/4);                        /* labmdas conform. to r*/

    deirm = dxirm[.,2+seqa(0,cols(x),cols(x))];       /* find income i-r */

    sige = (stdc(egr)^2)/sumc(dxirum^2);              /* shock variance */
    dxiruml  = dxirum*((sige)^(1/2));                 /* renorm. to var err=1*/

    varpih = zeros(rows(lambda),3);                   /* r x multi unim uniu */

    i = 1;
    do while i <= rows(lambda);
        lamtoj = lambda[i]^seqa(0,1,rows(deirm));
        varpih[i,1] = (sumc(sumc(deirm.*lamtoj)^2))^(1/2);
        varpih[i,2] = sumc(dxiruml.*lamtoj)^2;
        varpih[i,3] = sumc(dxiruu.*lamtoj)^2;
        i = i+1;
    endo;

    /* -------------- */
    /* print table    */
    /* -------------- */

     "Excess sensitivity estimates";
     "      std. dev. of consumption growth" stdc(cgr);
     "           std. dev. of income growth" stdc(egr);
     "        std. dev. of income from VAR "   (sumc(sumc(deirm^2)))^(1/2);
     "std. dev. of income from univ. <- VAR"   sumc(dxiruml^2)^(1/2);
     "std. dev. of income from univ. <-univ"   sumc(dxiruu^2)^(1/2);


     " ";

     table = ("        " ~ "VAR" ~ "univ." ~ "univ." )|
             ("int rate" ~ "VAR" ~ " VAR " ~ "univ." )|
             ( r ~ varpih );
     mask = zeros(2,4)|ones(rows(r),4);
     fmt = "*.*lf"~10~4;
     call printfm(table, mask, fmt);
endp;

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

/* test call */
/*
r = .05|.03|.01;
dxirum = dxir2;
dxiruu = dxir1;
dxirm = dxir;
call pihtest(dxirm,dxirum,dxiruu,r,x) ;
*/