%endogenous variables
var yt kt nt it gt nut vt xt qt mt ut kpt e_zt e_it Gammat pit lambdat ct e_bt
rkt qkt pwt at wt et mut chit wot st bt e_wt e_etat e_pt rt fyt e_rt e_gt thetat
% observables

dy dinve dc dw labobs pinfobs robs;

% exogenous variables

varexo innovation_bt innovation_it innovation_zt 

innovation_etat innovation_pt innovation_rt

innovation_gt ;

parameters alpha beta delta  rho yg s csi gammaz psinu etak h bb eta lambda lambdap gamma epsilonp 
epsilong rpi ry rhos rhob rhoi  rhoz rhoeta rhop rhor rhog sigmab sigmai sigmaz sigma
sigmap sigmar sigmag;

// fixed parameters
alpha = 0.33;
beta = 0.99;
delta = 0.025;
sigma = 0.5;
rho = 0.895;
yg= 0.2;
s = 0.95;
csi = 10;
%-----------------------estimated parameters prior
psinu = 0.5;
etak = 4;
h = 0.5;
bb = 0.5;
eta = 0.5;
lambda = 0.75;
lambdap = 0.66;
gamma = 0.5;
epsilonp = 1.15;
%to make it consistent with 0.2 ss gov't share
epsilong = 1.25;
rpi = 1.7;
ry = 0.125;
rhos = 0.75;
gammaz = 1.025;
rhob = 0.5;
rhoi = 0.5;
rhoz = 0.5;
rhoeta = 0.5;
rhop = 0.5;
rhor = 0.5;
rhog = 0.5;
sigmab = 0.15;
sigmai = 0.15;
sigmaz = 0.15;
sigmaeta = 0.15;
sigmap = 0.15;
sigmar = 0.15;
sigmag = 0.15;


model (linear);
#rk = (gammaz/beta)-1+delta; 
#x = 1-rho;
#n=s/(x+s);
#u = 1-n;
#qk = 1;
#k = ((rk/alpha)^(-1/(1-alpha)))*n;
#a = (1-alpha)*(k/n)^alpha;
#y = n*(k/n)^alpha;
#i = y*(1-(1-delta)/gammaz)*gammaz*(k/n)^(1-alpha);
#g = yg*y;
#mu = 1/(1-lambda*beta);
#e = 1/(1-rho*lambda*beta);
#chi = eta/(eta+(1-eta)*mu/e);
#yi = i/y;
#ynu = rk*k/y;
#htilda = h/gammaz;
#etanu = (1-psinu)/psinu; 
#betatilda = beta/gammaz;
#pw=1; %to check this
#kappa=(pw*a-chi*a-(1-lambda)*bb*pw*a)/((x+chi*beta*(x^2)/2)+chi*beta*s*x+(1-chi)*bb
*(beta/gammaz)*(1/2)*(x^2)+beta*x^2*(1/2)-beta*rho*x);
#b=bb*(pw*a+(beta/gammaz)*(kappa/2)*x^2);
#w=chi*(a+beta*(kappa/2)*(x^2)+beta*kappa*s*x)+(1-chi)*b;
#c = (1-(g/y+i/y+kappa/2*x^2*n/y))*y;
#Gamma = 1; %%to check this
#yc = c/y;
#yx = (kappa/2)*(x^2*n/y);
#kappaa = (kappa*x)^(-1)*pw*a;

#kappaw = (kappa*x)^(-1)*w;

#kappalambda = beta*(1+rho)/2;
#phia = chi*pw*a*w^(-1);




#H = chi/(1-chi)*kappa*x; % i take the initiative for finding H_upperbar, from equations  41 and B3 of the paper
#phis = (1-chi)*s*beta*H*w^(-1);
#phix = chi*beta*kappa*x^2*w^(-1);
#phib = (1-chi)*b*w^(-1);
#phichi = chi*(1-chi)^(-1)*kappa*x*w^(-1);
#phieta = phichi*(1-chi)*(1-eta)^(-1);

