%
% Fixed exchange rate & Flexible prices
% Country A pegs
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                 Parameters of the algorithm                 %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ncont=6;                  % # of controls
nbend=0;                  % # state, backward, endogenous
nshoc=4;                  % # shocks
nback=nbend+nshoc;        % # state backward
nforw=2;                  % # state forward
nstat=nback+nforw;        % # total state
Mcc=zeros(ncont,ncont);   % Matrixe Control-Control
Mcs=zeros(ncont,nstat);   % Matrice Control-State
Mss0=zeros(nstat,nstat);  % Matrice State-State t+1
Mss1=zeros(nstat,nstat);  % Matrice State-State t
Msc0=zeros(nstat,ncont);  % Matrice State- Control t+1
Msc1=zeros(nstat,ncont);  % Matrice State-Contol t
Mse=zeros(nstat,nshoc);   % Matrice State-Shocks
tol=1e-6;
nrep=20;

%%outfile='fixA.res';			% Output file name

par;							% Call the parameters

%
% Mcc
%
Mcc(1,1)=a11;
Mcc(1,5)=a21;
Mcc(2,2)=a12;
Mcc(2,5)=-a22;
Mcc(3,3)=b21;
Mcc(3,6)=-1;
Mcc(4,4)=b22;
Mcc(5,5)=1;
Mcc(6,3)=1;
Mcc(6,4)=-1;


%
% Mcs
%
Mcs(1,1)=-a31;
Mcs(2,2)=-a32;
Mcs(3,3)=-1;
Mcs(3,5)=-1;
Mcs(4,4)=-1;
Mcs(4,6)=-1;
Mcs(5,5)=-1;
Mcs(5,6)=1;



%
% Mss0
%
Mss0=eye(6);
%
% Mss1
%
Mss1(1,1)=-rg11;
Mss1(1,2)=-rg12;
Mss1(2,1)=-rg12;
Mss1(2,2)=-rg22;
Mss1(3,3)=-re11;
Mss1(3,4)=-re12;
Mss1(4,3)=-re12;
Mss1(4,4)=-re22;
Mss1(5,5)=-1;
Mss1(6,6)=-1;

%Mss1(5:7,5:7)=-eye(3);
%
% Msc0
%
% Empty !
%
% Msc1
%
Msc1(5,1)=-1;
Msc1(5,3)=1;
Msc1(6,2)=-1;
Msc1(6,4)=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%            Solution							                 %
%                                                             %
%           -1                                                %
% X(t) = Mcc  Mcs S(t)                                        %
%                                                             %
%                      -1             -1                      %
% S(t+1)=(Mss0-Msc0 Mcc  Mcs)(Msc1 Mcc  Mcs-Mss1)S(t)         %
%                       -1                                    %
%        +(Mss0-Msc0 Mcc  Mcs)e(t+1)                          %
%                                                             %
% S(t+1) = W S(t) + R e(t+1)                                  %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M0=inv(Mss0-Msc0*inv(Mcc)*Mcs);
M1=(Mss1-Msc1*inv(Mcc)*Mcs);
W=-M0*M1;
%
% Calculation of eigenvalues : MU -> Eigenvalues, P -> Eigenvectors
%
[P,MU] = eig(W);

% Sorting the eigenvalues

AMU=abs(MU);
%

%Rearrange system according to the absolute value of eignevalues

[MU,k] = sort(diag(AMU));
P=P(:,k);
Q=inv(P);
%
% Solution
%

L=-inv(Q(nback+1:nstat,nback+1:nstat))*Q(nback+1:nstat,1:nback);
M=[eye(nback);-inv(Q(nback+1:nstat,nback+1:nstat))*Q(nback+1:nstat,1:nback)];
MSS=W(1:nback,:)*M;
M2=M0*Mse;
MSE=M2(1:nback,1:nshoc);
PI=inv(Mcc)*Mcs*M;

PIP=L*MSS;

%Fixed prices +fixed exchange rates

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                 Parameters of the algorithm                 %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ncont=10;                 % # of controls (8 controls + 2 prices)
nbend=2;                  % # state, backward, endogenous
nshoc=4;                  % # shocks
nback=nbend+nshoc;        % # state backward
nforw=0;                  % # state forward
nstat=nback+nforw;        % # total state
Mcc=zeros(ncont,ncont);   % Matrixe Control-Control
Mcs=zeros(ncont,nstat);   % Matrice Control-State
Mss0=zeros(nstat,nstat);  % Matrice State-State t+1
Mss1=zeros(nstat,nstat);  % Matrice State-State t
Msc0=zeros(nstat,ncont);  % Matrice State- Control t+1
Msc1=zeros(nstat,ncont);  % Matrice State-Contol t
Mse=zeros(nstat,nshoc);   % Matrice State-Shocks
tol=1e-6;
nrep=20;


%
% Mcc
%
Mcc(1,1)=a11;
Mcc(1,5)=a21;
Mcc(1,6)=-1;
Mcc(1,7)=a41;

Mcc(2,2)=a12;
Mcc(2,5)=-a22;
Mcc(2,6)=a42;
Mcc(2,7)=-1;


Mcc(3,3)=b21;
Mcc(3,6)=b11;
Mcc(3,8)=-1;


Mcc(4,4)=b22;
Mcc(4,7)=b12;

Mcc(5,5)=1;

Mcc(6,1)=-1;
Mcc(6,3)=1;
Mcc(7,2)=-1;
Mcc(7,4)=1;

