/*
** VAR(1) System for M2, GNP, RTB
** Seemingly Unrelated Regression Estimation
** with impulse response functions
*/

use gpe2;

n=179;
load data[n,7]=tsinghua\usmq.txt;
M1=data[2:n,2];
M2=data[2:n,3];
M3=data[2:n,4];
RTB=data[2:n,5];
GNP=data[2:n,6];
PGDP=data[2:n,7];

RM1=ln((M1./PGDP)*100);
RM2=ln((M2./PGDP)*100);
RM3=ln((M3./PGDP)*100);
RGNP=ln((GNP./PGDP)*100);
RTB=ln(RTB);
/*
3 Equations VAR(1) System
RM2  = RM2(-1) RGNP(-1) RTB(-1) 1
RGNP = RM2(-1) RGNP(-1) RTB(-1) 1
RTB  = RM2(-1) RGNP(-1) RTB(-1) 1
*/
n=rows(rm2);
y=rm2~rgnp~rtb;

call reset;
_names={"RM2","RGNP","RTB"};
_dlags=1;
/*
_dlags={1 0 1,
        1 1 1,
        0 0 1};
@ _id={1 1 -1}; @
*/
_method=3;
_iter=50;
call estimate(y,0);

@ response parameters @
b=__pi[1:3,1:3];

print b;
print impulse(b,0,10);

end;

/* Impulse Response Functions */
@ b1 parameter matrix of lagged endogenous variables @
@ b2 parameter matrix of exogenous variables @
@ b1 must be a square matrix: equations (rows) = endogenous (cols) @
@ row index for # equations, col index for # endogenous variables @
@ impulse response plots include: @
@ (1) plots of all variables' response w.r.t all external shocks @
@ (2) plots of all variables' response w.r.t each individual shock @
proc impulse(b1,b2,lags);
    local n,k,l,i,j,s,bs,bsvec;
    
    n=rows(b1); @ #equations = #endogenous @
    if ismiss(b2) or b2==0; b2=eye(n); endif;
    k=rows(b2);
    l=seqa(0,1,lags+1);
    
    s=ftos(1,"%*.*lf",1,0); @ legend string @
    i=2;
    do until i>k;
        s=s $+ "\000" $+ ftos(i,"%*.*lf",1,0);
        i=i+1;
    endo;

    bsvec=ones(n*k,1);
    bs=b2*b1;
    j=1;
    do until j>lags;
        bsvec=bsvec~vec(bs');
        bs=bs*b1;
        j=j+1;
    endo;
    
    pqgwin many;
    library pgraph;
    graphset;
    _plegctl=1; @ show graph legends @
    _plegstr=s;
    call title("Plot of Impulse Response Functions");
    /*
    call xy(l,bsvec'); @ all IRF plots @
    */

    @ IRF plots of individual shock (for all variables) @
    i=1;
    do until i>k;
        j=(i-1)*n+1;
        call xy(l,bsvec[j:i*n,.]');
        i=i+1;
    endo;

    /*
    @ individual IRF plot @    
    i=1;
    do until i>n*k;
        call xy(l,bsvec[i,.]');
        i=i+1;
    endo;
    */
    retp(bsvec');
endp;