close all
clear all

% 
% This file solves the model presented in the paper "Technology Shocks
% and Employment in Open Economies" by Juha Tervala (juha.tervala@helsinki)
% published in Economics.
%
% If you run this file you can replicate Figure 1 presented in the paper.
%
%================================================
%
% The method developed by Klein (2000) and McCallum (2001) is used to solve the model.
% See Klein, P. (2000). Using the Generalized Schur Form to Solve a
% Multivariate Linear Rational Expectations Model.
% Journal of Economic Dynamics and Control 24 (10): 1405-1423
%
% The model is written in matrix form:
% A*E(t)x(t+1) = B*x(t) + C*Z(t)
% where x(t) is a vector of endogenous variables, x = [y(t)' k(t)']',
% the nk-dimensional vector, k(t) is pretermined and y(t) is
% non-predetermined;
% Z(t) is a vector of exogenous variables determined by an AR(1) process
%
%================================================
%
% This file uses the software deleloped by McCallum (2001).
% See McCallum, B. (2001). Software for RE Analysis. Computer software
% available at
% http://wpweb2.tepper.cmu.edu/faculty/mccallum/research.htm
%
%============= NOTE =============================
%
% The first part of the file solves the LCP model
% The second part of the file solves the PCP model
%
%================================================


%================================================
% The LCP version of the model
%================================================

% Define matrix dimensions

nx = 38; nz = 1; nu = nz;

% Define parameters

beta = 0.99; % the discount factor
epsilon = 1.18;  % 1/epsilon = the consumption elasticity of money demand
theta = 6.0; % the elasticity of substitution between differentiated goods
cnn = 0.5; % the relative country size of the home country
ark = 1.0; % AR(1) parameter, the persistence of technology shocks (1 = permanent)
gamma = 0.75; % the probability of not adjusting prices

ffkk = 1.0; %  Domestic technology shock (1 = one percent shock)
ffkkstar = 0.0; % Foreign technology shock (0 = no shock)

% Number of periods plotted
npts = 10;

% Create matrices

A = zeros(nx,nx);  B = zeros(nx,nx);  C = zeros(nx,nu);
phi = zeros(nu,nu);

% Domestic variables

ic = 1; % consumption
ip = 2; % price index 
ilam = 3; % Tau (See Appendix Page 4)
idelta = 4; % discount factor
ih = 5; %; % labour input
iw = 6; % nominal wage
ibh = 7; % bond holdings of domestic households
ix = 8; % domestic output sold in the foreign country
iv = 9; % domestic output sold in the home countryn
ipz = 10; % domestic currency price of domestic goods
is = 11; % nominal exchange rate
ipzstar = 12; % domestic currency price of Foreign products
iy = 13; % total domestic output
ik = 14;% domestic technology shock
ibz = 15; % Calvo-weighted domestic currency  price of domestuc goods
ibstarz = 16; % Calvo-weighted foriegn currency Calvo price of domestic goods
iens = 17; % exchange rate expectations

% Foreign variables

icstar = 18; % foreign consumption
ipstar = 19; % foreign price index
ilamstar = 20; % Taustar
ideltastar = 21; % foreign discount factor
ihstar = 22; %foreign labour input
iwstar = 23; %foreign wage
ixstar = 24; % foreign output sold in the foreign country
ivstar = 25; % % foreign output sold in the home country
iqzstar = 26; % foreign currency price of foreign goods
iqz = 27; % foreign currency price of domestic goods
iystar = 28; % total foreign output
ikstar = 29;% foreign technlogy shock
ibzstar = 30; % Calvo-weighted domestic currency price of foreign goods
ibstarzstar = 31; % Calco-weighted foriegn currency price of foreign goods

% Predetermined variables

iklag = 32;
iklagstar = 33;
ibhlag = 34;
ibzlag = 35;
ibstarzlag = 36;
ibzstarlag = 37;
ibstarzstarlag = 38;

nk = 7; 

% Shock

iz1 = 1; % Technology shock

% The equations of the model

% Domestic equations

% Define a -ilam =P(hat) + C(hat)
B(1,ic) = 1;
B(1,ip) = 1;
B(1,ilam) = 1;

% Money demand, see Appendix equation (A30)
B(2,ip) = 1-epsilon;
B(2,ilam) = 1/(1-beta);
A(2,ilam) = beta/(1-beta);

