clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% WARNING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Activate line 11 and 23 in final, before running the program %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

calibration % call calibration.m in order to solve the equations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Impulse Responses %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

shock = [1,0,0,0,0,0,0]'; % [Ax, An, r, f, g, Px, Ln], if Ax hikes 1%, [1,0,0,0,0,0,0]'
period = 25;

seris = [];
added = zeros(28,1);
seris = [seris added];
added = G1*added+impact*shock;
seris = [seris added];

for j= 2:period
   added = G1*added;
   seris = [seris added];
end

   uhat = seris(1,:);
   chat = seris(2,:);
   cxhat = seris(3,:);
   cmhat = seris(4,:);
   cnhat = seris(5,:);
   yxhat = seris(6,:);
   ynhat = seris(7,:);
   ixhat = seris(8,:);
   inhat = seris(9,:);
   lhat = seris(10,:);
   hxhat = seris(11,:);
   hnhat = seris(12,:);
   shat = seris(13,:);
   kxhat = seris(14,:);
   knhat = seris(15,:);
   bhat = seris(16,:);
   pnhat = seris(17,:);
   phat = seris(18,:);   
   lamdaonehat = seris(19,:);
   lamdathreehat = seris(20,:);
   lamdafourhat = seris(21,:);
   axhat = seris(22,:);
   anhat = seris(23,:);
   rhat = seris(24,:);
   fhat = seris(25,:);
   gnhat = seris(26,:);
   pxhat = seris(27,:);
   lnhat = seris(28,:);

   %aggregate variable in terms of import good prices: c,y,i,tb,ca
   
   yhat=(px*yx/y)*(pxhat+yxhat)+(pn*yn/y)*(pnhat+ynhat);
   cchat=(px*cx/cc)*(pxhat+cxhat)+(cm/cc)*cmhat+(pn*cn/cc)*(pnhat+cnhat); %p*c, aggregate nominal consumption
   ihat=(ix/i)*ixhat+(in/i)*inhat;
   nxhat=yhat-cc/y*cchat-g/y*(pnhat+gnhat)-(i/y)*ihat-(s/y)*shat;  %percentage from y
   cahat=nxhat(2:period+1)+px*yx/y*fhat(2:period+1)+bb/y*rhat(1:period)+px*yx*r/y*bhat(1:period);  %percentage from y
   cahat=[0,cahat];
   hhat=(hx/h)*hxhat+(hn/h)*hnhat;
   
   %aggregate price variables
   % RER and CPI both are phat.
    
   idefhat=(1+r)*(phat(2:period+1)-phat(1:period));
   idefhat=[0,idefhat];
   
   %sectoral varaibles
   
   pxyxhat=pxhat+yxhat;  %px*yx
   pxcxhat=pxhat+cxhat;  %px*cx
   pxixhat=ixhat;
   
   pnynhat=pnhat+ynhat; %pn*yn
   pncnhat=pnhat+cnhat; %pn*cn
   pninhat=inhat;

     
%Draw impulse responses
pp=[0:period];
zero       = zeros(size(pp));

figure;

subplot(3,2,1);
plot(pp,yhat,'-k',pp,yxhat,'--k', pp,ynhat,':k');
title('Output');
ylabel('Percent');
xlabel('Period');
line(pp,zero,'linestyle','-','LineWidth',0.5)
legend('Y', 'Y_x', 'Y_n');

subplot(3,2,2);
plot(pp,nxhat,'-k',pp,cahat,'--k');
title('External Balances');
ylabel('Percent');
xlabel('Period');
line(pp,zero,'linestyle','-','LineWidth',0.5)
legend('Nx', 'CA');

subplot(3,2,3);
plot(pp,cchat,'-k',pp,cxhat,'--k', pp,cnhat,':k');
title('Consumption');
ylabel('Percent');
xlabel('Period');
line(pp,zero,'linestyle','-','LineWidth',0.5)
legend('CC', 'C_x', 'C_n');

subplot(3,2,4);
plot(pp,phat,'-k',pp,pnhat,'--k');
title('Price');
ylabel('Percent');
xlabel('Period');
line(pp,zero,'linestyle','-','LineWidth',0.5)
legend('RER', 'p_n');

subplot(3,2,5);
plot(pp,ihat,'-k',pp,ixhat,'--k', pp,inhat,':k');
title('Investment');
ylabel('Percent');
xlabel('Period');
line(pp,zero,'linestyle','-','LineWidth',0.5)
legend('i', 'i_x', 'i_n');

subplot(3,2,6);
plot(pp,hhat,'-k',pp,hxhat,'--k', pp,hnhat,':k');
title('Labor');
ylabel('Percent');
xlabel('Period');
line(pp,zero,'linestyle','-','LineWidth',0.5)
legend('h', 'h_x', 'h_n');