Impulse response function with sign restrictions

Hello everyone!

I'm doing a IRF with sign restriction with three inputs gdp, cpi and discount rate. The first two are ok i have taken the growth rate and do the impulse response function with restriction and cumsum. The problem i have is with the third one. I know I have to take the discount rate in level and not no use the cumsum.

To be clear, the VAR for each region will be estimated with 3 series: i) output growth; ii) inflation; iii) discount rate in level and not in growth rate. So as the discount rate is in level, I should not take cumsum of the irf for the discount rate because the cumsum helps to get the irf in level. I'm really stuck in this part so if you can help me to do it it would be great. I thought that by using sir instead of sirc it was right but apparently i'm wrong.

clear all; close all; clc;

format long;

randn('seed',123456789); % Set the seed in order to replicate the same results whenever you run the code

[A,B]= xlsread('DataMemoIr','Lima'); % load data.

gdp1=A(2:end,2)

gdp2=A(1:end-1,2)

gdp =((gdp1./gdp2)-1)*100 % take the growth rates of the gdp

cpi = A(2:end,3)%ipc growth rate

ir =A(2:end,4) % discount rate

y_raw=([cpi gdp ir])

data=y_raw;

T=size(data,1);

%***********************************************************

% Preliminary transformation of the data (constant, linear or quadratic trend) of the data if needed

%***********************************************************

const=ones(size(data,1),1); % vector of ones for constant term

% you can also add trend to the model

trend=1:size(data,1); % linear trend term

qtrend=(1:size(data,1)).^2; % quadratic trend

% XXX=[const]; % put constant and trend terms together if any

% XXX=[const trend' qtrend']; % put constant and trend terms together if any

