
warning off;

var Cl Hl Ll laml Cb Hb b Lb lam1 lam2 Rf Rb Rhm Lam Omega n f d lam phi Qb Y YhH pH Ym theta1
    theta2 L K I pHr pm YH x1 x2 vp pFr pF YF x1r x2r vpr R pmF EX X YmF hm 
    Rr infl w qh Qf s A;

varexo eA eR eRr eX es;

parameters betab betal gammah Phi eta M kappa xi sigma zeta theta Delta gamma kappaI
           alpha delta epsilon Phip inflH inflF Pis Ys rhor rhopi rhoY v rhoX Xs 
           fracYH v1 v2 hms rhoRr Rrs Rs infls rhoA rhos ss EXs pmFs YmFs;
          

load param;

model;

%%%%%%%%%%%%%% Lender Households %%%%%%%%%%%%%%
% (1) Euler Equation for consumption 
exp(laml)=1/exp(Cl);

% (2) Euler equation for debt
exp(laml)=betal*exp(laml(+1))*exp(R)/exp(infl(+1));

% (3) Euler equation for housing
exp(laml)*exp(qh)-gammah/exp(Hl)=betal*exp(laml(+1))*exp(qh(+1));

% (4) labor supply
-Phi*exp(eta*Ll)+exp(laml)*exp(w)=0;

%%%%%%%%%%%%%%% Borrowers %%%%%%%%%%%%%%%%%%
% (5) Budget
%exp(Cb)+exp(b(-1))/exp(infl)+exp(qh)*(exp(Hb)-exp(Hb(-1)))=exp(w)*exp(Lb)+exp(Qb)*(exp(b)-kappa*exp(b(-1))/exp(infl));

% (5) Borrowing budget
exp(Qb)*(exp(b)-kappa*exp(b(-1))/exp(infl))=xi*exp(qh)*exp(Hb);

% (6) Euler equation for consumption
1/exp(Cb)-exp(lam1)=0;

% (7) Euler equation for housing
gammah/exp(Hb)-exp(lam1)*exp(qh)+betab*exp(lam1(+1))*exp(qh(+1))+xi*exp(lam2)*exp(qh)=0;

% (8) lam1 and lam2
(exp(lam1)-exp(lam2))*exp(Qb)+betab/exp(infl(+1))*(kappa*exp(lam2(+1))*exp(Qb(+1))-exp(lam1(+1))*(1+kappa*exp(Qb(+1))))=0;

% (9) Labor supply
-Phi*exp(eta*Lb)+exp(lam1)*exp(w)=0;

%%%%%%%%%%%%%%%%%%%%% Financial Market %%%%%%%%%%%%%%%%%%%%%
% (10) Rf
exp(Rf)=(1+kappa*exp(Qf))/exp(Qf(-1));

% (11) Rb
exp(Rb)=(1+kappa*exp(Qb))/exp(Qb(-1));

% (12) Budget
exp(Qb)*exp(b)+exp(Qf)*exp(f)=exp(d)+exp(n)+exp(s)*exp(hm);

% (13) n
exp(n)=sigma/exp(infl)*((exp(Rf)-exp(R(-1)))*exp(Qf(-1))*exp(f(-1))+(exp(Rb)-exp(R(-1)))*exp(Qb(-1))*exp(b(-1))+exp(R(-1))*exp(n(-1))-exp(Rhm(-1))*exp(s(-1))*exp(hm(-1)))+M;

% (14) Rhm
exp(Rhm)=exp(R)*exp(-zeta/hms*exp(s)*exp(hm));

% (15) Lam
exp(Lam)=betal*exp(laml)/exp(laml(-1));

% (16) Omega
exp(Omega)=1-sigma+sigma*theta*exp(phi);

% (17) 
exp(Lam(+1))*exp(Omega(+1))/exp(infl(+1))*(exp(Rf(+1))-exp(R))=exp(lam)/(1+exp(lam))*theta;