#psi = chi*beta*lambda*mu+(1-chi)*rho*beta*lambda*e; % to check

#tau = psi/(1+psi); %to check
#taup = 1+(epsilonp-1)*csi;
#smallcp = (1-lambdap)*(1-lambdap*beta)*(lambdap)^(-1);
#pi = 1;
#gammap = pi;
#fip = 1+beta*gammap;
#iotaf = beta*(fip)^(-1);
#iotao = (smallcp/taup)*(fip)^(-1);
#iotab = gammap*(fip)^(-1);
#G = (1-eta*x*beta*lambda*mu)*eta^(-1)*mu*kappaw;
#tau1 = (kappaw*mu*phix+phichi*(1-chi)*(x*beta*lambda)*(kappaw*mu)*mu*(rho*beta)+phis*G)*(1-tau);
#tau2 = -(kappaw*mu)*phichi*(1-chi)*(x*beta*lambda)*mu*(1-tau);
#smallc = (1-lambda)*(1-tau)*lambda^(-1);
#fi = (1+tau2)+smallc+(tau*lambda^(-1)-tau1);
#gammab = (1+tau2)*fi^(-1);
#gammao = smallc*fi^(-1);
#gammaf = (tau*lambda^(-1)-tau1)*fi^(-1);


%technology (1)
yt=alpha*kt+(1-alpha)*nt;
% resource constraint (2)
yt=yc*ct+yi*it+yg*gt+ynu*nut+yx*(2*xt+nt(-1));
% matching(3)
mt=sigma*ut+(1-sigma)*vt;
% employment dynamics (4)
nt=nt(-1)+(1-rho)*xt;
% transition probabilities (5)(6)
qt=mt-vt;
st=mt-ut;
% unemployment (7)
ut=-(n/u)*nt(-1);
% effective capital (8)
kt+e_zt=nut+kpt(-1);
% physical capital dynamics (9)
kpt=((1-delta)/gammaz)*(kpt(-1)-e_zt)+(1-(1-delta)/gammaz)*(it+e_it);
% aggregate vacancies (10)
xt=qt+vt-nt(-1);

% consumption-saving (11)
0=Gammat(+1)+(rt-pit(+1))-e_zt(+1);
% marginal utility (12)
(1-htilda)*(1-beta*htilda)*lambdat=htilda*(ct(-1)-e_zt)-
(1+beta*htilda^2)*ct+beta*htilda*(ct(+1)+e_zt(+1))+(1-htilda)*(e_bt-beta*htilda*e_bt(+1));
Gammat=lambdat/lambdat(-1);
% capital utilization (13)
nut=etanu*rkt;
% investment (14)
it=(1/(1+beta))*(it(-1)-e_zt)+(1/(etak*gammaz^2))*(1/(1+beta))*(qkt+e_it)
+(beta/(1+beta))*(it(+1)+e_zt(+1));
% capital renting (15)
pwt+yt-kt=rkt;

% tobin's q (16)    
qkt=betatilda*(1-delta)*qkt(+1)+(1-betatilda*(1-delta))*rkt(+1)-(rt-pit(+1));

% aggregate hiring rate  (17)
xt=kappaa*(pwt+at)-kappaw*wt+kappalambda*Gammat(+1)+beta*xt(+1);

% marginal product of labor (18)
at=yt-nt;


% weight in Nash bargaining (19) (20) (21)
chit=-(1-chi)*(mut-et);
et=(rho*lambda*beta)*(Gammat(+1)-pit(+1)+gamma*pit+et(+1)-e_zt(+1));
mut=(x*lambda*beta)*xt(+1)-(x*lambda*beta)*(kappaw*mu)*mu*
(wt+gamma*pit-pit(+1)-e_zt(+1)-wt(+1))+
(lambda*beta)*(mut(+1)+Gammat(+1)+gamma*pit-pit(+1)-e_zt(+1));
% spillover free target wage (22)
wot=phia*(pwt+at)+(phis+phix)*xt(+1)+phis*st(+1)+phib*bt+
(phis+phix/2)*Gammat(+1)+phix*(chit-(rho-s)*beta*chit(+1))+e_wt;
e_wt = phieta*(1-(rho-s)*beta*rhoeta)*e_etat;

