% Matlab File: plotimpflexiblerates.m

%Chapter 7

%MIUF model with flexible price and flexible exchange rates
%Endowment economy, two goods

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       First Version:June 2007
%       This Version(revised by AR):July 22 2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%This program is part of the book "Open Economy Macroeconomics in
%Developing Countries" (forthcoming MIT Press) by Carlos A. Vegh. Use of
%this program or modifications thereof should be acknowledged by citing
%this manuscript.  This program is a modification of the original
%King, Plosser, and Rebelo's (1988) programs. I am grateful to Guillermo
%Vuletin for his help in writing this program.  

close all

global shock


% This program is to be used after implmbs.m to plot impulse response functions

% Order of variables is: 
% States, exogenous, co-states, controls, and flows  (IR matrix)

%      | state    | exogenous     |  costates   |    endogenous                   |   flow                           
% Date    k         mu   yT   yN      lambda  m    cT  cN   epsilon  i    e     c      TB   
%  1      2          3   4    5         6     7     8   9    10      11   12    13     14

m=IR(:,7);   % this is just to check that the new steady state value of m is different from the initial one. In this particular case it is lower. 

%  Program was modified to include the zero value in period zero
%  in the plot

%clg
%time = 1:nir;

time = 0:nir;
[dim1 dim2]=size(IR);
PIR = [zeros(1,dim2) 
       IR];
   
       
%*****************************************************************************************
%We now add a dynamic flow variable (piN)

piN = mu; %this the steady-state value

[r c] = size(PIR);
PIRL = zeros(r,c);

%Here we define the lagged impulse response

PIRL(2:TT,:)=PIR(1:TT-1,:); %these are elements 2 and onwards

%note that the first element of PIRL is zero


piNhat = (epsilon/(1+epsilon))*((1+piN)/piN)*PIR(:,10) + ((1+piN)/piN)*PIRL(:,12)-((1+piN)/piN)*PIR(:,12);

piNhat(16,1) = piNhat (15,1); %we set last element to zero (i.e., equal to the next to last one); otherwise we get an incorrect last number
%**************************************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                        THIS IS THE NEW PART ADDED BY AR

%% Computing the path for E (nominal exchange rate)
%%


%% 1)   Computing the path for epsilon (NOT epsilon HAT, just epsilon)

eps=IR(:,10);  % this is epsilon hat (computed by the code) . I use this to compute piNhat (see below). (The code doesn't give me the 1st epsilon hat, but it is obviously zero)
e=IR(:,12);    % real exchange rate
E(1)=epsilon;
for i=1:size(IR,1)
E(i+1)=epsilon*(1+(eps(i)/100));     % (recall that epsilon is the steady state value for epsilon)
end


E1=(E(1)-epsilon)/epsilon;   % first element for the sequence of epsilon HAT. Notice that to get this I am using the first element of the sequence for epsilon.

epsh=[E1;eps];               % this is the whole sequence for epsilon HAT. 


%% 2)   Computing the path of E (NOT E HAT, but E)

% code notation.
% E is the rate of depreciation and EE is the actual level of the nominal exchange rate

%EE(1)=1;                
EE(1)=0;

%---------------- Second E (that is E_1) _ ( nominal exchange rate level)------------------------------
M0h=IR(1,3)/100;    %  this is mu
m0h=IR(1,7)/100;
E0h=M0h-m0h;    %  this is Ehat

EE(2)=log(E0h+1);   % this is E_1 (assuming that E_0=1)
%EE=[Ee EE(1:end-2)];
%------------------------------------------------------------------------------------------------------

Eaux=[E(2:end) E(end)];     % this is exogenous, we are choosing it.
for i=2:length(Eaux);       % This loops starts with the chosen value of the nominal exchange rate and iterates to get the whole path.
%EE(i+1)=(1+Eaux(i))*EE(i);
EE(i+1)=log(1+Eaux(i))+EE(i);
end


EE=[EE(1:end-1)];



%% 3)   Nominal exchange rate without stabilization

t=[0:length(EE)-1];
for i=1:length(t)
%Et(i)=(1+epsilon)^t(i);      % Notice that as above, here I am also starting at a nominal exchange rate level of 1 
Et(i)=t(i)*log(1+epsilon);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Notice that to make a 6 panel figure, we need to have the command
%'subplot(32x) where 32 refers to 3 rows of 2 panels each and x refers to
%the panel number

% Date  k   mu   yT   yN      lambda  m    cT  cN   epsilon i    e     c    TB   
%  1    2    3   4    5         6     7    8   9    10      11   12    13   14


figure(1)

if shock==1
subplot(421), plot(time,PIR(:,3),time,zeros(size(time)),'.')
title('money growth rate')
elseif shock==2
subplot(421), plot(time,PIR(:,4),time,zeros(size(time)),'.')
title('yT')
elseif shock==3
subplot(421), plot(time,PIR(:,5),time,zeros(size(time)),'.')
title('yN')
end


subplot(422), plot(time,PIR(:,11),time,zeros(size(time)),'.') 
title('Nominal intereset rate')


subplot(423), plot(time, PIR(:,8),time,zeros(size(time)),'.')
title('Consumption of tradables')


subplot(424), plot(time,PIR(:,7),time,zeros(size(time)),'.')
title('real money balances')


subplot(425), plot(time,PIR(:,14),time,zeros(size(time)),'.') 
title('Trade balance')


subplot(426), plot(time, PIR(:,12),time,zeros(size(time)),'.')
title('Real exhange rate')

%subplot(427),plot(time,piNhat_a,time,zeros(size(time)),'.')
subplot(427),plot(time,piNhat,time,zeros(size(time)),'.')
title('Inflation rate of Non-tradables')
xlabel('time')


subplot(428), 
plot(time,Et)
hold on
plot(time,EE,'ro-')
legend('E without stabilization','E with stabilization')
% legend('E without stabilization','E with stabilization','location','northwest')
xlabel('time')
title('Nominal Exchange rate')
xlabel('time')

% Super Title

[ax,h3]=suplabel('Figure 9.  Temporary reduction in money growth rate' ,'t');
set(h3,'FontSize',10)

print -dpsc2 FigCh7.ps

figure
plot(time,E)
title('Epsilon(NOT epsilon HAT)')
xlabel('time')

figure
subplot(211)
plot(time,EE)
title('level of E with stabilization')
xlabel('time')


subplot(212)
plot(time,Et)
title('level of E without stabilization')
xlabel('time')

figure
plot(time,Et)
hold on
plot(time,EE,'ro-')
legend('E without stabilization','E with stabilization','location','northwest')
xlabel('time')

for i=1:length(EE)
dif(i)=(EE(i)-Et(i));
end
%dif=log(dif);
%figure
%plot(time,dif,'k')
%title('difference between E without stabilization and E without')
%xlabel('time')