%dynflexiblerates.m

%Chapter 7

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

%June 2007

%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.  

% Dimension of control vector (nc), state and co-state vectors (ns),
% and forcing (exogenous state) vector.

clc
clear all
close all
format compact


nc = 6;        % cT, cN, i, e, epsilon, c
ns = 1;        % k
ncs = 2;       % lambda, m
ne = 3;        % mu, yT, yN                   

% ECONOMIC PARAMETER VALUES.

r = 0.015
Beta = 1/(1+r)
%k = 5
k=-5     %   !!! be aware that if using k negative, you should check lines 187 and 193 of this code !!!
yN = 1
yT = 1
mu = 0.5
%mu = 0.1
eta = 0.5 % share of c
%sigma = 1.5 %intertemporal elasticity
sigma = 1.5 %intertemporal elasticity
gamma = 0.5


% STEADY-STATE CALCULATIONS

disp('steady state')

cN = yN
cT = r*k + yT
epsilon = mu
e = (gamma/(1-gamma))*(yN/cT)
c = ((cT)^gamma)*(cN^(1-gamma))
i = (1+r)*(1+mu)-1
m = (cT/gamma)*((1-eta)/eta)*((1+i)/i) 
z = m*e^(1-gamma)
lambda = (((c^eta)*(z^(1-eta)))^(-1/sigma))*eta*((z^(1-eta))*(c^(eta-1)))*gamma*((cN/cT)^(1-gamma))
TB = yT - cT




disp('pause')
pause

%  Check steady state

error1 = cN - yN
error2 = e - (gamma/(1-gamma))*(cN/cT)
error3 = lambda -(((c^eta)*(z^(1-eta)))^(-1/sigma))*eta*((z^(1-eta))*(c^(eta-1)))*gamma*((cN/cT)^(1-gamma))
error4 = (1+i)-(1+r)*(1+epsilon)
error5 = r*k + yT - cT
error6 = yT - cT - TB
error7 = m -(cT/gamma)*((1-eta)/eta)*((1+i)/i)
error8 = epsilon - mu


disp('pause')

pause

% BASIC SYSTEM MATRICES:

% First, specify the matrices in control equations

Mcc = zeros(nc,nc);

%	cT	cN	  i    e   epsilon     c
%	1	2	  3     4    5         6

Mcc(1,1) = -(1-gamma);
Mcc(1,2) = (1-gamma);
Mcc(1,5) = (1-gamma)*(1-eta)*((1-sigma)/sigma);
Mcc(1,6) = -(1+eta*(((1-sigma)/sigma)));

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

Mcc(3,1) = -1; 
Mcc(3,4) = 1/(1+i);

Mcc(4,3) = -epsilon/(1+epsilon);
Mcc(4,4) = i/(1+i); 

Mcc(5,2) = 1; 

Mcc(6,1) = -gamma; 
Mcc(6,2) = -(1-gamma); 
Mcc(6,6) = 1; 


Mcs = zeros(nc,ns+ncs);
 
%       k        lambda    m       
%       1          2	   3          

Mcs(1,2) = 1; 
Mcs(1,3) = (1-eta)*((1-sigma)/sigma);

Mcs(3,3) = -1;


Mce = zeros(nc,ne);

Mce(5,3) = 1; 


%       k       lambda            
%       1       2	       


Mss0 = zeros(ns+ncs,ns+ncs);
Mss1 = zeros(ns+ncs,ns+ncs);

Mss0(1,2) = -1; 
Mss1(1,2) = 1; 

Mss0(2,3) = 1; 
Mss1(2,3) = -1; 

Mss0(3,1) = 1; 
Mss1(3,1) = -(1+r); 


%	cT	cN	   i   e
%	1	2	  3    4  


Msc0 = zeros(ns+ncs,nc);

Msc1 = zeros(ns+ncs,nc);


Msc1(2,3) = -epsilon/(1 + epsilon);

Msc1(3,1) = -cT/k;




%   epsilon   yT   yN     
%	1         2     3

Mse0 = zeros(ns+ncs,ne);

Mse1 = zeros(ns+ncs,ne);

Mse1(2,1)= mu/(1+mu); 

Mse1(3,2) = yT/k;





% Some additional matrices links other observable flow variables (nf) to
% controls, states, costates and exogeneous forcing variables. 

% Order of flow variables: TB


nf=1; %TB

FVc = zeros(nf,nc);

%FVc(1,1) = cT/TB;       % for k positive
FVc(1,1) = -cT/TB;       % for k negative

FVke= zeros(nf,ns+ne);

%FVke(1,3) = -yT/TB;     % for k positive
FVke(1,3) = yT/TB;       % for k negative