% (18)
exp(Lam(+1))*exp(Omega(+1))/exp(infl(+1))*(exp(Rb(+1))-exp(R))=exp(lam)/(1+exp(lam))*theta*Delta;

% (19)
exp(Qf)*exp(f)+Delta*exp(Qb)*exp(b)=exp(phi)*exp(n);

% (20)
%theta*exp(phi)=exp(Lam(+1))*exp(Omega(+1))/exp(infl(+1))*(exp(R)*exp(n)-exp(Rhm(-1))*exp(s(-1))*exp(hm(-1)));
exp(phi)=(exp(R)*exp(Lam(+1))*exp(Omega(+1))/exp(infl(+1)))  /  (theta-(exp(Rf(+1))-exp(R))*exp(Lam(+1))*exp(Omega(+1))/exp(infl(+1)));

%%%%%%%%%%%%%%%%% Final Product Producer %%%%%%%%%%%%%%%
% (21) Production function
exp(Y)=exp(YhH*gamma)*exp((1-gamma)*YF);

% (22)
exp(YhH)/exp(YF)=gamma/(1-gamma);

% (23) Domestic tradable goods
exp(YH)=exp(YhH)+exp(EX);

% (24) Final goods price
1=(exp(pH)/gamma)^gamma*(exp(pF)/(1-gamma))^(1-gamma);

%%%%%%%%%%%%%%% Intermediate goods%%%%%%%%%%%%%%
% (25) Production function
exp(Ym)=exp(A)*exp(alpha*(K(-1)))*exp((1-alpha)*L);

% (26) Investment budget
exp(I)=exp(Qf)*(exp(f)-kappa*exp(f(-1))/exp(infl));

% (27) Capital evolvment
exp(K)=(1-delta)*exp(K(-1))+exp(I)*(1-kappaI/2*(exp(I)/exp(I(-1))-1)^2);

% (28) FOC for L
(1-alpha)*exp(pm)*exp(A)*exp(alpha*K(-1))*exp(-alpha*L)=exp(w);

% (29) FOC for K
exp(Lam(+1))*(alpha*exp(pm(+1))*exp(A(+1))*exp((alpha-1)*K)*exp((1-alpha)*L(+1))+exp(theta2(+1))*(1-delta))=exp(theta2);

% (30) FOC for I
-(1+exp(theta1))+exp(theta2)*(1-kappaI/2*(exp(I)/exp(I(-1))-1)^2-kappaI*(exp(I)/exp(I(-1))-1)*(exp(I)/exp(I(-1))))+kappaI*exp(Lam(+1))*exp(theta2(+1))*(exp(I(+1))/exp(I)-1)*(exp(I(+1))/exp(I))^2=0;

% (31) FOC for f
exp(Qf)*(1+exp(theta1))-exp(Lam(+1))/exp(infl(+1))*(1+kappa*exp(Qf(+1))*(1+exp(theta1)))=0;

%%%%%%%%%%%%%%%%%%%%%%%%% Retailer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (32) x1
exp(x1)=(exp(pHr)/exp(pH))^(1-epsilon)*exp(YH)+Phip*inflH^(1-epsilon)*exp(Lam(+1))*(exp(pHr)/exp(pH(+1)))^(1-epsilon)*exp(x1(+1));

% (33) x2
exp(x2)=(exp(pHr)/exp(pH))^(-epsilon)*exp(pm)/exp(pH)*exp(YH)+Phip*inflH^(-epsilon)*exp(Lam(+1))*(exp(pHr)/exp(pH(+1)))^(-epsilon)*exp(x2(+1));

% (34) x1/x2
exp(x1)/exp(x2)=epsilon/(epsilon-1);

% (35) pHstar/pH
exp((1-epsilon)*pH)=(1-Phip)*exp((1-epsilon)*pHr)+Phip*inflH^(1-epsilon)*exp((1-epsilon)*pH(-1));

