%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Author:Rui Dang
% %DSGE Modelling
% %Data Generating Process 
% %Bayesian Estimation
% 1) Putting the model of Justiniano, Primiceri and Tambalotti (2009) in DSGE form
% 2) Solving for steady state
% 3) Solve the system of equations using Bayesian Inference(Metropolis-Hastings algorithm)and simulate Bayesian IRF
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set_dynare_seed(123)

%----------------------------------------------------------------
% 1. Defining variables and parameters
% 1.1 Variables
%y              % output
%k              % capital
%L              % hours
%Rk             % return on capital
%w              % real wages
%p              % inflation
%s              % marginal cost
%lambda         % multiplier
%c              % consumption
%R              % interst rate
%caplu          % capital utilization
%phi            % multiplier
%inv            % investment
%kbar           % physical capital
%wgap           % wage gap (w - mrs)
%gdp            % gdp
%ystar
%kstar
%Lstar
%Rkstar
%wstar
%mcstar
%lambdastar
%cstar
%iRstar
%caplustar
%phistar
%invstar
%kbarstar 
%wgapstar
%gdpstar
%----------------------------------------------------------------
% 1.2 shocks   
       
%ez             %productivity shock
%miu            %investment shock
%eg             %government spending shock  
%lambdap        %price mark up shock
%lambdaw        %wage mark up shock 
%eb             %intertemporal preference shock
%emp            %monetary policy shock
%----------------------------------------------------------------
% 1.3 Parameters
% alpha         %share of capital in the production function
% delta         %capital depreciation rate 
% indexp         %price indexation
% indexw          %wage indexation
% gamma100      %steady state growth rate of tech
% h             %habit formation
% lambdapss
% lambdawss
% Lss
% Pss100
% Fbeta
% gss
% niu
% xip           % price stickiness
% xiw           % wage stickiness
% chi           % elasticity of the capital utilization cost 
% fp            % reaction to inflation
% fy            %reaction to output gap 
% fdy           % reaction to output gap growth
% rhor          
% rhoz
% rhog
% rhomiu
% rholambda_p
% rholambda_w
% rhoemp
% ctrend
% conste_pinf
% conste_beta
% conste_L
% conste_R
%----------------------------------------------------------------
var y k L Rk w inf mc lamda c iR phi inv caplu kbar wgap gdp dy dc dinv dw pinf_obs r_obs L_obs ystar kstar Lstar Rkstar wstar mcstar lambdastar  cstar iRstar 
caplustar phistar  invstar  kbarstar wgapstar gdpstar;

varexo ez miu eg lambdap lambdaw eb emp;

parameters curvw cgy curvp constelab constepinf constebeta cmaw cmap calfa czcap csadjcost ctou csigma chabb ccs cinvs cfc 
cindw cprobw cindp cprobp csigl clandaw crdpi crpi crdy cry crr crhoa crhoas crhob crhog crhols crhoqs crhoms crhopinf crhow  ctrend cg ;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------
%fixed pamameters
beta    = 0.99;
delta   = 0.025;
rhoa    = 0.979;
rhog    = 0.9;
alpha   = 0.667;
eta     = 1; 
gy      = 0.2;
n_ss    = 0.3;

%estimated parameters initialisation


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

% Computation of the steady state
#kss=kLss*expLss;
#iss=(1-(1-delta)*exp(-gamma))*kss*exp(gamma);
#Rkss=exp(gamma)/beta-1+delta;
#sss=1/(1+lambdapss);
#wss=(sss*((1-alpha)^(1-alpha))/((alpha^(-alpha))*Rkss^alpha))^(1/(1-alpha));


#gamma=gamma100/100;
#beta=100/(Fbeta+100);
#rss=exp(gamma)/beta-1;
#rss100=rss*100;
#pss=pss100/100;
#gss=1/(1-gss);

#expLss=exp(Lss);
#kLss=(wss/Rkss)*alpha/(1-alpha);
#FLss=(kLss^alpha-Rkss*kLss-wss);
#yLss=kLss^alpha-FLss;
#F=FLss*expLss;

#yss=yLss*expLss;
#css=yss/gss-iss;






Model;

%flexible economy

0*(1-calfa)*a + 1*a =  calfa*rkf+(1-calfa)*(wf)  ;
	 
zcapf =  (1/(czcap/(1-czcap)))* rkf  ;
	      
rkf =  (wf)+labf-kf ;
	      
kf =  kpf(-1)+zcapf ;
	      
