close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------
periods 20000;
var c, lambda, mu, inf, i, vg, ug, f, mg, vb, ub, mb, tetag, tetab, qg, qb,pg, pb, ng, s, nb, u, n, qr, xg, a, hg, wg, xb, hb, wb, yg, yb, y, mc, pbb, pa, pstar, pbinf, painf;
varexo epsm epsa;

parameters beta, hbt, mm, Mgg, Mbb, rho, kappag, kappah, alfah, fi, kappab, alfa, ki, eta, gama, z, epsi, phi, tavi, alfapi, alfau, rr ubar;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

beta = 0.99;      
hbt = 0.6;         
mm = 0.4;          
Mgg = 0.6;          
Mbb = 0.6;          
rho = 0.1;         
kappag = 0.16;     
kappab = 0.04;     
alfah = 0.9;       
kappah = 1.2872181;      
fi = 10.0;           
alfa = 0.04;       
ki = 1.1;          
eta = 0.5;        
gama = 0.4;        
z = 0.4;          
epsi = 11.0;         
phi = 0.67;        
tavi = 0.9;        
alfapi = 1.5;      
alfau = 0.5;       
rr = 0.923;        
ubar = 0.12;

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

model;
  mu = 1/(exp(c)-hbt*exp(c(-1))); 
  mu-beta*hbt*mu(+1) = exp(lambda);
  1 = beta*exp(i)*exp(lambda(+1)-lambda)/inf(+1);
  exp(mg) = Mgg*exp(vg)^(1-mm)*exp(ug+f)^mm;
  exp(mb) = Mbb*exp(vb)^(1-mm)*exp(ub)^mm;
  exp(tetag) = exp(vg-(ug+f));
  exp(tetab) = exp(vb-ub);
  exp(qg) = Mgg*exp(tetag)^-mm;
  exp(qb) = Mbb*exp(tetab)^-mm;
  exp(pg) = Mgg*exp(tetag)^(1-mm);
  exp(pb) = Mbb*exp(tetab)^(1-mm);
  exp(ng) = (1-rho)*(exp(ng(-1))+exp(mg(-1)));
  exp(nb) = (1-rho)*(exp(nb(-1))+exp(mb(-1))-exp(s(-1))*exp(pg(-1))*exp(nb(-1)));
  exp(u) = exp(ug)+exp(ub);
  exp(u) = 1-exp(ng)-exp(nb);
  exp(n) = exp(ng)+exp(nb);
  exp(qr) = exp(s+nb+pg-n);
  exp(-qg) = (1/kappag)*(1-rho)*beta*exp(lambda(+1)-lambda)*(exp(xg(+1)+a(+1))*exp(hg(+1))^alfah-exp(wg(+1)+hg(+1))+kappag*exp(-qg(+1)));
  exp(-qb) = (1/kappab)*(1-rho)*beta*exp(lambda(+1)-lambda)*(exp(xb(+1)+a(+1))*exp(hb(+1))^alfah-exp(wb(+1)+hb(+1))+(1-exp(s(+1)+pg(+1)))*kappab*exp(-qb(+1)));
  alfa*ki*exp(s)^(ki-1) = (eta/(1-eta))*exp(pg)*(kappag*exp(-qg)-kappab*exp(-qb));
  exp(f) = exp(s+nb);
  kappag*exp(tetag) = kappab*exp(tetab); 
  exp(yg) = exp(a+ng)*exp(hg)^alfah;
  exp(yb) = exp(a+nb)*exp(hb)^alfah;
  exp(y) = exp(yg)^gama*exp(yb)^(1-gama);
  exp(mc) = (1/gama)^gama*(1/(1-gama))^(1-gama)*exp(xg)^gama*exp(xb)^(1-gama);
  exp(y) = exp(c)+kappag*exp(vg)+kappab*exp(vb);
  exp(wg) = eta*(exp(xg)*exp(a)*exp(hg)^(alfah-1)+kappag*exp(tetag)/exp(hg))+(1-eta)*((kappah/(1+fi))*exp(hg)^fi/exp(lambda)+(z/exp(hg)));
  exp(wb) = eta*(exp(xb+a)*exp(hb)^(alfah-1)+(1-exp(s))*kappag*exp(tetag-hb))+(1-eta)*((z/exp(hb))+(kappah/(1+fi))*exp(hb)^fi/exp(lambda)+alfa*exp(s)^ki/exp(hb));
  alfah*exp(xg+a+lambda) = ksi*exp(hg)^(1+fi-alfah);
  alfah*exp(xb+a+lambda) = ksi*exp(hb)^(1+fi-alfah);
  exp(yg-yb) = ((1-gama)/gama)*exp(xb-xg);
  exp(pstar) = (epsi/(epsi-1))*(pbb/pa);
  pbinf = pbb*inf^epsi;
  pbb = exp(lambda+mc+y)+phi*beta*pbinf(+1);
  painf = pa*inf^(epsi-1);
  pa = exp(lambda+y)+phi*beta*painf(+1);
  1 = phi/inf^(1-epsi)+(1-phi)*exp(pn)^(1-epsi);
  i = -(1-tavi)*log(beta)+tavi*i(-1)+alfapi*(1-tavi)*log(inf(+1))+alfau*(1-tavi)*(u-log(ubar))+epsm;
  a = rr*a(-1)+epsa;
