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

close all;

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

var y,c,inv,l,k,gk,gc,gi,cstar,w,rk,r,b,tau,tr,lambda,q,u,z,a,mu,uc,ul,ugi,ugc,utau,utr;
varexo ez,ea,ec,el,emu,egi,egc,etau,etr;

parameters beta,alpha,delta,delta0,delta1,delta2,kappa,h,eta,psi,
si,sgi,sgc,sb,
ybar,cbar,invbar,lbar,kbar,gkbar,gcbar,gibar,cstarbar,wbar,rkbar,rbar,bbar,taubar,trbar,lambdabar,
alpha_gc,alpha_gi,phi_gi,gamma_gi,phi_gc,gamma_gc,phi_tau,gamma_tau,phi_tr,gamma_tr,
rhoz,rhoa,rhoc,rhol,rhomu,rhogi,rhogc,rhotau,rhotr,
sigmaz,sigmaa,sigmac,sigmal,sigmamu,sigmagi,sigmagc,sigmatau,sigmatr;

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

beta=.98;
alpha=.5;
delta0=.025;
delta=delta0;
delta2=0.05;
h=0.5;
eta=2;
psi=20;
kappa=4;
sgi=0.02;
sgc=0.16;
sb=0.3;
taubar=0.15;

alpha_gc=-0.5;
alpha_gi=0.02;

phi_gi=0.05;
gamma_gi=0.05;
phi_gc=0.05;
gamma_gc=0.05;
phi_tau=0.05;
gamma_tau=0.05;
phi_tr=0.05;
gamma_tr=0.05;

rhoz=0.5;
rhoa=0.5;
rhoc=0.5;
rhol=0.5;
rhomu=0.5;
rhogi=0.5;
rhogc=0.5;
rhotau=0.5;
rhotr=0.5;
sigmaz=0.02;
sigmaa=0.02;
sigmac=0.02;
sigmal=0.02;
sigmamu=0.02;
sigmagi=0.02;
sigmagc=0.02;
sigmatau=0.02;
sigmatr=0.02;

rbar=1/beta;
rkbar=(1-beta*(1-delta0))/beta/(1-taubar);
delta1=rkbar*(1-taubar);
si=delta0*alpha/rkbar;
lbar=((1-beta*h)*(1-taubar)*(1-alpha)/psi/(1-h)/(1-si-sgi-(1-alpha_gc)*sgc))^(1/(eta+1));
ybar=(alpha/rkbar)^(alpha/(1-alpha-alpha_gi))*lbar^((1-alpha)/(1-alpha-alpha_gi))*(sgi/delta)^(alpha_gi/(1-alpha-alpha_gi));
kbar=alpha*ybar/rkbar;
gcbar=sgc*ybar;
gibar=sgi*ybar;
invbar=delta0*kbar;
gkbar=gibar/delta;
bbar=sb*ybar;
trbar=taubar*ybar-gcbar-gibar-(rbar-1)*bbar;
cbar=ybar-invbar-gcbar-gibar;
wbar=(1-alpha)*ybar/lbar;
lambdabar=psi*lbar^eta/wbar/(1-taubar);
cstarbar=cbar+alpha_gc*gcbar;


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

model(linear);
uc-(cstar-h*cstar(-1))/(1-h)-beta*h*uc(+1)+beta*h*(cstar(+1)-h*cstar)/(1-h)=(1-beta*h)*lambda;
cstarbar*cstar=cbar*c+alpha_gc*gcbar*gc;
ul+eta*l=lambda-taubar/(1-taubar)*tau+w-mu;
lambda+q=beta*rkbar*((1-taubar)*(lambda(+1)+rk(+1)+u(+1))-taubar*tau)+beta*(1-delta0)*(lambda(+1)+q(+1))-beta*delta1*u(+1);
rkbar*((1-taubar)*rk-taubar*tau)=delta1*q+delta2*u;
lambda=lambda+q+a-kappa*(inv-inv(-1))+beta*kappa*(inv(+1)-inv);
lambda=lambda(+1)+r;
k=(1-delta0)*k(-1)-delta1*u+delta0*(a+inv);
gk=(1-delta)*gk(-1)+delta*gi;
w=y-l;
rk=y-u-k(-1);
bbar*b+taubar*ybar*(tau+y)=gcbar*gc+gibar*gi+rbar*bbar*(r(-1)+b(-1))+trbar*tr;
ybar*y=cbar*c+invbar*inv+gcbar*gc+gibar*gi;
y=z+alpha*(u+k(-1))+(1-alpha)*l+alpha_gi*gk(-1);
gi=-phi_gi*y-gamma_gi*b(-1)+ugi;
gc=-phi_gc*y-gamma_gc*b(-1)+ugc;
tau=phi_tau*y+gamma_tau*b(-1)+utau;
tr=-phi_tr*y-gamma_tr*b(-1)+utr;

z=rhoz*z(-1)+sigmaz*ez;
a=rhoa*a(-1)+sigmaa*ea;
uc=rhoc*uc(-1)+sigmac*ec;
ul=rhol*ul(-1)+sigmal*el;
mu=rhomu*mu(-1)+sigmamu*emu;
ugi=rhogi*ugi(-1)+sigmagi*egi;
ugc=rhogc*ugc(-1)+sigmagc*egc;
utau=rhotau*utau(-1)+sigmatau*etau;
utr=rhotr*utr(-1)+sigmatr*etr;


end;

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

steady;

check;

shocks;
var ez=1;
var ea=1;
var ec=1;
var el=1;
var emu=1;
var egi=1;
var egc=1;
var etau=1;
var etr=1;
end;

%----------------------------------------------------------------
% 4. Estimation
%----------------------------------------------------------------

estimated_params;
h,beta_pdf,0.5,0.2;
eta,normal_pdf,2,0.75;
delta2,gamma_pdf,0.7,0.5;
kappa,gamma_pdf,4,1;

alpha_gc,uniform_pdf,,,-1,0.33;
alpha_gi,inv_gamma_pdf,0.03,0.1;

phi_gi,gamma_pdf,0.05,0.05;
gamma_gi,gamma_pdf,0.05,0.05;
phi_gc,gamma_pdf,0.05,0.05;
gamma_gc,gamma_pdf,0.05,0.05;
phi_tau,gamma_pdf,0.05,0.05;
gamma_tau,gamma_pdf,0.05,0.05;
phi_tr,gamma_pdf,0.05,0.05;
gamma_tr,gamma_pdf,0.05,0.05;

rhoz,beta_pdf,0.5,0.2;
rhoa,beta_pdf,0.5,0.2;
rhoc,beta_pdf,0.5,0.2;
rhol,beta_pdf,0.5,0.2;
rhomu,beta_pdf,0.5,0.2;
rhogi,beta_pdf,0.5,0.2;
rhogc,beta_pdf,0.5,0.2;
rhotau,beta_pdf,0.5,0.2;
rhotr,beta_pdf,0.5,0.2;
sigmaz,inv_gamma_pdf,0.01,2;
sigmaa,inv_gamma_pdf,0.01,2;
sigmac,inv_gamma_pdf,0.01,2;
sigmal,inv_gamma_pdf,0.01,2;
sigmamu,inv_gamma_pdf,0.01,2;
sigmagi,inv_gamma_pdf,0.01,2;
sigmagc,inv_gamma_pdf,0.01,2;
sigmatau,inv_gamma_pdf,0.01,2;
sigmatr,inv_gamma_pdf,0.01,2;
end;

varobs y c gc gi tau l w;

estimation(datafile=data,mh_replic=1000);



%--------------------------------------------------------------------------