%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THIS CODE SOLVES THE QUANTITATIVE MODEL FROM SECTION 4 OF "EXCHANGE RATE   %%
%% DISCONNECT IN GENERAL EQUILIBRIUM" BY O.ITSKHOKI AND D.MUKHIN AND COMPUTES %%
%% SIMULATED MOMENTS FROM TABLES 1-4 FOR DIFFERENT CALIBRATIONS OF THE MODEL  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%% PARAMETERS
clear
sigma = 2; %inverse of IES
phi = 1; %inverse of Frisch elasticity
f1 = 0.5; %intermediate share
theta = 1.5; %elasticity of substitution between H and F goods
alpha = 0.4; %strategic complementarity
gamma = 0.07; %Home trade openness 
gamma_f = 0.07; %Foreign trade openness 
eps = 4; %elasticity of labor substitution

beta = 0.99; %discount factor
rho = 0.97; %persistence of psi-shocks
rho_a = 0.97; %persistence of a-shocks
rho_xi = 0.97; %persistence of xi-shocks
rho_r = 0.95; %interest rate smoothing
delta_p = 2.15; %Taylor rule coefficient

f2 = (1-f1)*0.3; %capital share
kappa = 6.8; %capital adjustment costs
delta = 0.02; %capital depreciation

lambda_w = 0.85; %Calvo parameter for wages (prob. of non-adjustment)
lambda_p = 0.75; %Calvo parameter for prices (prob. of non-adjustment)
currency = [2;2]; %export invoicing of H and F (1,1 - PCP, 2,2 - LCP, 2&1 - DCP)

xi = 0.001; %financial sector parameter

vol = [0;2.5;1]; %relative std of shocks (m,a,xi,psi)
crl = [0.56;0.37]; %cross-country correlation for (m,a)

B = 1000; %number of simulations
TT = 120; %length of simulations



%% STATIC PART
MS1 = eye(16,16); %system for {p,p*,mc,mc*,l,l*,y,y*,yh,yh*,yf,yf*,rk,rk*,r,r*}
MS2 = zeros(16,26);
MS3 = zeros(16,5);

d = 1-delta*f2/(1-f1)/(1/beta-1+delta); %SS share of C in GDP

MS2(1,[1 3]) = [1-gamma gamma]; %p

MS2(2,[2 4]) = [gamma_f 1-gamma_f]; %p*

MS1(3,[1 13]) = [-f1-f2 -f2]; %mc
MS2(3,[5 24]) = [1-f1-f2 -1];

MS1(4,[2 14]) = [-f1-f2 -f2]; %mc*
MS2(4,[6 25]) = [1-f1-f2 -1];

MS1(5,[3 5 7]) = [1 -1 1]; %l
MS2(5,5) = 1;

MS1(6,[4 6 8]) = [1 -1 1]; %l*
MS2(6,6) = 1;

MS1(7,9:10) = -[1-gamma gamma]; %y

MS1(8,11:12) = -[gamma_f 1-gamma_f]; %y*

MS1(9,[1 3 7]) = -[theta-f1 f1 f1]; %yh
MS2(9,[1 7 9]) = [-theta (1-f1)*d (1-f1)*(1-d)];

MS1(10,[2 4 8]) = -[theta-f1 f1 f1]; %yh*
MS2(10,[2 8 10]) = [-theta (1-f1)*d (1-f1)*(1-d)];

MS1(11,[1 3 7]) = -[theta-f1 f1 f1]; %yf
MS2(11,[3 7 9]) = [-theta (1-f1)*d (1-f1)*(1-d)];

MS1(12,[2 4 8]) = -[theta-f1 f1 f1]; %yf*
MS2(12,[4 8 10]) = [-theta (1-f1)*d (1-f1)*(1-d)];

MS1(13,[1 3 7 13]) = [-1 1 1 -1]; %rk
MS2(13,12) = 1;

MS1(14,[2 4 8 14]) = [-1 1 1 -1]; %rk*
MS2(14,13) = 1;

MS1(15,15) = 1; %r
MS2(15,15) = rho_r;
MS2(15,[1 3 17 19]) = [(1-rho_r)*delta_p*[1-gamma gamma -(1-gamma) -gamma]];
MS3(15,1) = 1;

