% Program to simulate, estimate and identify a structural VAR(1)
clear all
clc
%--------------------------------------------------------------------------
% Define true model
%--------------------------------------------------------------------------
A0=[1,0;
    -.2,1];
A1=[0.9,-.2;
    .2,.75;];
C=inv(A0);
FI=inv(A0)*A1;
OMEGA=inv(A0)*(inv(A0))';
%--------------------------------------------------------------------------
% Simulate data
%--------------------------------------------------------------------------
T=100;
Y=zeros(2,T);
Y(:,1)=[2;-3;];
U=zeros(2,T);
for tt=2:T;
   U(:,tt)=randn(2,1);
   Y(:,tt)=FI*Y(:,tt-1) + C*U(:,tt);
end
%--------------------------------------------------------------------------
% Plot time series
%--------------------------------------------------------------------------
figure(1)
subplot(2,1,1);
plot(Y(1,:),'k','linewidth',2);
title('Y1','fontsize',18);
subplot(2,1,2);
plot(Y(2,:),'k','linewidth',2);
title('Y2','fontsize',18);

%--------------------------------------------------------------------------
% Estimate reduced form
%--------------------------------------------------------------------------
y=Y(:,2:end);
y1=Y(:,1:end-1);

FIhat=y*y1'*inv(y1*y1');
OMEGAhat=inv(T-1)*(y-FIhat*y1)*(y-FIhat*y1)';

%--------------------------------------------------------------------------
% Find estimated structural form using a12=0
%--------------------------------------------------------------------------
Chat=chol(OMEGAhat)';
A0hat=inv(Chat);
A1hat=inv(Chat)*FIhat;

%--------------------------------------------------------------------------
% Compute impulse response to structural shocks
%--------------------------------------------------------------------------
S=20;
IRF1=zeros(2,S);
IRF2=zeros(2,S);
for ss=1:S;
    IRF1(:,ss)=FIhat^(ss-1)*Chat(:,1);
    IRF2(:,ss)=FIhat^(ss-1)*Chat(:,2);
end
%--------------------------------------------------------------------------
% Plot impulse response to structural shocks
%--------------------------------------------------------------------------
figure(2)

subplot(2,2,1);hold on;
plot(IRF1(1,:),'k','linewidth',2);
title('Y1 to U1','fontsize',18);

subplot(2,2,2);hold on;
plot(IRF1(2,:),'k','linewidth',2);
title('Y2 to U1','fontsize',18);

subplot(2,2,3);hold on;
plot(IRF2(1,:),'k','linewidth',2);
title('Y1 to U2','fontsize',18);

subplot(2,2,4);hold on;
plot(IRF2(2,:),'k','linewidth',2);
title('Y2 to U2','fontsize',18);

%--------------------------------------------------------------------------
% Compute unconditional variance
%--------------------------------------------------------------------------
diff=1;
tol=1e-6;
SIGYYst=zeros(2,2);

while diff > tol;
    SIGYY=FI*SIGYYst*FI'+C*C';
    diff=max(max(abs(SIGYY-SIGYYst)));
    SIGYYst=SIGYY;
end

%--------------------------------------------------------------------------
% Compute variance due to individual shocks
%--------------------------------------------------------------------------
 
diff=1;
SIGYY1st=zeros(2,2);

while diff > tol;
    SIGYY1=FI*SIGYY1st*FI'+C(:,1)*C(:,1)';
    diff=max(max(abs(SIGYY1-SIGYY1st)));
    SIGYY1st=SIGYY1;
end

diff=1;
SIGYY2st=zeros(2,2);

while diff > tol;
    SIGYY2=FI*SIGYY2st*FI'+C(:,2)*C(:,2)';
    diff=max(max(abs(SIGYY2-SIGYY2st)));
    SIGYY2st=SIGYY2;
end
    
frac1=diag(SIGYY1)./diag(SIGYY); 
frac2=diag(SIGYY2)./diag(SIGYY);