%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var
//The first order conditions of the household
   C lambda P pi Rstar H S Q M R I
//Optimal inputs
   MC A W w r K
//Domestic firm's price setting 
    pid Y pid1 Pd Pd1 e1 e2 
//import firm's price setting 
  re Pstar Pm Pm1 pim pim1 f1 f2 IM
  X Ystar Px//Exports
  T G B vd vm
za zg zm zs zystar zpstar zrstar; 
varexo epsa epsg epsm epss epsystar epspstar epsrstar;

parameters 
//Foc of the household
  beta h sigma epsilon PI nu delta omega
//Domestic firm's price setting 
  alpha xid  PId PId1 gamma chid
//import relative price
 xim  PIm PIm1 theta chim phi//Exports
//Government Resource constraint
 a1 a2 a3 a4
//Money Policy Rule
d1 d2 d3 
//market clearing
omicron b1 b2 b3 eta b4 b5 
rhoa rhom rhos rhog rhoystar rhopstar rhorstar;


%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------
 
 beta= 0.99;  
 h=0.65;
 sigma=2.1;
 epsilon=3;
 PI=1.03;
 nu=5; 
 delta=0.35;
 omega=1.58;
 alpha   = 0.33;
 xid=0.85;
 PId=1.02;
 PId1=1.01;
 gamma=1.7; 
 chid=0.65;
 xim=0.83; 
 PIm=1.035;
 PIm1=1.01;
 theta=3.5;
 chim=0.65;
 phi=0.15; 
 a1=0.91;
 a2=0.95;
 a3=0.8;
 a4=2.82;
 d1=0.34;
 d2=0.38;
 d3=0.5;
 omicron=0.33;
 b1=0.39;
 b2=0.43;
 b3=0.14; 
 eta=3.5; 
 b4=0.123;
 b5=0.35;
rhoa=0.75; 
rhom=0.59; 
rhos=0.75; 
rhog=0.75;
rhoystar=0.75;
rhopstar=0.75; 
rhorstar=0.75;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 
//The first order conditions of the household

  beta*h*C(+1)-(1+beta*h^2)*C+h*C-(((1-h)*(1-beta*h))/sigma)*lambda=0;
  (beta/PI)*lambda(+1)-(beta/PI)*pi(+1)-(1-beta/PI)*epsilon*M+(1-beta/PI)
  *epsilon*P-lambda=0;
  lambda+W-P-nu*H=0;
  lambda(+1)-lambda+R-pi(+1)=0;
  lambda(+1)-lambda+Rstar-pi(+1)+S(+1)-S=0;
  (1-beta*(1-delta))*lambda(+1)+(1-beta*(1-delta))*r(+1)+beta*(1-delta)*Q(+1)-Q=0;
  omega*beta*I-(1-beta)*omega*I+Q+omega*I-lambda=0;

//Risk-adjusted uncovered interest parity 

  R-Rstar-S(+1)+S=0;
  pi=P-P(-1);

//Optimal inputs

  K-H-w+r=0;
  MC-(1-alpha)*w-alpha*r+A=0;
  w-W+P=0;

//Domestic firm's price setting 

  (1-beta*xid*PId^gamma*(1-chid))*(lambda+MC+Y)+beta*xid*PId^gamma*(1-chid)
*(gamma*pid(+1)-gamma*chid*pid+e1(+1))-e1=0;
 (1-beta*xid*PId^(1-gamma)*(chid-1))*(lambda+pid1+Y)+beta*xid
*PId^(gamma-1)*(1-chid)*(e2(+1)+pid1-pid1(+1)-(1-gamma)*(pid(+1)-chid*pid))-e2=0;
e1=e2;
pid1=Pd1-Pd;
pid=Pd-Pd(1);
// the domestic price level evolves

((xid*PId^(1-gamma)*(chid-1))/((1-xid)*PId1^(1-gamma)))*(pid-chid*pid(-1))-pid1=0;

//Real exchange rate

  re=S+Pstar-P;
  pim=Pm-Pm(-1);
  pim1=Pm1-Pm;