% Optimal consumption, see Appendix equation (A29)
B(3,ilam) = -1;
B(3,idelta) = -1;
A(3,ilam) = -1;

% Optimal labour supply (A3) (Makes use of the definition on itau)
B(4,ih) = -1;
B(4,iw) = 1;
B(4,ilam) = 1;

% Consolidated budget contraint of the domestic economy (A27)
B(5,ibh) = -beta;
B(5,ibz) = cnn;
B(5,ix) = cnn;
B(5,ip) = -cnn-(1-cnn);
B(5,iv) = 1-cnn;
B(5,ibstarz) = 1-cnn;
B(5,is) =  (1-cnn);
B(5,ibhlag) =  1;
B(5,ic) =  -1;

% Production function (A7)
B(6,iy) = -1;
B(6,ih) = 1;
B(6,ik) = 1;

% Demand equation (A9)
B(7,ix) = -1;
B(7,ibz) = -theta;
B(7,ip) = theta;
B(7,ic) = 1;

% (A10)
B(8,iv) = -1;
B(8,ibstarz) = -theta;
B(8,ipstar) = theta;
B(8,icstar) = 1;

% Composition of total output (A13)
B(9,iy) = -1;
B(9,ix) = cnn;
B(9,iv) = 1-cnn;

% Domestic Calvo-weighted price indes (A17)
B(10,ip) = -1;
B(10,ibz) = cnn;
B(10,ibzstar) = 1-cnn;

% Uncovered interest rate parity (A28)
B(11,ideltastar) = -1;
B(11,idelta) = 1;
B(11,iens) = 1;
B(11,is) = -1;

% Techology shock (A15)
B(12,ik) = -1;
B(12,iklag) = ark;
C(12,iz1) = ffkk;

% Optinal pricing rule (A23)
B(13,ipz) = -1;
A(13,ipz) = -beta*gamma;
B(13,iw) = 1-beta*gamma;
B(13,ik) = -(1-beta*gamma);

% Calvo-weighted price index of domestic goods sold in the domestic market
% (A19)
B(14,ibz) = -1;
B(14,ibzlag) = gamma;
B(14,ipz) = 1-gamma;

% Optimal pricing rule (A24)
B(15,iqz) = -1;
A(15,iqz) = -beta*gamma;
B(15,iw) = 1-beta*gamma;
B(15,is) = -(1-beta*gamma);
B(15,ik) = -(1-beta*gamma);

% (A20)
B(16,ibstarz) = -1;
B(16,ibstarzlag) = gamma;
B(16,iqz) = 1-gamma;

% Foreign equations

% Define a -itaustar = ipstar)t + icstar
B(17,icstar) = 1;
B(17,ipstar) = 1;
B(17,ilamstar) = 1;

% Money demand
B(18,ipstar) = 1-epsilon;
B(18,ilamstar) = 1/(1-beta);
A(18,ilamstar) = beta/(1-beta);

% Optimal consumption
B(19,ilamstar) = -1;
B(19,ideltastar) = -1;
A(19,ilamstar) = -1;

% Optimal labour supply (A4)
B(20,ihstar) = -1;
B(20,iwstar) = 1;
B(20,ilamstar) = 1;

% Production function (A8)
B(21,iystar) = -1;
B(21,ihstar) = 1;
B(21,ikstar) = 1;

% (A12)
B(22,ixstar) = -1;
B(22,ibstarzstar) = -theta;
B(22,ipstar) = theta;
B(22,icstar) = 1;

% (A11)
B(23,ivstar) = -1;
B(23,ibzstar) = -theta;
B(23,ip) = theta;
B(23,ic) = 1;

% (A14)
B(24,iystar) = -1;
B(24,ixstar) = 1-cnn;
B(24,ivstar) = cnn;

% (A18)
B(25,ipstar) = -1;
B(25,ibstarzstar) = 1-cnn;
B(25,ibstarz) = cnn;

% (A16)
B(26,ikstar) = -1;
B(26,iklagstar) = ark;
C(26,iz1) = ffkkstar;

% (A25)
B(27,iqzstar) = -1;
A(27,iqzstar) = -beta*gamma;
B(27,iwstar) = 1-beta*gamma;
B(27,ikstar) = -(1-beta*gamma);