XXX=[const trend']; % put constant and trend terms together if any

% residuals: regressing data on a constant and saving residuals

% this is de-trending data for constant and trend terms such that

% coefficient estimated later for VAR are directly the AR coefficients

Bols_det=inv(XXX'*XXX)*XXX'*data; %OLS estimate on deterministic terms

res1=data-XXX*Bols_det; % de-trending data for constant and trend

%------------------------------

% define some useful quantities

%-------------------------------

p = 4; % number of lags; We can change here by selecting the optimal lag found for your data

nPeriods = 16; % number of impulse periods

c = 0 ; % use no constant in the VAR ounce we have already removed these from the data

n=size(data,2); % number of variables in the VAR

ho = 12; % number of months over which sign restrictions are imposed

Adraws=500; % maximum number of accpeted draws for the restrictions

% estimation of the var system

[Y,X] = VAR_str(res1,c,p); % creates the variables X and the Y

Bet=inv(X'*X)*X'*Y; % beta_ols

b1=Bet';

% residuals and their covariance matrix

res=Y-X*Bet;

Sigma=cov(res); % sigma_ols

%------------------------------

% COMPUTE FIRST CholeskY IRF

%-------------------------------

wimpu=[];

cimpu=[];

BB=companion2(Bet,p,n,c); % construct companion form of the VAR

I = eye(n*p);

BETA = [b1;I(1:end-n,1:end)]; % construct companion form of the VAR

for j=1:nPeriods

if j==1

wimpu(:,:,j)=eye(size(BB,1),size(BB,2)); % (non-orthogonal) Wold impulses

elseif j>1

wimpu(:,:,j)=BB^(j-1); %(non-orthogonal) Wold impulses

end

cimpu(:,:,j)=wimpu(1:n,1:n,j)*chol(Sigma)'; % (orthogonal) Choleski IRF

end

%---------------------------------------------

%Define some useful parameters to be used below

%---------------------------------------------

dr=1; % to control for total number of draws

gg=0; % to control for the accepted draws

irf=[]; % initialize the IRF

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

%Identify the structural shocks with Sign restrictions

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

while gg<=Adraws

%---------------------------------------------------------------

%1. draw a random matrix and take appropriate tra,formations

%-------------------------------------------------------------

a=normrnd(0,1,n,n);

[q r]=qr(a); % Compute qr decomposition of random matrix a

for ii=1:n;

if r(ii,ii)<0

q(:,ii)=-q(:,ii); % flip the sign

end

end

%-------------------------------------------

% 2. Compute candidate IRF implied by draw in q

%-------------------------------------------

for j=1:n %shock

for i = 1:nPeriods

S(:,j,i) = cimpu(:,:,i)*q(:,j);

end

end

%-----------------------------------

% 3. compute cumulated candidate IRF

%------------------------------------

A=zeros(size(S,1),size(S,2),nPeriods);

for k=1:nPeriods % horizon

if k==1 % impact period

A(:,:,1)=S(:,:,1);

else

for c=1:size(S,2),

for r=1:size(S,1),

A(c,r,k)=A(c,r,k-1)+S(c,r,k);

end;

end;

end;

end;

%--------------------------------------------------------------------------------------------------------------

% 4. check if the draw satisfies the sign restrictions implied by ADAS model are satisfied for responses to shocks at horizons 1 to ho

%---------------------------------------------------------------------------------------------------------------

z1=[squeeze(A(1,1,1:ho))'>=0 squeeze(A(2,1,1:ho))'<=0 squeeze(A(1,2,1:ho))'>=0 squeeze(A(2,2,1:ho))'>=0 ];

zz1=sum(z1);

if zz1==size(z1,2) % check if the restrictions are satisfied

disp('I found one accepeted draw')

gg=gg+1;

for j=1:n % for the shocks

for i = 1:nPeriods

irf(:,j,i,gg) = cimpu(:,:,i)*q(:,j); % Now compute the identified IRF for difference series

end

end

end

dr=dr+1;

end

disp('number of draws made and accepted draws')

[dr-1, gg]

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

% Now compute valid cumulated resp i.e. IRF for series in levels

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Ac=zeros(n,n,nPeriods,gg);

for jjj=1:gg % do this for each accepeted draw

for k=1:nPeriods,

if k==1,

Ac(:,:,1,jjj)=irf(:,:,1,jjj);

else

for c=1:size(S,2),

for r=1:size(S,1),

Ac(c,r,k,jjj)=Ac(c,r,k-1,jjj)+irf(c,r,k,jjj);

end;

end;

end;

end;

end;

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

% Now plot IFRs

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

sir=sort(irf,4); % irf in difference

sirc=sort(Ac,4);% level

% estract 68 intervals for the responses using identification uncertainty

mec11=squeeze(sirc(1,1,:,fix(gg*0.50)));

upc11=squeeze(sirc(1,1,:,fix(gg*0.84)));

loc11=squeeze(sirc(1,1,:,fix(gg*0.16)));

mec21=squeeze(sirc(2,1,:,fix(gg*0.50)));

upc21=squeeze(sirc(2,1,:,fix(gg*0.84)));

loc21=squeeze(sirc(2,1,:,fix(gg*0.16)));

mec31=squeeze(sir(3,1,:,fix(gg*0.50)));

upc31=squeeze(sir(3,1,:,fix(gg*0.84)));

loc31=squeeze(sir(3,1,:,fix(gg*0.16)));

mec12=squeeze(sirc(1,2,:,fix(gg*0.50)));

upc12=squeeze(sirc(1,2,:,fix(gg*0.84)));

loc12=squeeze(sirc(1,2,:,fix(gg*0.16)));

mec22=squeeze(sirc(2,2,:,fix(gg*0.50)));

upc22=squeeze(sirc(2,2,:,fix(gg*0.84)));

loc22=squeeze(sirc(2,2,:,fix(gg*0.16)));

mec32=squeeze(sir(3,2,:,fix(gg*0.50)));

upc32=squeeze(sir(3,2,:,fix(gg*0.84)));

loc32=squeeze(sir(3,2,:,fix(gg*0.16)));

mer=[mec11,mec12,mec21,mec22,mec31,mec32];

figure;

subplot(2,2,1),plot([mec11 upc11 loc11]),axis tight;

legend('Output')

grid on;

title('Supply shock','FontWeight','bold','FontSize',10)

subplot(2,2,2),plot([mec21 upc21 loc21]),axis tight;

legend('Price')

grid on;

title('Supply shock','FontWeight','bold','FontSize',10)

subplot(2,2,3),plot([mec12 upc12 loc12]),axis tight;

legend('Output')

grid on;

title('Demand shock','FontWeight','bold','FontSize',10)

subplot(2,2,4),plot([mec22 upc22 loc22]),axis tight;

legend('Price')

grid on;

title('Demand shock',...

'FontWeight','bold','FontSize',10)

xrom=1:nPeriods; %this is used for Xasis

yrom=zeros(1,nPeriods);

figure;

plot(1:nPeriods,loc11,'--b','LineWidth',2),axis tight;hold on;

plot(1:nPeriods,upc11,'--b','LineWidth',2),axis tight;hold on;

plot(1:nPeriods,mec11,'--r','LineWidth',2);hold off

legend('Output'), grid on;line(xrom,yrom);

title('Supply shock','FontWeight','bold','FontSize',10)

set(gca,'XLim',[1 nPeriods],'FontSize',15);xlabel('\bf monthly','FontSize',20)

figure;

plot(1:nPeriods,loc12,'--b','LineWidth',2),axis tight;hold on;

plot(1:nPeriods,upc12,'--b','LineWidth',2),axis tight;hold on;

plot(1:nPeriods,mec12,'--r','LineWidth',2);hold off

legend('Output'), grid on;line(xrom,yrom);

title('Demand shock','FontWeight','bold','FontSize',10)

set(gca,'XLim',[1 nPeriods],'FontSize',15);xlabel('\bf monthly','FontSize',20)

figure;

plot(1:nPeriods,loc21,'--b','LineWidth',2),hold on;

plot(1:nPeriods,upc21,'--b','LineWidth',2),hold on;

plot(1:nPeriods,mec21,'--r','LineWidth',2);hold off

legend('Prices'), grid on;line(xrom,yrom);

title('Supply shock','FontWeight','bold','FontSize',10)

set(gca,'XLim',[1 nPeriods],'FontSize',15);xlabel('\bf monthly','FontSize',20)

figure;

plot(1:nPeriods,loc22,'--b','LineWidth',2),hold on;

plot(1:nPeriods,upc22,'--b','LineWidth',2),hold on;

plot(1:nPeriods,mec22,'--r','LineWidth',2);hold off

legend('Prices'), grid on;line(xrom,yrom);

title('Demand shock','FontWeight','bold','FontSize',10)

set(gca,'XLim',[1 nPeriods],'FontSize',15);xlabel('\bf monthly','FontSize',20)

figure;

plot(1:nPeriods,loc31,'--b','LineWidth',2),hold on;

plot(1:nPeriods,upc31,'--b','LineWidth',2),hold on;

plot(1:nPeriods,mec31,'--r','LineWidth',2);hold off

legend('Discount rate'), grid on;line(xrom,yrom);

title('Supply shock','FontWeight','bold','FontSize',10)

set(gca,'XLim',[1 nPeriods],'FontSize',15);xlabel('\bf monthly','FontSize',20)

figure;

plot(1:nPeriods,loc32,'--b','LineWidth',2),hold on;

plot(1:nPeriods,upc32,'--b','LineWidth',2),hold on;

plot(1:nPeriods,mec32,'--r','LineWidth',2);hold off

legend('Discount rate'), grid on;line(xrom,yrom);

title('Demand shock','FontWeight','bold','FontSize',10)

set(gca,'XLim',[1 nPeriods],'FontSize',15);xlabel('\bf monthly','FontSize',20)