FVl = zeros(nf,ncs);


% This completes the specification of the model.


% FUNDAMENTAL STATE-COSTATE DIFFERENCE EQUATION


% Derive fundamental difference equation for states and co-states
% The memnonic MS is used to indicate the M* matrices in the text
% Notice that the symbol \ involves a particular concept of matrix
% division (see PCMatlab manual, p. 2-18) in which X=A\B is the solution
% to A*X=B.

MSss0 = Mss0 - Msc0*(Mcc\Mcs);

MSss1 = Mss1 - Msc1*(Mcc\Mcs);


MSse0 = Mse0 + Msc0*(Mcc\Mce);

MSse1 = Mse1 + Msc1*(Mcc\Mce);

% Put the fundamental difference equation in normal form

W = -(MSss0\MSss1);
R = MSss0\MSse0;
Q = MSss0\MSse1;

% Housekeeping

% clear Mss0; clear Mss1; clear Msc0; clear Msc1; clear Mse0; clear Mse1;
% clear MSss0; clear MSss1; clear MSse0; clear MSse1; TO BE UNCOMMENT


% EIGENVECTOR-EIGENVALUE DECOMPOSITION OF STATE TRANSITION MATRIX:

% Compute the eigenvalues and eigenvectors of the matrix W.

[P,MU] = eig(W);
MU

% Order the eigenvalues, in order of ascending absolute value.

AMU=abs(MU);

% Sort by the absolute value of the eigenvectors, then reorder
% the columns of the eigenvector matrix (P) by this sort. (See
% the discussion of sorting in the PCMATLAB manual, p. 3-113,
% where they use this example).

[MU,k] = sort(diag(AMU));
P=P(:,k);

% Next, we need to order the eigenvalue (diagonal) matrix in
% order of ascending absolute value.  As it was unclear how to
% accomplish this, the matrix is simply computed from the standard
% diagonalization formula given the sorted eigenvalue matrix and
% its inverse, denoted PS (matches the technical appendix notation,
% where the inverse of P is denoted with a star).

PS=inv(P);
MU=PS*W*P;


% PARTITIONING MATRICES

% The P, PS, MU, R and Q matrices must be partitioned, as discussed in
% the technical appendix section B.2, for the purpose of computing the
% different components of the solution.  Partitioning
% is undertaken with the colon notation in PCMATLAB, as discussed in
% the manual in the section starting p. 2-31.  In view of potential
% multiple state variable problems, the partitioning uses ns as a
% parameter.

MU1 = MU(1:ns,1:ns);
MU2 = MU(ns+1:ns+ncs,ns+1:ns+ncs);

P11 = P(1:ns,1:ns);
P12 = P(1:ns,ns+1:ns+ncs);
P21 = P(ns+1:ns+ncs,1:ns);
P22 = P(ns+1:ns+ncs,ns+1:ns+ncs);


PS11 = PS(1:ns,1:ns);
PS12 = PS(1:ns,ns+1:ns+ncs);
PS21 = PS(ns+1:ns+ncs,1:ns);
PS22 = PS(ns+1:ns+ncs,ns+1:ns+ncs);


Rke = R(1:ns,1:ne);
Rle = R(ns+1:ns+ncs,1:ne);
Qke = Q(1:ns,1:ne);
Qle = Q(ns+1:ns+ncs,1:ne);

% COMPOSITE EXPRESSIONS

% The solution for the shadow price involves some composite expressions.
% These are:

SP1 = - MU2\(PS21*Rke+PS22*Rle);
SP2 = - MU2\(PS21*Qke+PS22*Qle);

KLK = P11*MU1/P11;
KTL = (P11*MU1*PS12+P12*MU2*PS22)/PS22;

% Note that these expressions involve the use of matrix right division
% in the multiple state variable case (PCMATLAB manual, p. 2-18,19).
% Thus, the formulas appear slightly different in these expressions
% than in King, Plosser, and Rebelo's technical appendix.

% This terminates the program dynmbs.m:  It is necessary at this stage 
% to make decisions concerning the type of experiment that will be studied.
% The output of this program must be

% nc,ns,ne,nf
% Mcc,Mcs,Mce
% FVc,FVke,FVl
% W,R,Q
% MU2,PS21,PS22
% Rke,Qke,
% KTL,KEC,KLK
% SP1,SP2

% Thus, we eliminate all other material from the workspace

clear P; clear PS; clear MU; clear AMU;
clear k;

% and then save the worskpace contents on disk.  This enables multiple uses
% of the outputs of this program, since subsequent programs do housekeeping.

disp('run pathflexiblerates to compute perfect foresight paths')

format