clc;

close all;

%----------------------------------------------------------------
% Defining variables
%----------------------------------------------------------------

% endo 41
var Rss Rstarss A C I G R Pd W S T Pmstar Rstar ldr X M Yd Yx Kd Kx K Md Mx Ld Lx L Qd Px Qx Qk Pm Ql Rk NGDP B Lambda Pxstar Rl dPd dW dS dPxstar dPmstar;

% exo 12
varexo ea ec ei eg er epd ew es et epmstar erstar eldr;

% param 28
parameters alpha beta delta eta gammadl gammadm gammaxl gammaxm;
parameters mud muw phil pie piestar psi sigma chi kappa rhoa rhog;
parameters rhor rhorstar rhot tau upsilon xib xid xii xiw; 

%----------------------------------------------------------------
% Calibration
%----------------------------------------------------------------

alpha = log(1.024)/4;beta = 1.03^(-0.25);delta = 0.0105;eta = 1/0.33;
gammadl = 0.70;gammadm = 0.15;gammaxl = 0.64;gammaxm = 0.18;mud = 1.2; 
muw = 1.05;phil = 1;pie = log(1.03)/4;piestar = log(1.03)/4;psi = 0.25; 
sigma = 0.2;chi = 0.85;kappa = 25;rhoa = 0.8;rhog = 0.8;rhor = 0.85; 
rhorstar = 0.8;rhot = 0.8;tau = 0.5;upsilon = log(1.02)/4; 
xib = 0.4;xid = 10;xii = 1;xiw = 10;


%----------------------------------------------------------------
% Model 
%----------------------------------------------------------------

model; 

Pd*Lambda = (1-chi)/(C-chi*exp(alpha+ec)*C(-1)); 
Lambda = beta*exp(R)*Lambda(1); 
phil*L^eta = Lambda*Ql; 
log(muw*Ql)-log(W)=xiw*[(dW-dW(-1))-beta*(dW(1)-dW)-ew];  
Qk = beta*Lambda(1)/Lambda*(Rk(1)+(1-delta)*Qk(1)); 
Qk/Pd = exp(ei)*1/(1-xii*(I/I(-1))*(I/I(-1)-(1+alpha))-0.5*xii*(I/I(-1)-(1+alpha))^2)*(1-beta*xii*(Lambda(1)/Lambda)*(Qk(1)/Pd)*(I(1)/I)^2*(I(1)/I-(1+alpha)));  
K = (1-delta)*K(-1) + I; 

Yd = (Kd)^(1-gammadl-gammadm)*(A*Ld)^gammadl*(Md)^gammadm;  
gammadl*Qd*Yd = (1+Rl)*W*Ld;
gammadm*Qd*Yd = Pm*Md;
(1-gammadl-gammadm)*Qd*Yd = Rk*Kd;
log(mud*Qd)-log(Pd)=xid*[(dPd-dPd(-1))-beta*(dPd(1)-dPd)-epd];

Yx  = (Kx)^(1-gammaxl-gammaxm)*(A*Lx)^gammaxl*(Mx)^gammaxm;   
gammaxl*Qx*Yx =W*Lx;
gammaxm*Qx*Yx = Pm*Mx;
(1-gammaxl-gammaxm)*Qx*Yx = Rk*Kx;
Qx=Px ;  

ldr = ((Yd/(Yd(-1)*(1+alpha)))^tau)*exp(eldr);
Rl = (1/ldr)*R; 

B = B(-1)*exp(Rstar(-1)+dS+xib*(B(-1)/(4*NGDP(-1))-psi)+upsilon)-(Px*X-Pm*M);  

R-Rstar = dS(1)+xib*(B/(4*NGDP)-psi)+upsilon+es; 

L = Ld + Lx;
M = Md + Mx;
K(-1) = Kd + Kx;  
Yd = C+G+I;
Yx = X;

log(Pd*G)=rhog*log(Pd(-1)*G(-1))+(1-rhog)*log(sigma*NGDP)+eg;
 
R=rhor*R(-1)+(1-rhor)*(Rss+kappa*(dPd(1)-pie))+er;

Px = Pxstar*S;
Pm = Pmstar*S;
Pxstar=T*Pmstar;

dPmstar=piestar+epmstar;
log(T)=rhot*log(T(-1))+et;
  
Rstar = rhorstar*Rstar(-1)+(1-rhorstar)*Rstarss + erstar;    

dPd = log(Pd)- log(Pd(-1));
dS = log(S)- log(S(-1));
dW = log(W)- log(W(-1));
NGDP = Pd*(C+G+I)+Px*X-Pm*M; 
dPxstar=log(Pxstar)-log(Pxstar(-1));
dPmstar=log(Pmstar)-log(Pmstar(-1));

log(A)=log(A(-1))+alpha+ea;

Rss = (-log(beta)+ alpha + pie);
Rstarss = (-log(beta)+piestar+alpha-upsilon);    

end;


steady;
check;