/*


                             multir2

   finds implied r2 as a function of horizon from orthogonalized
   impulse response function.

   x(t) = D(L)e(t)

   x(t)+x(t+1) = (D1 e(t-1) + D2 e(t-2) + ... ) + (D2 e(t-1) + D3 e(t-2) + ...)
                 + D0 e(t) + D0 e(t+1) + D1 e(t).

   x(t)+x(t+1) = (D1+D2) e(t-1) + (D2+D3) e(t-2) +...+ (D0+D1) e(t) + D0 e(t+1)

   E(e e') = I

   R2(first variable) = (D1+D2)(D1+D2)' + (D2+D3)(D2+D3)' + ...  (11 element)
                        ------------------------------------------------------
                           (all that + D0D0' + (D0+D1)(D0+D1)' ) (11 element)

*/

proc(2) = multir2(dxir,max);
clearg j,dsum,topmat,i,topmat2,i,botmat,botmat2,r2c,r2y;
r2c = zeros(max,1);
r2y = zeros(max,1);
j = 1;
do while j <= max;
    dsum = sumc(dxir[2:2+j-1,.])';
    topmat = reshape(dsum,2,2)';
    topmat = topmat*topmat';                  /* (D1+D2+..)(D1+D2+..)' */
    i = 1;
    do while 2+i+j-1 <= rows(dxir);
        dsum = sumc(dxir[2+i:2+i+j-1,.])';
        topmat2 = reshape(dsum,2,2)';
        topmat  = topmat + topmat2*topmat2';  /* (D2+D3..)(D2+D3+..)' */
        i = i+1;
    endo;
    dsum = dxir[1,.];
    botmat = reshape(dsum,2,2)';
    botmat = botmat*botmat';
    i = 2;
    do while i <= j;
        dsum = sumc(dxir[1:i,.])';
        botmat2 = reshape(dsum,2,2)';
        botmat = botmat + botmat2*botmat2';
        i = i+1;
    endo;
    botmat = botmat+topmat;
    r2c[j] = topmat[1,1]/botmat[1,1];
    r2y[j] = topmat[2,2]/botmat[2,2];
    j = j+1;
endo;
retp(r2c,r2y);
endp;

proc(1) = lhvardec(dxir,max);
local answer,j,bigdx,bigmat,i;
answer = zeros(max,4);
j = 1;
do while j <= max;
    if j == 1;
        bigmat = dxir;
    endif;
    if j > 1;
        bigdx = zeros(j-1,4)|dxir;
        bigmat = bigdx;
        i = j;
        do while i <= rows(bigmat);
            bigmat[i,.] = sumc(bigdx[i-j+1:i,.])';
            i = i+1;
        endo;
    endif;
    answer[j,.] = sumc(bigmat^2)';
    j = j+1;
endo;
retp(answer);
endp;


/* test calls */

output file=multir2.out reset;
max = 20;
{r2c , r2y} = multir2(dxir,max);

table = "horiz"~"R2c"~"R2y";
table = table|(seqa(1,1,max)~r2c~r2y);
mask = (0~0~0)|ones(rows(table)-1,3);
fmt = "*.*lf"~9~2;
call printfm(table,mask,fmt);



answer = lhvardec(dxir,max);
"variance decomposition by horizon";
"totals";
"";
table = ("horiz"~"c,c shok"~"y,c shok"~"c,y shok"~"y,y shok")|
        (seqa(1,1,max)~answer);
mask = zeros(1,5)|ones(rows(answer),5);
fmt = "*.*lf"~12~3;
call printfm(table,mask,fmt);

"";
"percents";
"";
total = (1~1).*.((answer[.,1]+answer[.,3])~(answer[.,2]+answer[.,4]));
table = ("horiz"~"c,c shok"~"y,c shok"~"c,y shok"~"y,y shok")|
        (seqa(1,1,max)~(100*answer./total));
mask = zeros(1,5)|ones(rows(answer),5);
fmt = "*.*lf"~12~1;
call printfm(table,mask,fmt);

output off;