MS1(16,16) = 1; %r*
MS2(16,16) = rho_r;
MS2(16,[4 2 20 18]) = [(1-rho_r)*delta_p*[1-gamma gamma -(1-gamma) -gamma]];
MS3(16,2) = 1;

MS = MS1\MS2; %solution to static variables
MSin = MS1\MS3;

X = zeros(12,26);
X(1,11) = 1; %NER
X(2,:) = X(1,:)+MS(2,:)-MS(1,:); %CPI-RER
X(3,:) = MS(1,:)-MS(2,:); %p-p*
X(4,[2 3 11]) = [-1 1 -1]; %ToT
X(5,7:8) = [1 -1]; %Ó-Ó*
X(10,9:10) = [1 -1]; %i-i*
X(6,:) = d*X(5,:)+(1-d)*X(10,:)+2*X(7,:); %gdp-gdp*
X(7,:) = gamma/(1-f1)*(MS(10,:)-MS(11,:)-X(4,:)); %nx
X(8,:) = -gamma_f/(1-f1)*(MS(10,:)-MS(11,:)-X(4,:)); %nx*
X(9,15:16) = [1 -1]; %r-r* (previous period)
X(11,:) = X(6,:)-d*X(5,:); %s-s* (savings)
X(12,9) = 1; %i
X(13,[7 9]) = [d 1-d]; %gdp
X(13,:) = X(13,:)+X(7,:); %gdp
X(14,7) = 1; %c
X(15,10) = 1; %i*
X(16,[8 10]) = [d 1-d]; %gdp*
X(16,:) = X(16,:)-X(7,:); %gdp*
X(17,8) = 1; %c*
X(18,:) = MS(5,:); %l
X(19,:) = MS(6,:); %l*
X(20,14) = 1; %b
X(21,:) = MS(1,:); %p
X(22,[5 6 11]) = [-1 1 1]; %wage-RER
X(23,[1 4 11]) = [-1 1 1]; %PPI-RER



%% DYNAMIC PART
MD1 = eye(26); %system for x={ph,ph*,pf,pf*,w,w*,c,c*,i,i*,e,k,k*,b,r_,r*_,ph_,ph*_,pf_,pf*_,w_,w*_,e_,a,a*,psi} 
MD2 = zeros(26,26); %MD1*x_{t+1} = MD2*x_{t}+MD3*epsilon_{t+1}
MD3 = zeros(26,5); %epsilon={em,em*,ea,ea*,epsi} - current {em,em*}, but future {ea,ea*,epsi}!

kp = (1-beta*lambda_p)*(1-lambda_p)/lambda_p;
kw = (1-beta*lambda_w)*(1-lambda_w)/lambda_w/(1+eps/phi); 

MD1(1,1) = -beta; %ph 
MD2(1,[1 17]) = [-(1+beta+kp) 1]; 
MD2(1,:) = MD2(1,:)+kp*((1-alpha)*MS(3,:)+alpha*MS(1,:)); 

if currency(1) == 1
    MD1(2,[2 11]) = -beta*[1 1]; %ph* (PCP)
    MD2(2,[2 11 18 23]) = [-(1+beta+kp) -(1+beta+kp*(1-alpha)) 1 1];
    MD2(2,:) = MD2(2,:)+kp*((1-alpha)*MS(3,:)+alpha*MS(2,:));
else
    MD1(2,2) = -beta; %ph* (LCP/DCP)
    MD2(2,[2 11 18]) = [-(1+beta+kp) -kp*(1-alpha) 1];
    MD2(2,:) = MD2(2,:)+kp*((1-alpha)*MS(3,:)+alpha*MS(2,:));
end 

if currency(2) == 1
    MD1(3,[3 11]) = -beta*[1 -1]; %pf (PCP/DCP)
    MD2(3,[3 11 19 23]) = [-(1+beta+kp) (1+beta+kp*(1-alpha)) 1 -1];
    MD2(3,:) = MD2(3,:)+kp*((1-alpha)*MS(4,:)+alpha*MS(1,:));
else 
    MD1(3,3) = -beta; %pf (LCP)
    MD2(3,[3 11 19 23]) = [-(1+beta+kp) kp*(1-alpha) 1 0];
    MD2(3,:) = MD2(3,:)+kp*((1-alpha)*MS(4,:)+alpha*MS(3,:));
end