Mcc(8,3)=1;
Mcc(8,4)=-1;


Mcc(9,9)=1;
Mcc(10,10)=1;

%
% Mcs
%
Mcs(1,3)=-a31;
Mcs(2,4)=-a32;

Mcs(3,1)=-1;
Mcs(3,5)=-1;

Mcs(4,2)=-1;
Mcs(4,6)=-1;

Mcs(5,1)=-1;
Mcs(5,2)=1;

Mcs(6:7,1:2)=-eye(2);
Mcs(6:7,3:6)=PIP;

Mcs(9,1)=1;
Mcs(10,2)=1;


%
% Mss0
%
Mss0=eye(6);
%
% Mss1
%
Mss1(1:2,3:6)=-PIP;
Mss1(3,3)=-rg11;
Mss1(3,4)=-rg12;
Mss1(4,3)=-rg12;
Mss1(4,4)=-rg22;
Mss1(5,5)=-re11;
Mss1(5,6)=-re12;
Mss1(6,5)=-re12;
Mss1(6,6)=-re22;

%
% Msc0
%
% Empty !
%
% Msc1
%
% Empty !
%
% Mse
%
Mse(3:6,1:4)=eye(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%            Solution							                 %
%                                                             %
%           -1                                                %
% X(t) = Mcc  Mcs S(t)                                        %
%                                                             %
%                      -1             -1                      %
% S(t+1)=(Mss0-Msc0 Mcc  Mcs)(Msc1 Mcc  Mcs-Mss1)S(t)         %
%                       -1                                    %
%        +(Mss0-Msc0 Mcc  Mcs)e(t+1)                          %
%                                                             %
% S(t+1) = W S(t) + R e(t+1)                                  %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M0=inv(Mss0-Msc0*inv(Mcc)*Mcs);
M1=(Mss1-Msc1*inv(Mcc)*Mcs);
W=-M0*M1;
%
% Calculation of eigenvalues : MU -> Eigenvalues, P -> Eigenvectors
%
[P,MU] = eig(W);

% Sorting the eigenvalues

AMU=abs(MU);
%

%Rearrange system according to the absolute value of eignevalues

[MU,k] = sort(diag(AMU));
P=P(:,k);
Q=inv(P);
%
% Solution
% Z = backward states
% Z(t+1) = MSS*Z(t) + MSE*e(t+1)
%
M=[eye(nback);-inv(Q(nback+1:nstat,nback+1:nstat))*Q(nback+1:nstat,1:nback)];

MSS=W(1:nback,:)*M;
M2=M0*Mse;
MSE=M2(1:nback,1:nshoc);
% Solution for controls (+ prices)
PIC=inv(Mcc)*Mcs*M;

VARIAB=[
   'r1 ';
   'r2 ';
   'R1 ';
   'R2 ';
   'q  ';
   'y1 ';
   'y2 ';   
   'm  ';
   'p1 ';
   'p2 ';
   ];

%%eval(['delete ',outfile])
%%eval(['diary ',outfile])


disp('Flexible exchange rate & fixed prices ')
disp(' ')
disp('      p1      p2       g1     g2        e1      e2 ')
disp('     -----------------------------------------------------')

for i=1:ncont
   disp([sprintf('%s: ',VARIAB(i,:)) sprintf('%5.4f ',PIC(i,:))])
end;


SEE=MSE*SEE*SEE*MSE';
VAR=PIC*reshape(inv(eye(size(kron(MSSC,MSSC)))-kron(MSSC,MSSC))*SEE(:),nback,nback)*PIC';
VY=diag(VAR);

%%select=[7:18];
%%VY=diag(VAR(select,select));
%%CV=VAR(select,select)./sqrt(VY*VY');
%%disp(' ')
%%disp(' q12    q13    q23     y1     y2     y3    s12    s13    s23     p1     p2     p3')
%%disp('----------------------------------------------------------------------------------')
%%disp(sprintf('%5.5f  ',VY'))

disp(' ')
disp(' r1    r2     R1      R2      q       y1       y2   m    p1   p2' )
disp('-----------------------------------------------------------------')
disp(sprintf('%5.5f  ',VY'))

diary off

ans=1;
while ans~=5;
	ans = menu('Impulse Response Function','fiscal shock 1','fiscal shock 2', ...
      'Money shock 1','Money shock 2','Quit');   
   close
	if ans~=5;
		CA=zeros(nshoc,1);
      CA(ans)=1;									% initial value of the impulse				
      REPSA=zeros(nback,nrep);				% IRF matrix of state variables
		REPSA(:,1)=MSE*CA;						% Initialize the IRF matrix of state variables
		for i=2:nrep;								%
		    REPSA(:,i)=MSSC*REPSA(:,i-1);		% Main loop for IRF
		end;											%
      REPXA=PIC*REPSA;							% IRF matrix of controls
		subplot(321);plot(REPXA(6:7,:)');title('Output');legend('Y^1','Y^2')
		subplot(322);plot(REPSA(1:2,:)');title('Price');legend('P^1','P^2')
      subplot(323);plot(REPXA(3:4,:)');title('Nominal Interest Rate');legend('R^1','R^2')
      subplot(324);plot(REPXA(5,:)');title('Terms of trade');legend('q')
		subplot(325);plot(REPXA(1:2,:)');title('Real interest rate');legend('r^1','r^2')
		subplot(326);plot(REPXA(8,:)');title('money');legend('m')
	end;
end;