% (A22)
B(28,ibstarzstar) = -1;
B(28,ibstarzstarlag) = gamma;
B(28,iqzstar) = 1-gamma;

% (A26)
B(29,ipzstar) = -1;
A(29,ipzstar) = -beta*gamma;
B(29,iwstar) = 1-beta*gamma;
B(29,is) = (1-beta*gamma);
B(29,ikstar) = -(1-beta*gamma);

% (A21)
B(30,ibzstar) = -1;
B(30,ibzstarlag) = gamma;
B(30,ipzstar) = 1-gamma;

% Pretermined variables
% For more information how to specify predetermined variables
% see McCallum (2001): Software for RE Analysis, Appendix A, page 6.

% Definition of k(t-1)
A(31,iklag) = 1;
B(31,ik) = 1;

% Definition of k*(t-1)
A(32,iklagstar) = 1;
B(32,ikstar) = 1;

% Define Bh*(t-1)
A(33,ibhlag) = 1;
B(33,ibh) = 1;

% Define b(z)*(t-1)
A(34,ibzlag) = 1;
B(34,ibz) = 1;

% Define b*(z)*(t-1)
A(35,ibstarzlag) = 1;
B(35,ibstarz) = 1;

% Define b*(z*)*(t-1)
A(36,ibstarzstarlag) = 1;
B(36,ibstarzstar) = 1;

% Define b(z*)*(t-1)
A(37,ibzstarlag) = 1;
B(37,ibzstar) = 1;

% Exchange rate expectations (rational expectations)
B(38,iens) = 1;
A(38,is) = 1;

% Coefficients of the AR(1) matrix
phi(iz1,iz1) = 0.0;

% Use solvek.m written by McCallum (2001)

[m,n,p,q,z22h,s,t,lambda] = solvek(A,B,C,phi,nk);

bigmn = [m n];
bigpq = [p q];
bigpsi = eye(nz,nu);

% Calculate impulse response functions
% Use dimpulse.n developed by McCallum 

A = [p q;zeros(nz,nk) phi];
C = [m n];
D = zeros(nx-nk,nu);
B = [zeros(nk,nu);bigpsi];

jj=[0:npts];

% Define shock (iz1 = technology shock)

ishock = iz1;

% Use dimpulse.m

[Y,X]=dimpulse(A,B,C,D,ishock,npts+1);

%================================================
% Plot impulse response functions
%================================================

% Solve the Calvo-weighted domestic terms of trade =
% the relative price of domestic imports in terms of
% domestic exports.
% And note that iw is not needed anymore, use variable iw.
% Overwrite the result of iw

Y(:,iw) = -Y(:,ibstarz)-Y(:,is)+Y(:,ibzstar);

% Solve the welfare effect of a technology shock.
% idelta is not needed anymore, use idelta.
% Domestic welfare effect:

Y(:,idelta) = Y(:,ic)-((theta-1)/theta)*Y(:,ih);

% Foreign welfare effect:

Y(:,ideltastar) = Y(:,icstar)-((theta-1)/theta)*Y(:,ihstar);

% Plot the effects of a domestic technology shock

figure(1)
subplot(5,2,2)
plot(jj,Y(:,ic),'k--')
title('(b) Domestic consumption')

hold on

subplot(5,2,1)
plot(jj,Y(:,iy),'k--')
title('(a) Domestic output')

hold on

subplot(5,2,4)
plot(jj,Y(:,icstar),'k--')
title('(c) Foreign consumption')

hold on

subplot(5,2,3)
plot(jj, Y(:,iystar),'k--')
title('(d) Foreign output')

hold on

subplot(5,2,5)
plot(jj,Y(:,ih),'k--')
title('(e) Domestic employment')

hold on

subplot(5,2,8)
plot(jj,Y(:,iw),'k--')
title('(h) Domestic terms of trade')

hold on
 
subplot(5,2,7)
plot(jj,Y(:,ibh),'k--') 
title('(g) Bond holdings of domestic households')

hold on

subplot(5,2,6)
plot(jj,Y(:,is),'k--') 
title('(f) Nominal exchange rate')

hold on

subplot(5,2,9)
plot(jj,Y(:,idelta),'k--')
title('(i) Domestic welfare')

hold on

subplot(5,2,10)
plot(jj,Y(:,ideltastar),'k--')
title('(j) Foreign welfare')

hold on

clear all

