%dynerbsflexprices.m

%Chapter 7

%CIA model with predetermined exchange rates

%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
format compact


nc = 4;        % cT, cN, i, e
ns = 1;        % k
ncs = 1;       % lambda
ne = 1;        % epsilon                   

% ECONOMIC PARAMETER VALUES.

r = 0.015
Beta = 1/(1+r)
k = 5
yN = 1
yT = 1
epsilon = 0.5


% STEADY-STATE CALCULATIONS

disp('steady state')

cN = yN
cT = r*k + yT
piN = epsilon
e = yN/cT
i = (1+r)*(1+epsilon)-1
lambda =1/(cT*(1+i))
TB = yT - cT


disp('pause')
pause

%  Check steady state

error1 = cN - yN
error2 = e- cN/cT
error3 = (1/cT) - lambda*(1+i) 
error4 = (1+i)-(1+r)*(1+epsilon)
error5 = r*k + yT - cT
error6 = yT - cT - TB


disp('pause')

pause

% BASIC SYSTEM MATRICES:

% First, specify the matrices in control equations

Mcc = zeros(nc,nc);

%	cT	cN	  i    e
%	1	2	  3     4 

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

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

Mcc(3,3) = 1;

Mcc(4,2) = 1; 


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

Mcs(1,2) = 1;


Mce = zeros(nc,ne);

Mce(3,1) = (epsilon/i)*(1+r);



%       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,1) = 1; 
Mss1(2,1) = -(1+r);


%	cT	cN	   i   e
%	1	2	  3    4  


Msc0 = zeros(ns+ncs,nc);

Msc1 = zeros(ns+ncs,nc);


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




%   epsilon        
%	1


Mse0 = zeros(ns+ncs,ne);

Mse1 = zeros(ns+ncs,ne);


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

% Order of flow variables: TB, rd


nf=1; %TB, rd

%	cT	cN	i   
%	1	2	3	


FVc = zeros(nf,nc);

FVc(1,1) = cT/TB;


%       lambda    m     piN     
%       1         2	     3	       

FVke= zeros(nf,ns+ne);

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; 


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

format