MD1(4,4) = -beta; %pf*
MD2(4,[4 20]) = [-(1+beta+kp) 1]; 
MD2(4,:) = MD2(4,:)+kp*((1-alpha)*MS(4,:)+alpha*MS(2,:)); 

MD1(5,5) = -beta; %w 
MD2(5,[5 7 21]) = [-(1+beta+kw) kw*sigma 1]; 
MD2(5,:) = MD2(5,:)+kw*(1/phi*MS(5,:)+MS(1,:)); 

MD1(6,6) = -beta; %w*
MD2(6,[6 8 22]) = [-(1+beta+kw) kw*sigma 1]; 
MD2(6,:) = MD2(6,:)+kw*(1/phi*MS(6,:)+MS(2,:)); 

MD1(7,7) = sigma; %c
MD1(7,:) = MD1(7,:)+MS(1,:);
MD2(7,7) = sigma; 
MD2(7,:) = MD2(7,:)+MS(1,:)+MS(15,:);
MD3(7,:) = MSin(15,:);

MD1(8,8) = sigma; %c*
MD1(8,:) = MD1(8,:)+MS(2,:);
MD2(8,8) = sigma; 
MD2(8,:) = MD2(8,:)+MS(2,:)+MS(16,:);
MD3(8,:) = MSin(16,:);

MD1(9,[7 9 12]) = [-sigma beta*kappa*delta -beta*kappa*delta]; %i
MD1(9,:) = MD1(9,:)+beta*(1/beta-1+delta)*MS(13,:);
MD2(9,[7 9 12]) = [-sigma kappa*delta -kappa*delta]; 

MD1(10,[8 10 13]) = [-sigma beta*kappa*delta -beta*kappa*delta]; %i*
MD1(10,:) = MD1(10,:)+beta*(1/beta-1+delta)*MS(14,:);
MD2(10,[8 10 13]) = [-sigma kappa*delta -kappa*delta]; 

MD1(11,[11 14]) = [1 -xi]; %e
MD2(11,[11 26]) = [1 -1]; 
MD2(11,:) = MD2(11,:)+MS(15,:)-MS(16,:);
MD3(11,:) = MSin(15,:)-MSin(16,:);

MD1(12,12) = 1; %k
MD2(12,[9 12]) = [delta 1-delta]; 

MD1(13,13) = 1; %k*
MD2(13,[10 13]) = [delta 1-delta]; 

MD1(14,14) = beta; %b
MD2(14,14) = 1;
MD2(14,:) = MD2(14,:)+X(7,:); 

MD1(15,15) = 1; %r
MD2(15,:) = MS(15,:); 
MD3(15,:) = MSin(15,:);

MD1(16,16) = 1; %r*
MD2(16,:) = MS(16,:); 
MD3(16,:) = MSin(16,:);

MD1(17:23,17:23) = eye(7); %lagged prices, wages and ER
MD2(17:22,1:6) = eye(6);
if sum(currency)==4
    MD2(23,23) = 1; %redundant state under LCP
else
    MD2(23,11) = 1;
end

MD1(24:26,24:26) = eye(3); %shocks
MD2(24:25,24:25) = rho_a*eye(2);
MD2(26,26) = rho;
MD3(24:26,3:5) = eye(3); 



%% BLANCHARD-KAHN
M1 = MD1\MD2; %x_{t+1} = M1*x_{t}+M2*epsilon_{t+1}
M2 = MD1\MD3;
[EVec,EVal] = eig(M1);
ev = diag(abs(EVal));
ind = find(diag(EVal)>1); %find eigenvalues>=1
IEVec = inv(EVec);
M3 = EVal\IEVec*M2; 
Vec = zeros(11,26);
RSI = zeros(11,5);
i = 1;
while i <= 11
    if imag(EVal(ind(i),ind(i))) == 0
        Vec(i,:) = IEVec(ind(i),:);
        RSI(i,:) = M3(ind(i),:);
        i = i+1;
    else
        Vec(i,:) = real(IEVec(ind(i),:));
        Vec(i+1,:) = imag(IEVec(ind(i),:));
        RSI(i,:) = real(M3(ind(i),:));
        RSI(i+1,:) = imag(M3(ind(i),:));
        i = i+2;
    end
end
Vec = real(Vec);
RSI = real(RSI);
RSI(:,3:5) = 0;

