%      PROGRAM:       STUDYGR.M
%
%      DESCRIPTION:   solves a two-country sticky prices
%                     DSGE monetary model
%
%      ASSUMPTIONS:   i) Home and Foreign goods are aggregated by means of
%                        the Cobb-Douglas consumption index
%                     ii) Calvo-type staggered price setting
%                     iii) nominal money growth rates are AR(1) in both countries
%
%      FURTHER READING:    Obstfeld and Rogoff's Redux model
%

% Setting the parameters:

alpha=0.7 ;      % Probability of not changing prices
beta=0.99 ;      % Discount factor
epsilon=1.2 ;    % Elasticity of marginal utility of money
k=2 ;            % Parameter of utility function
theta=7.88 ;     % Degree of monopolistic comp. in the labour mkt
n=0.5 ;          % Size of the Home country

rho=0.9 ;        % Autocorrelation coefficient Home
srho=0.9 ;       % Autocorrelation coefficient Foreign
eps=0.9 ;          % Home money shock. Units: percent
seps=0 ;         % Foreign money shock. Units: percent


% Steady state level of consumption:

C0=((theta-1)/(theta*k))^(1/2);
   % in this model C0=Y0 (output). We could parametrise k so that the steady state level of hours
   % is equal to 0.25, but this will not affect the properties of the model, because we could 
   % have normalised the level of bonds holdings (dB=dB/C0) and C0 would have disappered from
   % the system.


% The system is:
%     PSI*y(t+1)=PHI*y(t)
%
% y(t)=[dB(t-1) C(t) SC(t) M(t-1) SM(t-1) PI(t) SPI(t) I(t) SI(t) R(t-1) ...
%       Ex_rate(t-1) X(t) T(t-1) MU(t) SMU(t) R(t)]'
% X is defined as follows: X(t+1)=Ex_rate(t+1)
%
% 8 States: dB(t-1), M(t-1), SM(t-1), R(t-1), Ex_rate(t-1), T(t-1), MU(t), SMU(t)
% 8 Free variables: C(t), SC(t), PI(t), SPI(t), I(t), SI(t), X(t), R(t)
      
   
% Declaring the matrices:

PSI=zeros(16,16);
PSI(1,1)=1;
PSI(2,2)=1;
PSI(3,3)=1;
PSI(4,4)=1;
PSI(8,4)=1;
PSI(5,5)=1;
PSI(9,5)=1;
PSI(6,6)=n;
PSI(7,6)=n;
PSI(11,6)=-alpha/((1-alpha)*(1-alpha*beta))*beta;
PSI(6,7)=1-n;
PSI(7,7)=1-n;
PSI(12,7)=-alpha/((1-alpha)*(1-alpha*beta))*beta;
PSI(2,10)=beta-1;
PSI(3,10)=beta-1;
PSI(6,10)=1-beta;
PSI(7,10)=1-beta;
PSI(6,11)=n-1;
PSI(7,11)=n;
PSI(8,11)=1-n;
PSI(9,11)=-n;
PSI(10,11)=1;
PSI(13,11)=-1;
PSI(6,12)=1-n;
PSI(7,12)=-n;
PSI(11,13)=-2*(1-n)/(1+theta);
PSI(12,13)=2*n/(1+theta);
PSI(13,13)=1;
PSI(14,14)=1;
PSI(15,15)=1;
PSI(16,10)=1;


PHI=zeros(16,16);
PHI(1,1)=1/beta;
PHI(1,2)=-(1-n)*C0;
PHI(2,2)=1;
PHI(4,2)=1/epsilon;
PHI(11,2)=2/(1+theta);
PHI(1,3)=(1-n)*C0;
PHI(3,3)=1;
PHI(5,3)=1/epsilon;
PHI(12,3)=2/(1+theta);
PHI(8,4)=1;
PHI(9,5)=1;
PHI(8,6)=-n;
PHI(9,6)=-n;
PHI(11,6)=-alpha/((1-alpha)*(1-alpha*beta));
PHI(13,6)=-1;
PHI(8,7)=n-1;
PHI(9,7)=n-1;
PHI(12,7)=-alpha/((1-alpha)*(1-alpha*beta));
PHI(13,7)=1;
PHI(4,8)=-beta/epsilon;
PHI(6,8)=1-beta;
PHI(5,9)=-beta/epsilon;
PHI(7,9)=1-beta;
PHI(8,11)=1-n;
PHI(9,11)=-n;
PHI(13,11)=-1;
PHI(10,12)=1;
PHI(13,13)=1;
PHI(8,14)=1;
PHI(14,14)=rho;
PHI(9,15)=1;
PHI(15,15)=srho;
PHI(16,16)=1;


[VV,DD] = eig(PHI,PSI) ;

% You should get 2 eigenvalues =1, 6 eigenvalues <1 and 8 eigenvalues >1


[DD_sort,index]=sort(abs(diag(DD)));
VV_inv=inv(VV) ;
VV_sort=VV_inv(index,1:16) ;


% Now we find a policy function for C(t), SC(t), PI(t), SPI(t),
%                                   I(t), SI(t), X(t), R(t)
% by solving:
% AA*[C(t) SC(t) PI(t) SPI(t) I(t) SI(t) X(t) R(t)]'=
% BB*[dB(t-1) M(t-1) SM(t-1) R(t-1) Ex_rate(t-1) T(t-1) MU(t) SMU(t)]'
% where:
% X(t+1)=Ex_rate(t+1)