% aggregate wage (23)
wt=gammab*(wt(-1)-pit+gamma*pit(-1)-e_zt)
+gammao*wot+gammaf*(wt(+1)+pit(+1)-gamma*pit+e_zt(+1));

% Philips curve (24)
pit=iotab*pit(-1)+iotao*(pwt+e_pt)+iotaf*pit(+1);

% Taylor rule (25)
rt=rhos*rt(-1)+(1-rhos)*(rpi*pit+ry*(yt-fyt))+e_rt;

% government spending (26)
gt=yt+(1-yg)/yg*e_gt;

% market tightness (27)
thetat=vt-ut;

% benefits (28)
bt=kpt;

% FLEX ECONOMY   30
fyt = 0;    

%------------------- shocks

% preferences
e_bt=rhob*e_bt(-1)+innovation_bt;

% investment 
e_it=rhoi*e_it(-1)+innovation_it;

% labour productivity    33
e_zt=rhoz*e_zt(-1)+innovation_zt;

% bargaining power
e_etat=rhoeta*e_etat(-1)+innovation_etat;

% markup   35
e_pt=rhop*e_pt(-1)+innovation_pt;

% monetary policy
e_rt=rhor*e_rt(-1)+innovation_rt;

% govenrment spendings     37
e_gt=rhog*e_gt(-1)+innovation_gt;
% ------------ observation equations


dy=100*(gammaz-1)+yt-yt(-1)+e_zt;

dc=100*(gammaz-1)+ct-ct(-1)+e_zt;

dinve=100*(gammaz-1)+it-it(-1)+e_zt;

dw=100*(gammaz-1)+wt-wt(-1)+e_zt;

labobs=nt;

pinfobs=pit;

robs=rt;


end;
varobs  dy dc dinve dw labobs pinfobs robs;
%%%%%%%%%%%%%%%%%%%%%%%%%%% ESTIMATION
 
estimated_params;

rhoz, beta_pdf, 0.5, 0.2;
rhob, beta_pdf, 0.5, 0.2;
rhoi, beta_pdf, 0.5, 0.2;
rhop, beta_pdf, 0.5, 0.2;
rhoeta, beta_pdf, 0.5, 0.2;
rhog, beta_pdf, 0.5, 0.2;
rhor, beta_pdf, 0.5, 0.2;

stderr innovation_zt, inv_gamma_pdf, 0.15, 0.15;
stderr innovation_bt, inv_gamma_pdf, 0.15, 0.15;
stderr innovation_it, inv_gamma_pdf, 0.15, 0.15;
stderr innovation_pt, inv_gamma_pdf, 0.15, 0.15;
stderr innovation_etat, inv_gamma_pdf, 0.15, 0.15;
stderr innovation_gt, inv_gamma_pdf, 0.15, 0.15;
stderr innovation_rt, inv_gamma_pdf, 0.15, 0.15;

psinu, beta_pdf, 0.5, 0.1;
etak, normal_pdf, 4, 1.5;
h, beta_pdf, 0.5, 0.1;
eta, beta_pdf, 0.5, 0.1;
bb, beta_pdf, 0.5, 0.1;
lambda, beta_pdf, 0.75, 0.1;
lambdap, beta_pdf, 0.66, 0.1;
gamma, uniform_pdf, , , 0, 1;
epsilonp, normal_pdf, 1.15, 0.05;
rpi, normal_pdf, 1.7, 0.3;
ry, gamma_pdf, 0.125, 0.1;
rhos, beta_pdf, 0.75, 0.1;
gammaz, uniform_pdf, , , 1, 1.5;


end;


check;
resid(1);
steady;


estimation(order=1,datafile=trigari_dataUS, plot_priors=0, mh_jscale=0.3, mh_nblocks=2, mode_compute=6);