invef = (1/(1+cbetabar*cgamma))* (  invef(-1) + cbetabar*cgamma*invef(1)+(1/(cgamma^2*csadjcost))*pkf ) +qs ;
          
pkf = -rrf-0*b+(1/((1-chabb/cgamma)/(csigma*(1+chabb/cgamma))))*b +(crk/(crk+(1-ctou)))*rkf(1) +  ((1-ctou)/(crk+(1-ctou)))*pkf(1) ;
	      
cf = (chabb/cgamma)/(1+chabb/cgamma)*cf(-1) + (1/(1+chabb/cgamma))*cf(+1) +((csigma-1)*cwhlc/(csigma*(1+chabb/cgamma)))*(labf-labf(+1)) - (1-chabb/cgamma)/(csigma*(1+chabb/cgamma))*(rrf+0*b) + b ;
	      
yf = ccy*cf+ciy*invef+g  +  crkky*zcapf ;
	      
yf = cfc*( calfa*kf+(1-calfa)*labf +a );
	      
wf = csigl*labf 	+(1/(1-chabb/cgamma))*cf - (chabb/cgamma)/(1-chabb/cgamma)*cf(-1) ;
	      
kpf =  (1-cikbar)*kpf(-1)+(cikbar)*invef + (cikbar)*(cgamma^2*csadjcost)*qs ;








%sticky price - wage economy

y=
k=w+L-k




mc =  calfa*rk+(1-calfa)*(w) - 1*a - 0*(1-calfa)*a ;

zcap =  (1/(czcap/(1-czcap)))* rk ;

rk =  w+lab-k  ;

k =  kp(-1)+zcap ;

inve = (1/(1+cbetabar*cgamma))* (  inve(-1) + cbetabar*cgamma*inve(1)+(1/(cgamma^2*csadjcost))*pk ) +qs ;
          
pk = -r+pinf(1)-0*b +(1/((1-chabb/cgamma)/(csigma*(1+chabb/cgamma))))*b + (crk/(crk+(1-ctou)))*rk(1) +  ((1-ctou)/(crk+(1-ctou)))*pk(1) ;
	      
c = (chabb/cgamma)/(1+chabb/cgamma)*c(-1) + (1/(1+chabb/cgamma))*c(+1) +((csigma-1)*cwhlc/(csigma*(1+chabb/cgamma)))*(lab-lab(+1)) - (1-chabb/cgamma)/(csigma*(1+chabb/cgamma))*(r-pinf(+1) + 0*b) +b ;
	      
y = ccy*c+ciy*inve+g  +  1*crkky*zcap ;
	      
y = cfc*( calfa*k+(1-calfa)*lab +a );
	      
pinf =  (1/(1+cbetabar*cgamma*cindp)) * ( cbetabar*cgamma*pinf(1) +cindp*pinf(-1) 
               +((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc-1)*curvp+1)*(mc)  )  + spinf ; 
	      
w =  (1/(1+cbetabar*cgamma))*w(-1)
               +(cbetabar*cgamma/(1+cbetabar*cgamma))*w(1)
               +(cindw/(1+cbetabar*cgamma))*pinf(-1)
               -(1+cbetabar*cgamma*cindw)/(1+cbetabar*cgamma)*pinf
               +(cbetabar*cgamma)/(1+cbetabar*cgamma)*pinf(1)
               +(1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))*
               (csigl*lab + (1/(1-chabb/cgamma))*c - ((chabb/cgamma)/(1-chabb/cgamma))*c(-1) -w) 
               + 1*sw ;
 r =  crpi*(1-crr)*pinf +cry*(1-crr)*(y-yf)     
     +crdy*(y-yf-y(-1)+yf(-1))
               +crr*r(-1)
               +ms  ;
               
 log(a) = crhoa*log(a(-1)) + ea;
	      
 log(b) = crhob*log(b(-1)) + eb;
	      
 log(g) = crhog*log(g(-1)) + eg + cgy*ea;
	      
 qs = crhoqs*qs(-1) + eqs;
	      
 ms = crhoms*ms(-1) + em;
	      
 spinf = crhopinf*spinf(-1) + epinfma - cmap*epinfma(-1);
	          
 epinfma=epinf;
	      
 sw = crhow*sw(-1) + ewma - cmaw*ewma(-1) ;
	          
 ewma=ew; 
	      
 kp =  (1-cikbar)*kp(-1)+cikbar*inve + cikbar*cgamma^2*csadjcost*qs ;

 
% measurment equations 
dy=y-y(-1)+ctrend;
dc=c-c(-1)+ctrend;
dinv=inv-inv(-1)+ctrend;
dw=w-w(-1)+ctrend;