RS = -Vec(:,1:11)\Vec(:,12:26); %controls in terms of state variables
RSI = -Vec(:,1:11)\RSI; %controls in terms of innovations
RD = M1(12:26,1:11)*RS+M1(12:26,12:26); %transition matrix for states
RDI = M1(12:26,1:11)*RSI+M2(12:26,:); %transition matrix for innovations



%% SIMULATIONS
Omega = zeros(5); %matrix for shocks
Omega(5,5) = vol(3);
for i = 1:2
    a = sqrt((1-crl(i))/(1+crl(i)));
    Omega(2*(i-1)+(1:2),2*(i-1)+(1:2)) = vol(i)*[1 a;1 -a]*sqrt((1+crl(i))/2); 
end
moments = zeros(22,B);
moments2 = zeros(4,B);
moments3 = zeros(5,B);
moments4 = zeros(5,B);
for b = 1:B
    S_in = Omega*randn(5,2*TT); %innovations
    S_in(:,1:2) = 0; 
    S = zeros(26,2*TT); %state
    Sh = zeros(26,2*TT); %dynamics under financial shocks only
    
    for t = 2:2*TT
        S(12:26,t) = RD*S(12:26,t-1)+RDI*S_in(:,t-1);
        S(1:11,t) = RS*S(12:26,t)+RSI*S_in(:,t);
        
        Sh(12:26,t) = RD*Sh(12:26,t-1)+RDI(:,5)*S_in(5,t-1); 
        Sh(1:11,t) = RS*Sh(12:26,t)+RSI(:,5)*S_in(5,t);
    end
    S_in = S_in(:,TT+1:2*TT); %burn first half of the sample
    S = S(:,TT+1:2*TT); 
    Sh = Sh(:,TT+1:2*TT);
    
    dss = S(:,2:TT)-S(:,1:TT-1); %take time-difference
    ds = Sh(:,2:TT)-Sh(:,1:TT-1); 
    
    dss_cov = cov(dss'); %covariance matrix
    x_cov = X*dss_cov*X';
    
    
    % TABLE 1
    ac = autocorr(X(1,:)*dss);
    moments(1,b) = ac(2); %autocorrelation of d(e)
    moments(2,b) = sqrt(x_cov(1,1)/x_cov(13,13)); %std(e)/std(gdp)
    moments(3,b) = sqrt(x_cov(1,1)/x_cov(14,14)); %std(e)/std(c)
    
    ac = autocorr(X(2,:)*S);
    moments(4,b) = ac(2); %autocorrelation of q
    moments(5,b) = (x_cov(2,2)/x_cov(1,1))^(1/2); %std(q)/std(e)
    moments(6,b) = x_cov(2,1)/(x_cov(1,1)*x_cov(2,2))^(1/2); %cor(q,e)
    moments(7,b) = (x_cov(22,22)/x_cov(1,1))^(1/2); %std(q_w)/std(e)
    moments(8,b) = x_cov(22,1)/(x_cov(1,1)*x_cov(22,22))^(1/2); %cor(q_w,e)
    
    moments(9,b) = x_cov(2,5)/(x_cov(5,5)*x_cov(2,2))^(1/2); %cor(q,c-c*)
    [b_ols,~,~,~,stat_ols] = regress((X(1,:)*dss)',[ones(TT-1,1) (S(15,2:TT)-S(16,2:TT))']);
    moments(10,b) = b_ols(2); %UIP beta
    moments(11,b) = stat_ols(1); %UIP R-squared
    returns = (S(26,1:TT-1)-xi*S(14,2:TT)).*(S(15,2:TT)-S(16,2:TT)-X(1,:)*dss);
    moments(12,b) = mean(returns)/std(returns); %Sharpe ratio
    moments(13,b) = std(S(15,2:TT)-S(16,2:TT))/std(X(1,:)*dss); %std(r-r*)/std(de)
    ac = autocorr(S(15,2:TT)-S(16,2:TT));
    moments(14,b) = ac(2); %autocorrelation of r-r*
    ac = autocorr(S(15,:));
    moments(15,b) = ac(2); %autocorrelation of r
    
    moments(16,b) = sqrt(x_cov(14,14)/x_cov(13,13)); %std(c)/std(gdp)
    moments(17,b) = x_cov(13,14)/sqrt(x_cov(13,13)*x_cov(14,14)); %cor(c,gdp)
    moments(18,b) = x_cov(13,12)/sqrt(x_cov(13,13)*x_cov(12,12)); %cor(i,gdp)
    moments(19,b) = x_cov(13,16)/sqrt(x_cov(13,13)*x_cov(16,16)); %cor(gdp,gdp*)
    moments(20,b) = x_cov(14,17)/sqrt(x_cov(14,14)*x_cov(17,17)); %cor(c,c*)
    moments(21,b) = x_cov(12,15)/sqrt(x_cov(12,12)*x_cov(15,15)); %cor(i,i*)
    moments(22,b) = sqrt(x_cov(12,12)/x_cov(13,13)); %std(i)/std(gdp)
    
    
    % TABLE 2
    b_ols = regress((X(1,:)*ds(:,:))',[ones(TT-1,1) (X(1,:)*dss)']);
    moments2(1,b) = b_ols(2); %psi component in NER
    b_ols = regress((X(2,:)*ds(:,:))',[ones(TT-1,1) (X(2,:)*dss)']);
    moments2(2,b) = b_ols(2); %psi component in RER
    b_ols = regress((X(14,:)*ds(:,:))',[ones(TT-1,1) (X(14,:)*dss)']);
    moments2(3,b) = b_ols(2); %psi component in C
    b_ols = regress((X(13,:)*ds(:,:))',[ones(TT-1,1) (X(13,:)*dss)']);
    moments2(4,b) = b_ols(2); %psi component in GDP
    
    
    % TABLE 3    
    moments3(1,b) = (x_cov(2,2)/x_cov(1,1))^(1/2); %std(q)/std(e)
    moments3(2,b) = (x_cov(4,4)/x_cov(1,1))^(1/2); %std(s)/std(e)
    moments3(3,b) = x_cov(4,1)/(x_cov(1,1)*x_cov(4,4))^(1/2); %cor(s,e)
    moments3(4,b) = (x_cov(7,7)/x_cov(2,2))^(1/2); %std(nx)/std(q)
    moments3(5,b) = x_cov(2,7)/(x_cov(7,7)*x_cov(2,2))^(1/2); %cor(q,nx)
    
    
    % TABLE 4    
    ac = autocorr(X(2,:)*S);
    moments4(1,b) = ac(2); %autocorrelation of q
    moments4(2,b) = (x_cov(2,2)/x_cov(1,1))^(1/2); %std(q)/std(e)
    moments4(3,b) = x_cov(2,5)/(x_cov(5,5)*x_cov(2,2))^(1/2); %cor(q,c-c*)
    moments4(4,b) = sqrt(x_cov(1,1)/x_cov(13,13)); %std(e)/std(gdp)
    moments4(5,b) = (x_cov(7,7)/x_cov(2,2))^(1/2); %std(nx)/std(q)