//import  firm's price setting 
 (1-beta*xim*PIm^theta*(1-chim))*(lambda+re+IM)+beta*xim*PIm^(theta*(1-chim))
*(theta*pim(+1)-theta*chim*pim+f1(+1))-f1=0;

(1-beta*xim*PIm^(1-theta)*(chim-1))*(lambda+pim1+IM)+beta*xim
*PIm^(theta-1)*(1-chim)*((1-theta)*chim*pim-(1-theta)*pim(+1)+f2(+1))-f2=0;
 
// the import price level evolves

((xim*PId^(1-theta)*(chim-1))/((1-xim)*PIm1^(1-theta)))*(pim-chim*pim(-1))-pim1=0;

// Exports (in foreign currency)
  X=phi*(Pstar-Px)+Ystar;
  Px-P-S=0;
//Government Resource constraint

P+G+a1*(R+B(-1))-a2*T-a3*B-a4*(M-M(-1))=0;
 
//Money Policy Rule
  M=M(-1)-d1*pi-d2*S-d3*Y;
//market clearing
  Y=A-alpha*K+(1-alpha)*H-vd;
Y=(1-omicron)*(b1*C+b2*I+b3*G+(b1+b2+b3)*eta*(P-Pd))+(re+X)/b4;
K(+1)=(1-delta)*K+b5*I;
P=(1-omicron)*Pd+omicron*Pm;
vd=xid*PId^gamma*(1-chid)*(vd(-1)+gamma*(pid-chid*pid(-1)))-gamma*(1-xid*PId^(gamma*(1-chid)))*pid1;
vm=xim*PIm^gamma*(1-chim)*(vm(-1)+theta*(pim-chim*pim(-1)))-theta*(1-xim*PIm^(theta*(1-chim)))*pim1;

//exogenous growth

  za=A-A(-1);
  zg=G-G(-1);
  zm=M-M(-1);
  zs=S-S(-1);
  zystar=Ystar-Ystar(-1);
  zpstar=Pstar-Pstar(-1);
  zrstar=Rstar-Rstar(-1);

//exogenous processes

  za= rhoa*za(-1) + epsa;
  zg= rhog*zg(-1) + epsg;
  zm= rhom*zm(-1) + epsm;
  zs= rhos*zs(-1) + epss;
  zystar= rhoystar*zystar(-1) + epsystar;
  zpstar= rhopstar*zpstar(-1) + epspstar;
  zrstar= rhorstar*zrstar(-1) + epsrstar;

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
 C=0;
 lambda=0;
 P=0;
 pi=0;
 Rstar=0;
 H=0;
 S=0;
 Q=0;
 M=0; 
 R=0;
 I=0;
MC=0;                            
 A=0;
 W=0;
 w=0;
 r=0;
 K=0;
 pid=0;
 Y=0;
 pid1=0;
 Pd=0;
 Pd1=0;
 e1=0;
 e2=0; 
 re=0;
 Pstar=0;
 Pm=0;
 Pm1=0;
 pim=0;
 pim1=0;
 f1=0;
 f2=0;
 IM=0;
 X=0;
 Ystar=0;
 Px=0;
 T=0;
 G=0;
 B=0;
 vd=0;
 vm=0;
 za=0;
 zg=0;
 zm=0;
 zs=0;
 zystar=0;
 zpstar=0;
 zrstar=0; 
 epsa=0;
 epsg=0;
 epsm=0;
 epss=0;
 epsystar=0;
 epspstar=0;
 epsrstar=0;
end;

shocks;
var epsa;              stderr  0.01 ;   
var epsg;              stderr  0.01 ; 
var epsm;              stderr  0.01 ; 
var epss;              stderr  0.01 ;     
var epsystar;          stderr  0.0046 ;    
var epspstar;          stderr  0.01 ; 
var epsrstar;          stderr  0.0195 ;      
end;

steady;

model_diagnostics;

stoch_simul(periods=2000, drop=200);