pinfobs = 1*(pinf) + constepinf;
robs =    1*(r) + conster;
L = L + constelab; 



end;

%----------------------------------------------------------------
% 4. Shocks
%----------------------------------------------------------------

shocks;
var ez;          % Technology shock
stderr 0.4618;  

var miu;            %  Investment Shock 
stderr 1.8513;

var eg;            % shock for government spending process
stderr 0.6090;

var eb;           %intertemporal preference shock
stderr 0.6017;

var emp;            % Monetary Policy Shock
stderr 0.2397;

var lamdap;         % Price Mark-Up Shock
stderr 0.1455;

var lamdaw;            % Wage Mark-Up Shock
stderr 0.2089;
end;


estimated_params;
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF
stderr ea,0.4618,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr eb,0.1818513,0.025,5,INV_GAMMA_PDF,0.1,2;
stderr eg,0.6090,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr eqs,0.46017,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr em,0.2397,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr epinf,0.1455,0.01,3,INV_GAMMA_PDF,0.1,2;
stderr ew,0.2089,0.01,3,INV_GAMMA_PDF,0.1,2;
crhoa,.9676 ,.01,.9999,BETA_PDF,0.5,0.20;
crhob,.2703,.01,.9999,BETA_PDF,0.5,0.20;
crhog,.9930,.01,.9999,BETA_PDF,0.5,0.20;
crhoqs,.5724,.01,.9999,BETA_PDF,0.5,0.20;
crhoms,.3,.01,.9999,BETA_PDF,0.5,0.20;
crhopinf,.8692,.01,.9999,BETA_PDF,0.5,0.20;
crhow,.9546,.001,.9999,BETA_PDF,0.5,0.20;
cmap,.7652,0.01,.9999,BETA_PDF,0.5,0.2;
cmaw,.8936,0.01,.9999,BETA_PDF,0.5,0.2;
csadjcost,6.3325,2,15,NORMAL_PDF,4,1.5;
csigma,1.2312,0.25,3,NORMAL_PDF,1.50,0.375;
chabb,0.7205,0.001,0.99,BETA_PDF,0.7,0.1;
cprobw,0.7937,0.3,0.95,BETA_PDF,0.5,0.1;
csigl,2.8401,0.25,10,NORMAL_PDF,2,0.75;
cprobp,0.7813,0.5,0.95,BETA_PDF,0.5,0.10;
cindw,0.4425,0.01,0.99,BETA_PDF,0.5,0.15;
cindp,0.3291,0.01,0.99,BETA_PDF,0.5,0.15;
czcap,0.2648,0.01,1,BETA_PDF,0.5,0.15;
cfc,1.4672,1.0,3,NORMAL_PDF,1.25,0.125;
crpi,1.7985,1.0,3,NORMAL_PDF,1.5,0.25;
crr,0.8258,0.5,0.975,BETA_PDF,0.75,0.10;
cry,0.0893,0.001,0.5,NORMAL_PDF,0.125,0.05;
crdy,0.2239,0.001,0.5,NORMAL_PDF,0.125,0.05;
constepinf,0.7,0.1,2.0,GAMMA_PDF,0.625,0.1;//20;
constebeta,0.7420,0.01,2.0,GAMMA_PDF,0.25,0.1;//0.20;
constelab,1.2918,-10.0,10.0,NORMAL_PDF,0.0,2.0;
ctrend,0.3982,0.1,0.8,NORMAL_PDF,0.4,0.10;
cgy,0.05,0.01,2.0,NORMAL_PDF,0.5,0.25;
calfa,0.24,0.01,1.0,NORMAL_PDF,0.3,0.05;
end;

varobs dy dc dinv L_obs pinf_obs dw r_obs;

estimation(optim=('MaxIter',200),datafile=usmodel_data,mode_compute=0,mode_file=usmodel_mode,first_obs=71,presample=4,lik_init=2,prefilter=0,mh_replic=0,mh_nblocks=2,mh_jscale=0.20,mh_drop=0.2);
shock_decomposition y;

estimation(optim=('MaxIter',200),datafile=usmodel_data,mode_compute=0,mode_file=usmodel_mode,first_obs=71,presample=4,lik_init=2,prefilter=0,mh_replic=250000,mh_nblocks=2,mh_jscale=0.20,mh_drop=0.2,bayesian_irf,irf=20,moments_varendo,filtered_vars,smoother);

// Generate IRfunctions and Moments of output, inflation and interest rate
   stoch_simul(irf=20) dy pinfobs robs ;