%================================================
%================================================
% The PCP version of the model
%================================================
%================================================

% Define matrix dimensions

nx = 28; nz = 1; nu = nz;

% Define parameters

beta = 0.99;
epsilon = 1.18; 
theta = 6.0;
cnn = 0.5;
ark = 1.0;
gamma = 0.75;

ffkk = 1.0; %  Domestic technology shock (1 = 1 percent shock)
ffkkstar = 0.0; % Foreign technology shock (0 = no shock)

% Number of periods plotted
npts = 10;

% Create matrices

A = zeros(nx,nx);  B = zeros(nx,nx);  C = zeros(nx,nu);
phi = zeros(nu,nu);

% Domestic variables

ic = 1;
ip = 2;
ilam = 3;
idelta = 4;
ih = 5;
iw = 6;
ibh = 7;
ipz = 8; % domestic currency price of domestic goods
is = 9;
iy = 10; % domestic (total) output
ik = 11;
ibz = 12; % Calvo-weighted domestic currency price of domestic goods
iens = 13; % exchange rate expectations

% Foreign variables

icstar = 14;
ipstar = 15;
ilamstar = 16;
ideltastar = 17;
ihstar = 18;
iwstar = 19;
ipzstar = 20; % foreign currency price of foreign goods
iystar = 21; % foreign (total) output
ikstar = 22;
ibstarzstar = 23; % Calvo-weighted foriegn currency price of foreign goods

% Predetermined variables

iklag = 24;
iklagstar = 25;
ibhlag = 26;
ibzlag = 27;
ibstarzstarlag = 28;

nk = 5; 

% Shock

iz1 = 1; % Technology shock

% The equations of the model

% Domestic equations

% Define a -itau =P(hat)t + C(hat)
B(1,ic) = 1;
B(1,ip) = 1;
B(1,ilam) = 1;

% Money demand (A30)
B(2,ip) = 1-epsilon;
B(2,ilam) = 1/(1-beta);
A(2,ilam) = beta/(1-beta);

% Optimal consumption (A29)
B(3,ilam) = -1;
B(3,idelta) = -1;
A(3,ilam) = -1;

% Labour supply (A32), makes use of the definition of tau.
B(4,ih) = -1;
B(4,iw) = 1;
B(4,ilam) = 1;

% Consolidated budget contraint (A49)
B(5,ibh) = -beta;
B(5,ibz) = 1;
B(5,ip) = -1;
B(5,iy) = 1;
B(5,ibhlag) =  1;
B(5,ic) =  -1;

% Production function (A7)
B(6,iy) = -1;
B(6,ih) = 1;
B(6,ik) = 1;

% Demand for domestic goods (A39)
B(7,iy) = -1;
B(7,ibz) = -theta;
B(7,ip) = theta;
B(7,ic) = cnn;
B(7,icstar) = 1-cnn;

% Aggregate price level (A43)
B(8,ip) = -1;
B(8,ibz) = cnn;
B(8,ibstarzstar) = 1-cnn;
B(8,is) = 1-cnn;

% Uncovered interest rate parity (A50)
B(9,ideltastar) = -1;
B(9,idelta) = 1;
B(9,iens) = 1;
B(9,is) = -1;

% Technology (A15)
B(10,ik) = -1;
B(10,iklag) = ark;
C(10,iz1) = ffkk;

% Optimal pricing rule for a domestic good (A47)
B(11,ipz) = -1;
A(11,ipz) = -beta*gamma;
B(11,iw) = 1-beta*gamma;
B(11,ik) = -(1-beta*gamma);

% Calvo-weighted price index of domestic goods (A45)
B(12,ibz) = -1;
B(12,ibzlag) = gamma;
B(12,ipz) = 1-gamma;

% Foreign equations

% Define -itaustar
B(13,icstar) = 1;
B(13,ipstar) = 1;
B(13,ilamstar) = 1;

% Money demand (A35), uses the defintion of ilam (=tau in the appendix)
B(14,ipstar) = 1-epsilon;
B(14,ilamstar) = 1/(1-beta);
A(14,ilamstar) = beta/(1-beta);

% Optimal consumption (A32), uses ilamstar
B(15,ilamstar) = -1;
B(15,ideltastar) = -1;
A(15,ilamstar) = -1;

% (A34)
B(16,ihstar) = -1;
B(16,iwstar) = 1;
B(16,ilamstar) = 1;

