/* ============================================ 
    Stock & Watson, " Vector Autoregressions "(2001)

Recursive VAR(4) : Backward Looking Monetary policy
Sample Period :  1960:1~ 2001:3, monthly.

Two methods to estimate invB;
Method 1: Cholessky decomposition
Method 2: Maximum likelihood Estimation

This program is based on the Method 1 !!.

Note that
e1 : Inflation shock
e2 : Unemployment rate shock
e3 : Monetary policy Shock
============================================= */

New;
library pgraph;

cls;
t0=167;      @ Population period, 1960:1 ~ 2001:3 @
k=3;           @ Dimension of VAR @
p=4;           @ # of Lag @
j=24;          @ Length of horizon for impulse response @

/*======== loading data  ==============*/

load data[t0,3] =c:\data\watson.txt;
data=data[1:t0,.]; @ sample @

inf=data[.,1];      @ inflation rate @
un=data[.,2];      @ Un @
ffr=data[.,3];       @ ffr @

y=inf~un~ffr;
y0=y[p+1:t0,.]; @ t0-p by k, Y(t) @

y1=y[4:t0-1,.]; @ Y(t-1) @
y2=y[3:t0-2,.]; @ Y(t-2) @
y3=y[2:t0-3,.]; @ Y(t-3) @
y4=y[1:t0-4,.]; @ Y(t-4) @

y_lag=y1~y2~y3~y4;

@ Demean @
y0=y0-meanc(y0)';
y_lag=y_lag-meanc(y_lag)';

/*===== Estimation of Phi and Omega_hat ====*/
phi_hat=inv(y_lag'y_lag)*y_lag'y0; @ p*k by k @
 
F=phi_hat'|(eye((p-1)*k)~zeros(k*(p-1),k));  @ p*k by p*k @

T=rows(y0); @ rows of Y @

@ Omega_hat @
u_hat=y0-y_lag*phi_hat; @ t-p by k @
Omeg=u_hat'*u_hat/T;  @ variance- covariance matrix,  k by k @


/* ===== Estimation of inv(B) using Cholesky decomposition========*/
@ calculating the Inverse matrix of B by Cholesky Decomposition @ 
inv_B=chol(Omeg)' ; @ Lower triangular matrix @

@======= Impulse Response ===============@
   theta_i=zeros(j+1,k); @ Impluse response of inflation to each shock @
   theta_u=zeros(j+1,k); @ Impluse response of Un to each shock @
   theta_f=zeros(j+1,k); @ Impluse response of FFR to each shock @

  @ theta_f[.,1] = to e1t,theta_f[.,2] = to e2t,theta_f[.,3] = to e3t @

   FF = eye(p*k);
   itr = 1; 
   do until itr>(j+1);

              psi_j = FF[1:k,1:k]; @ k by k @
              theta = psi_j*inv_B; @ k by k @

              theta_i[itr,.]=theta[1,.]; @ 1 by k, Inf @
              theta_u[itr,.]=theta[2,.]; @ 1 by k, Un @
              theta_f[itr,.]=theta[3,.]; @ 1 by k, FFR @
       
              FF = FF*F;
              itr = itr+1;
   endo;



/* ======PLOT IMPULSE RESPONSES ============ */

 graphset;
   begwind;
  window(3,3,0);

    setwind(1);
    title("Inflation shock to inflation");
    xy( seqa(0,1,j+1), theta_i[.,1]~zeros(j+1,1));

    setwind(2);
    title("Inflation shock to Un");
    xy(seqa(0,1,j+1),theta_u[.,1]~zeros(j+1,1));

    setwind(3);
    title("Inflation shock to FFR");
    xy(seqa(0,1,j+1), theta_f[.,1]~zeros(j+1,1));

    setwind(4);
    title("Un shock to Inflation");
    xy(seqa(0,1,j+1),theta_i[.,2]~zeros(j+1,1));

   setwind(5);
   title("Un shock to Un");
   xy(seqa(0,1,j+1),theta_u[.,2]~zeros(j+1,1));

    setwind(6);
    title("Un shock to FFR");
    xy(seqa(0,1,j+1),theta_f[.,2]~zeros(j+1,1));

    setwind(7);
    title("FFR shock to Inflation");  
    xy(seqa(0,1,j+1),theta_i[.,3]~zeros(j+1,1));

    setwind(8);
    title("FFR shock to Un");  
    xy(seqa(0,1,j+1),theta_u[.,3]~zeros(j+1,1));

    setwind(9);
    title("FFR shock to FFR");  
    xy(seqa(0,1,j+1),theta_f[.,3]~zeros(j+1,1));

   endwind;

end;