AA=[VV_sort(9:16,2:3), VV_sort(9:16,6:9), VV_sort(9:16,12), VV_sort(9:16,16)];
BB=-[VV_sort(9:16,1), VV_sort(9:16,4:5), VV_sort(9:16,10:11), VV_sort(9:16,13:15)];
PP=inv(AA)*BB;    


% IMPULSE RESPONSES

PER=100;
dF=zeros(PER+1,1);           % starts from t-1
C=zeros(PER,1);              % starts from t
SC=zeros(PER,1);             % starts from t
M=zeros(PER+1,1);            % starts from t-1
SM=zeros(PER+1,1);           % starts from t-1
PI=zeros(PER,1);             % starts from t
SPI=zeros(PER,1);            % starts from t
I=zeros(PER,1);              % starts from t
SI=zeros(PER,1);             % starts from t
R=zeros(PER+1,1);            % starts from t-1
Ex_rate=zeros(PER+1,1);      % starts from t-1
T=zeros(PER+1,1);            % starts from t-1
MU=zeros(PER+1,1);           % starts from t
SMU=zeros(PER+1,1);          % starts from t

MU(1)=eps;
SMU(1)=seps;

t=1;
for t=1:PER,
   y=[dF(t), M(t), SM(t), R(t), Ex_rate(t), T(t), MU(t), SMU(t)]';
   z=PP*y;
   C(t)=z(1);
   SC(t)=z(2);
   PI(t)=z(3);
   SPI(t)=z(4);
   I(t)=z(5);
   SI(t)=z(6);
   Ex_rate(t+1)=z(7);
   dF(t+1)=1/beta*dF(t)+(1-n)*C0*(SC(t)-C(t));
   T(t+1)=T(t)+Ex_rate(t+1)-Ex_rate(t)-PI(t)+SPI(t);
   M(t+1)=M(t)-n*PI(t)-(1-n)*SPI(t)-(1-n)*(Ex_rate(t+1)-Ex_rate(t))+MU(t);
   SM(t+1)=SM(t)-n*PI(t)-(1-n)*SPI(t)+n*(Ex_rate(t+1)-Ex_rate(t))+SMU(t);
   R(t+1)=z(8);
   MU(t+1)=rho*MU(t);
   SMU(t+1)=srho*SMU(t);
end;

Time=(1:PER)';
Response=[dF(2:(PER+1)), C, SC, M(2:(PER+1)), SM(2:(PER+1)), ...
          PI, SPI, I, SI, R(2:(PER+1)), Ex_rate(2:(PER+1)), ...
          T(2:(PER+1)), MU(1:PER), SMU(1:PER)];
    
  
names=['MU '
       'SMU'];
full=[': Nominal Money Growth Rate Home   '
      ': Nominal Money Growth Rate Foreign'];
plot(Time, Response(:,13:14)), title(['Impulse responses']), ...
   xlabel('quarters'), ylabel('% deviations');
for k=1:2, text(4,Response(3,k+12),names(k,:)); end;
disp([' ']);
disp(['Variable names:']);
disp([names, full]);  
pause;
  
names=['Bond'
       'C H '
       'C F '];
full=['                     '
      ': Home Consumption   '
      ': Foreign Consumption'];
plot(Time, Response(:,1:3)), title(['Impulse responses']), ...
   xlabel('quarters'), ylabel('% deviations');
for k=1:3, text(4,Response(3,k),names(k,:)); end;
disp([' ']);
disp([names, full]);
pause;

names=['RMoney H'
       'RMoney F'];
full=[': Real Money Balances Home   '
      ': Real Money Balances Foreign'];
plot(Time, Response(:,4:5)), title(['Impulse responses']), ...
   xlabel('quarters'), ylabel('% deviations');
for k=1:2, text(4,Response(3,k+3),names(k,:)); end;
disp([' ']);
disp([names, full]);
pause;

names=['Infl H'
       'Infl F'];
full=[': Inflation Home   '
      ': Inflation Foreign'];
plot(Time, Response(:,6:7)), title(['Impulse responses']), ...
   xlabel('quarters'), ylabel('% deviations');
for k=1:2, text(4,Response(3,k+5),names(k,:)); end;
disp([' ']);
disp([names, full]);
pause;

names=['Exrate'
       'TT    '];
full=[': Nominal Exchange Rate'
      ': Terms of Trade       '];
plot(Time, Response(:,11:12)), title(['Impulse responses']), ...
   xlabel('quarters'), ylabel('% deviations');
for k=1:2, text(4,Response(3,k+10),names(k,:)); end;
disp([' ']);
disp([names, full]);
pause;

names=['i_t H'
       'i_t F'
       'r_t  '];
full=[': Nominal Interest Rate Home   '
      ': Nominal Interest Rate Foreign'
      ': Real Interest Rate           '];
plot(Time, Response(:,8:10)), title(['Impulse responses']), ...
   xlabel('quarters'), ylabel('% deviations');
for k=1:3, text(10,Response(3,k+7),names(k,:)); end;
disp([' ']);
disp([names, full]);




% Laura Povoledo
% 23 January 2003