/*
** Lesson 13.1: Klein's Model I
** Simultaneous Equation System of Klein's Model I
*/
use gpe2;
output file = gpe\output13.1 reset;
load data[23,10] = gpe\klein.txt;

a=data[2:23,1]-1931;      @ time trend: 1931 = 0 @
c=data[2:23,2];           @ consumption @
p=data[2:23,3];           @ profit income @
w1=data[2:23,4];          @ private wage income @
i=data[2:23,5];           @ investment @
k1=data[2:23,6];          @ lagged capital stock @
x=data[2:23,7];           @ private total income @
w2=data[2:23,8];          @ public wage income @
g=data[2:23,9];           @ government spending @
t=data[2:23,10];          @ tax @

k=k1[2:22]|209.4;         @ capital stock @
w=w1+w2;                  @ total wage income @

yvar=c~i~w1~x~p~k~w;
xvar=lag1(x~p~k)~w2~a~g~t;

call reset;

_names={"c","i","w1","x","p","k","w",
        "x-1","p-1","k-1","w2","a","g","t"};
_vcov=1;
        @ C  I W1  X  P  K  W XL PL KL W2  A  G  T  1 @
_eq =   {-1  0  0  0  1  0  1  0  1  0  0  0  0  0,
          0 -1  0  0  1  0  0  0  1  1  0  0  0  0,
          0  0 -1  1  0  0  0  1  0  0  0  1  0  0};
_id =   { 1  1  0 -1  0  0  0  0  0  0  0  0  1  0,
          0  0 -1  1 -1  0  0  0  0  0  0  0  0 -1,
          0  1  0  0  0 -1  0  0  0  1  0  0  0  0,
          0  0  1  0  0  0 -1  0  0  0  1  0  0  0};

_begin=2;

_method=0;        @ OLS estimation @
call estimate(yvar,xvar);
_method=1;        @ LIML estimation @
call estimate(yvar,xvar);
_method=2;        @ 2SLS estimation @
call estimate(yvar,xvar);
/*
_iter=100;
_method=3;        @ 3SLS estimation (iterative) @
call estimate(yvar,xvar);
_method=4;        @ FIML estimation @
call estimate(yvar,xvar);
*/
print __pi;       
@ 8x7 matrix (3 equations + 4 identities) @
@ (3 lagged endogenous + 5 exogenous, incl. constant) @
p1=zeros(3,7)|__pi[1:3,1:7]|zeros(1,7); @ lag1: C I W1 X P K W2 @
p2=__pi[6:7,1:7]; @ 6:G; 7:T @

print impulse(p1,p2,10);

print p2*inv(eye(7)-p1); @ equilibrium or lone-run effects @

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;