% (A38)
B(17,iystar) = -1;
B(17,ihstar) = 1;
B(17,ikstar) = 1;

% (A40)
B(18,iystar) = -1;
B(18,ibstarzstar) = -theta;
B(18,ipstar) = theta;
B(18,icstar) = 1-cnn;
B(18,ic) = cnn;

% (A44)
B(19,ipstar) = -1;
B(19,ibstarzstar) = 1-cnn;
B(19,is) = -cnn;
B(19,ibz) = cnn;

% Technology shock (A42)
B(20,ikstar) = -1;
B(20,iklagstar) = ark;
C(20,iz1) = ffkkstar;

% Pricing rule for a foreign good (A48)
B(21,ipzstar) = -1;
A(21,ipzstar) = -beta*gamma;
B(21,iwstar) = 1-beta*gamma;
B(21,ikstar) = -(beta*gamma);

% Calvo-weighted price index of foreign goods (A46)
B(22,ibstarzstar) = -1;
B(22,ibstarzstarlag) = gamma;
B(22,ipzstar) = 1-gamma;

% Pretermined variables

% Definition of k(t-1)
A(23,iklag) = 1;
B(23,ik) = 1;

% Definition of k*(t-1)
A(24,iklagstar) = 1;
B(24,ikstar) = 1;

% Define Bh*(t-1)
A(25,ibhlag) = 1;
B(25,ibh) = 1;

% Define b(z)*(t-1)
A(26,ibzlag) = 1;
B(26,ibz) = 1;

% Define b*(z*)*(t-1)
A(27,ibstarzstarlag) = 1;
B(27,ibstarzstar) = 1;

% Exchange rate expectations
B(28,iens) = 1;
A(28,is) = 1;

% Coefficients of the AR(1) matrix
phi(iz1,iz1) = 0.0;

% Use solvek.m written by McCallum (2001)

[m,n,p,q,z22h,s,t,lambda] = solvek(A,B,C,phi,nk);

bigmn = [m n];
bigpq = [p q];
bigpsi = eye(nz,nu);


% Calculate impulse response functions
% Use dimpulse.n developed by McCallum 

A = [p q;zeros(nz,nk) phi];
C = [m n];
D = zeros(nx-nk,nu);
B = [zeros(nk,nu);bigpsi];

jj=[0:npts];

% Define shock

ishock = iz1;

% dimpulse.m

[Y,X]=dimpulse(A,B,C,D,ishock,npts+1);

%================================================
% Plot impulse response functions
%================================================

% Solve the Calvo-weighted domestic terms of trade =
% the relative price of domestic imports in terms of
% domestic exports.
% And note that iw is not needed anymore, use variable ix.
% Overwrite the result of iw

Y(:,iw) = -Y(:,ibz)+Y(:,is)+Y(:,ibstarzstar);

% Solve the welfare effect of a technology shock
% idelta is not needed anymore, use idelta.
% Domestic welfare effect:

Y(:,idelta) = Y(:,ic)-((theta-1)/theta)*Y(:,ih);

% Foreign welfare effect:

Y(:,ideltastar) = Y(:,icstar)-((theta-1)/theta)*Y(:,ihstar);

% Plot the effects of a domestic technology shock

figure(1)
subplot(5,2,2)
plot(jj,Y(:,ic),'k-')
legend('LCP','PCP') 

subplot(5,2,1)
plot(jj,Y(:,iy),'k-')
legend('LCP','PCP')

subplot(5,2,4)
plot(jj,Y(:,icstar),'k-')
legend('LCP','PCP')

subplot(5,2,3)
plot(jj, Y(:,iystar),'k-')
legend('LCP','PCP')

subplot(5,2,5)
plot(jj,Y(:,ih),'k-')
legend('LCP','PCP')

subplot(5,2,8)
plot(jj,Y(:,iw),'k-')
legend('LCP','PCP')

subplot(5,2,7)
plot(jj,Y(:,ibh),'k-') 
legend('LCP','PCP')

subplot(5,2,6)
plot(jj,Y(:,is),'k-') 
legend('LCP','PCP')

subplot(5,2,9)
plot(jj,Y(:,idelta),'k-')
legend('LCP','PCP')

subplot(5,2,10)
plot(jj,Y(:,ideltastar),'k-')
legend('LCP','PCP')