end



%% TABLES
MM1 = moments;
med1 = median(MM1,2);
std1 = std(MM1,[],2);
rows1 = {'ac(de)';'std(e)/std(gdp)';'std(e)/std(c)';'ac(q)';'std(q)/std(e)';'cor(q,e)';'std(qw)/std(e)';'cor(qw,e)';'cor(q,c-c*)';'beta';'R2';'SR';'std(r-r*)/std(de)';'ac(r-r*)';'ac(r)';'std(c)/std(gdp)';'cor(gdp,c)';'cor(gdp,i)';'cor(gdp,gdp*)';'cor(c,c*)';'cor(i,i*)';'std(i)/std(gdp)'};
table(rows1,med1,std1) %moments from Table 1

MM2 = moments2;
med2 = median(MM2,2);
std2 = std(MM2,[],2);
rows2 = {'NER:psi';'RER:psi';'C:psi';'GDP:psi'};
table(rows2,med2,std2) %moments from Table 2

MM3 = moments3;
med3 = median(MM3,2);
std3 = std(MM3,[],2);
rows3 = {'std(q)/std(e)';'std(s)/std(e)';'cor(s,e)';'std(nx)/std(q)';'cor(q,nx)'};
table(rows3,med3,std3) %moments from Table 3

MM4 = moments4;
med4 = median(MM4,2);
std4 = std(MM4,[],2);
rows4 = {'ac(q)';'std(q)/std(e)';'cor(q,c-c*)';'std(e)/std(gdp)';'std(nx)/std(q)'};
table(rows4,med4,std4) %moments from Table 4