end;

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

initval;
  s = log(0.36);
  tetag = log(alfa)+log(ki)+log(1-eta)-log(eta)-log(kappag-kappab*(kappag/kappab)^mm)+(ki-1)*s;
  pg = log(Mgg)+(1-mm)*tetag;
  qg = log(Mgg)-mm*tetag;
  tetab = log(kappag)-log(kappab)+tetag;
  pb = log(Mbb)+(1-mm)*tetab;
  qb = log(Mbb)-mm*tetab;
  nb = log(rho*(1-ubar)-(1-rho)*exp(pg)*ubar)-log(1-(kappab/kappag)^(1-mm))-log(rho+(1-rho)*exp(s)*exp(pg));
  ng = log(1-exp(ng)-ubar);
  mg = log(rho)+ng-log(1-rho);
  vg = -log(Mgg)+mg+mm*tetag;
  mb = -log(1-rho)+log(rho+(1-rho)*exp(s)*exp(pg))+nb;
  vb = -log(Mbb)+mb+mm*tetab;
  f = s+nb;
  ub = -log(Mbb)+mb+(mm-1)*tetab;
  ug = log(ubar-exp(ub));
  u = log(0.12);
  n = log(0.88);
  qr = s+pg+nb-n;
  hg = 0;
  hb = -0.0029733;
  xg = log(kappah)-log(alfah)+((1+fi-alfah)*hg)-lambda-a;
  xb = log(kappah)-log(alfah)+(1+fi-alfah)*hb-lambda-a;
  wg = log(eta*(exp(xg)*exp(a)*exp(hg)^(alfah-1)+kappag*exp(tetag-hg))+(1-eta)*(kappah*exp(hg)^fi/((1+fi)*exp(lambda))+(z*exp(-hg))));
  wb = log(eta*(exp(xb)*exp(a)*exp(hb)^(alfah-1)+(1-exp(s))*exp(tetag)*kappag*exp(-hb))+(1-eta)*(z*exp(-hb)+(kappah*exp(hb)^fi)/((1+fi)*exp(lambda))+alfa*exp(s)^ki*exp(-hb)));
  mc = gama*(xg-log(gama))+(1-gama)*(xb-log(1-gama));
  yg = a+ng+alfah*hg;
  yb = a+nb+alfah*hb;
  y = gama*yg+(1-gama)*yb;
  i = -log(beta);
  inf = 1;
  pbb = exp(lambda+y)/(1-phi*beta);
  pbinf = pbb;
  pa = exp(lambda+y)/(1-phi*beta);
  painf = pa;
  pstar = 0;
  c = log(exp(y)-kappag*exp(vg)-kappab*exp(vb));
  mu = -log((1-hbt)*exp(c));
  lambda = log(1-beta*hbt)-log((1-hbt)*exp(c));
  epsm = 0;
  epsa = 0;
  a = 0;
end;







shocks;
var epsa = 0.000049;
var epsa, epsm = 0;
var epsm = 0.000006;
end;


stoch_simul(IRF=24, order = 1) y n u mc inf vg vb ;