% (36) Ym/YH
exp(Ym)=exp(vp)*exp(YH);

% (37) vp
exp(vp)=(1-Phip)*(exp(pHr)/exp(pH))^(-epsilon)+Phip*(exp(pH(-1))*inflH/exp(pH))^(-epsilon)*exp(vp(-1));

%%%%%%%%%%%%%%%%%%%%%% Retailers for foreign tradable %%%%%%%%%%%%%%%%%
% (38) x1r
exp(x1r)=(exp(pFr)/exp(pF))^(1-epsilon)*exp(YF)+Phip*inflF^(1-epsilon)*exp(Lam(+1))*(exp(pFr)/exp(pF(+1)))^(1-epsilon)*exp(x1r(+1));

% (39)x2r
exp(x2r)=(exp(pFr)/exp(pF))^(-epsilon)*exp(pmF)/exp(pF)*exp(YF)+
Phip*inflF^(-epsilon)*exp(Lam(+1))*(exp(pFr)/exp(pF(+1)))^(-epsilon)*exp(x2r);

% (40) x1r/x2r
exp(x1r)/exp(x2r)=epsilon/(epsilon-1);

% (41) pFstar/pF
exp((1-epsilon)*pF)=Phip*(exp(pF(-1))*inflF)^(1-epsilon)+(1-Phip)*exp((1-epsilon)*pFr);

% (42) Ym/YF
exp(YmF)=exp(vpr)*exp(YF);

% (43) vpr
exp(vpr)=(1-Phip)*(exp(pFr)/exp(pF))^(-epsilon)+Phip*(exp(pF(-1))*inflF/exp(pF))^(-epsilon)*exp(vpr(-1));

%%%%%%%%%%%%%%%% Monetary Policy %%%%%%%%%%%%%%%%%
% (44) R
R=rhor*R(-1)+(1-rhor)*(log(Rs)+rhopi*(infl-log(Pis))+rhoY*(Y-log(Ys)))+eR;

%%%%%%%%%%%%%%%%%% Current account %%%%%%%%%%%%%%%%%%%
% (45) PmF
exp(pmF)=exp(s);

% (46) EX
exp(EX)=(exp(pH)/exp(s))^(-v)*exp(X);

% (47) X
X=(1-rhoX)*log(Xs)+rhoX*X(-1)+eX;

%%%%%%%%%%%%%%%%% Capital account%%%%%%%%%%%%%%%%%
% (48) hm
exp(hm)=hms*(exp(s)*exp(R)/(exp(s(+1))*exp(Rr)))^v1*(exp(Y(+1))/exp(Y))^v2;

% (49) Rr
Rr=(1-rhoRr)*log(Rrs)+rhoRr*Rr(-1)+eRr;

%%%%%%%%%%%%%%%%%% Clearing conditions %%%%%%%%%%%%%%%%%%%%%%
% (50) Final goods clearing
exp(Y)=exp(Cl)+exp(Cb)+exp(I);

% (51) Housing clearing
exp(Hl)+exp(Hb)=1;

% (52) Labor market clearing
exp(Ll)+exp(Lb)=exp(L);

% (53) International budget
exp(EX)-exp(pmF)*exp(YmF)-exp(Rhm(-1))*exp(s)*exp(hm(-1))=-exp(s)*(exp(hm)-exp(hm(-1)));

% (54) Exchange rate policy
ln(exp(s)/ss)=rhos*ln(  (exp(EX)-exp(pmF)*exp(YmF))  /   (EXs-pmFs*YmFs))+es;

% (55) A
A=rhoA*A(-1)+eA;


end;

%options_.noprint=1;

steady;

shocks;
var eA=1;
var eR=1;
var eRr=1;
var eX=1;
var es=1;
end;

stoch_simul(order=1,irf=50,nograph,ar=0,qz_zero_threshold=1e-8) qh Y;