%% %%%%%%%%%Dynare code for CCF's (2012) paper %%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        V A R I A B L E S, P A R A M E T E R S  &  S H O C K S  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%% Declaration of variables, parameters & shocks for exog. eq's %%%%%%
%Exogenous process for z(t)
parameters rho_z;
var z;
varexo eps_z;
%Exogenous processes for varphi's
parameters rho_varphi;
var varphi_u varphi_r;
varexo eps_varphi_u eps_varphi_r;
%Process for mu
var mu;
parameters rho_mu;
varexo eps_mu;
%Processes for b hats
parameters rho_b;
var bh_u bh_r;
varexo eps_b_u eps_b_r;
%Process for G_hat
parameters rho_g;
var Gh;
varexo eps_g;
%Process for risk premia shocks
parameters rho_zeta;
var eps_zeta;
varexo nu_zeta;
%Cost-push shock
var lambda_f_t;
varexo eps_lambda;
%Other shocks
varexo eps_m eps_T eps_B;

%%%%%%%%%%%%%%%%%%%%% Definition of  model parameters %%%%%%%%%%%%%%%%%%%%% 
parameters  alpha beta_u beta_r betabar zeta_p omega_u omega_r lambda_f
            Csi_ratio delta gamma zeta h adp Sdp sigma_u sigma_r zeta_w
            nu lambda_w chi_wu csi_pu qu phi_T rho_B rho_r phi_pi phi_y
            zetaprime;



%Steady-state variables
parameters  rk RL kappa b_u b_r varphi_u_ss varphi_r_ss B_L B G T lnPi Y C_u C_r 
            I Kbar coef1 coef2;






%%%%%%%%%%%%%%%%%%%%%% Declaration of model variables %%%%%%%%%%%%%%%%%%%%%
%Marginal cost
var mch rhk wh;
%Capital demand
var Kh Lh;
%Technological constraint
var Yh;
%Pricing equations
var Csih_u Csih_r Xh_pnu Xh_pnr Xh_pdu Xh_pdr;
%Prices LOM
var pi;
%Effective capital
var Kbh uh;
%Capital LOM
var Ih;
%Capital utilisation
%Q LOM
var qh;
%Marginal utilities for each type of consumer
var Ch_u Ch_r;
%Euler equations
var r zetah rL;
%Wage setting
var Xh_wnu Xh_wnr Xh_wdu Xh_wdr;
%Budget constraint
var Bh Bh_L Th;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          P A R A M E T E R S 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 0.33;
delta = 0.025;

gamma = 1.9922/400;
lnPi = 2.2404/400;
beta_u = 1/(1+0.4943/400);
zeta = 0.4761/400;

beta_r = beta_u/(1+zeta);


Sdp = 4.4277;
adp = 0.2093;
h = 0.8370;
sigma_u = 3.0151;
sigma_r = 1.5635;
zetaprime = 0.2420/100;
omega_u = 0.9832;
omega_r = 1 - omega_u;

Csi_ratio = 0.7917;
C_ratio = 0.7641;
chi_wu = 0.5170;
nu = 1.6814;
zeta_w = 0.6860;
zeta_p = 0.9260;
phi_T = 1.2660;
rho_r = 0.8547;
phi_pi = 1.5644;
phi_y = 0.2975;

rho_z = 0.1286;
rho_mu = 0.8293;
rho_b = 0.9603;
rho_varphi = 0.52;
rho_B = 0.9773;
rho_zeta = 0.9614;
rho_g = 0.7581;

b_u = 10;
b_r = 5;
varphi_u_ss = 20;
varphi_r_ss = 2;

R = beta_u^(-1)*exp(gamma)*exp(lnPi);
betabar = (omega_u*beta_u*Csi_ratio+omega_r*beta_r)/(omega_u*Csi_ratio+omega_r);

%steady-state values
Y = 1;
RL = R*(1+zeta);
rk = betabar^(-1)*exp(gamma)-(1-delta);
lambda_f = 1.5;
lambda_w = 1.5;
kappa = RL-RL/30;


B = 0.16;
B_L = 0.16;
G = 0.2;
T = G - (1-beta_u^(-1))*B - (1/(RL-kappa)-RL/(RL-kappa)*1/(exp(gamma)*exp(lnPi)))*B_L;



mc = 1/(1+lambda_f);
w_til = (1+lambda_f)^(-1/(1-alpha))*alpha^(alpha/(1-alpha))*(1-alpha);
w = w_til * rk^(-alpha/(1-alpha));
K_til = alpha/(1+lambda_f);
K = K_til * rk^(-1);
Kbar = exp(gamma) * K_til/rk;
I = (exp(gamma)-(1-delta))*K_til/rk;
C_r = (1-I-G)/(omega_u * C_ratio + omega_r);
C_u = C_ratio * C_r;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                C O M P O S I T E      P A R A M E T E R S
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Composite parameter definition...
qu = (betabar/beta_r-1)*zeta^(-1);
csi_pu = omega_u/(omega_u + omega_r * 1/Csi_ratio * (1-beta_u*zeta_p)/(1-beta_r*zeta_p));
%chi_wu = omega_u/(omega_u + omega_r *(b_u*varphi_u_ss/(b_r*varphi_r_ss) * Csi_ratio)^(1/(lambda_w+(1+lambda_w)*nu)));
coef1 = (B_L/B)/(RL-kappa); 
coef2 = (1-exp(-gamma)*(exp(lnPi))^(-1)*kappa)*RL/(RL-kappa);


%predetermined_var Kh;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           M   O   D   E   L   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;
%%%%%%%%%%%%%%%%%%%% Equations for exogenous processes %%%%%%%%%%%%%%%%%%%%
    z = rho_z*z(-1)+eps_z;
    varphi_u = rho_varphi * varphi_u(-1) + eps_varphi_u;
    varphi_r = rho_varphi * varphi_r(-1) + eps_varphi_r;
    mu = rho_mu*mu(-1)+eps_mu;
    bh_u = rho_b*bh_u(-1) + eps_b_u;
    bh_r = rho_b*bh_r(-1) + eps_b_r;
    Gh = rho_g*Gh(-1) + eps_g;
    eps_zeta = rho_zeta*eps_zeta(-1) + nu_zeta;
    lambda_f_t = eps_lambda;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Model equations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    mch = alpha * rhk + (1-alpha) * wh; %ok
    Kh = wh - rhk + Lh; %ok
    Yh = alpha*Kh + (1-alpha)*Lh; %ok
    Xh_pnu = (1-beta_u*zeta_p) * (Csih_u + Yh + lambda_f_t + mch)  + beta_u*zeta_p * ( (1+lambda_f)/lambda_f * EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(Xh_pnu(+1)) ); %ok
    Xh_pnr = (1-beta_r*zeta_p) * (Csih_r + Yh + lambda_f_t + mch)  + beta_r*zeta_p * ( (1+lambda_f)/lambda_f * EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(Xh_pnr(+1)) ); %ok
    Xh_pdu = (1-beta_u*zeta_p) * (Csih_u + Yh)  + beta_u*zeta_p * ( 1/lambda_f * EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(Xh_pdu(+1)) ); %ok
    Xh_pdr = (1-beta_r*zeta_p) * (Csih_r + Yh)  + beta_r*zeta_p * ( 1/lambda_f * EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(Xh_pdr(+1)) ); %ok
    
    pi = (1-zeta_p)/zeta_p * (csi_pu * (Xh_pnu - Xh_pdu) + (1-csi_pu) * (Xh_pnr - Xh_pdr)); %ok
    Kh = -  z + uh + Kbh(-1); %ok
    Kbh = (1-delta)*exp(-gamma)*(Kbh(-1)-z) + (1-(1-delta)*exp(-gamma))*(mu+Ih); %ok
    rhk = adp/rk * uh; %ok
    %%qh = - EXPECTATION(0)(z(+1)) + betabar * exp(-gamma)*(rk*EXPECTATION(0)(rhk(+1)) + (1-delta) * EXPECTATION(0)(qh(+1))) + qu * ((1+zeta)/(1+zeta*qu) * EXPECTATION(0)(Csih_u(+1))-Csih_u) + (1-qu) * (1/(1+qu*zeta) * EXPECTATION(0)(Csih_r(+1))-Csih_r);     
    qh = - EXPECTATION(0)(z(+1)) + betabar * exp(-gamma)*(rk*EXPECTATION(0)(rhk(+1)) + (1-delta) * EXPECTATION(0)(qh(+1))) + (omega_u*beta_u*Csi_ratio/(omega_u*beta_u*Csi_ratio+omega_r*beta_r)) * (EXPECTATION(0)(Csih_u(+1))-betabar/beta_u * Csih_u) + (omega_r*beta_r/(omega_u*beta_u*Csi_ratio+omega_r*beta_r)) * (EXPECTATION(0)(Csih_r(+1))-betabar/beta_r * Csih_r);     
    %0 = qh + mu - Sdp * (z + Ih - Ih(-1)) + betabar * Sdp * exp(2*gamma) * (EXPECTATION(0)(z(+1)) + EXPECTATION(0)(Ih(+1)) - Ih); %ok
    0 = qh + mu - exp(2*gamma) * Sdp *(z+Ih-Ih(-1)) - betabar*exp(2*gamma)*Sdp*(EXPECTATION(0)(z(+1)) + EXPECTATION(0)(Ih(+1)) - Ih);
    Csih_u = 1/(1-beta_u*h) * (bh_u - beta_u * h * EXPECTATION(0)(bh_u(+1)) - sigma_u/(1-h) * ((1+beta_u*h^2)* Ch_u - beta_u*h*EXPECTATION(0)(Ch_u(+1)) - h*Ch_u(-1)));
    Csih_r = 1/(1-beta_r*h) * (bh_r - beta_r * h * EXPECTATION(0)(bh_r(+1)) - sigma_r/(1-h) * ((1+beta_r*h^2)* Ch_r - beta_r*h*EXPECTATION(0)(Ch_r(+1)) - h*Ch_r(-1)));
    
    Csih_u = r + EXPECTATION(0)(Csih_u(+1)) - EXPECTATION(0)(z(+1)) - EXPECTATION(0)(pi(+1)); %ok
    zetah + Csih_u = (RL/(RL-kappa))*rL + (EXPECTATION(0)(Csih_u(+1)) - EXPECTATION(0)(z(+1)) - EXPECTATION(0)(pi(+1)) - kappa/(RL-kappa) * EXPECTATION(0)(rL(+1)));
            Csih_r = (RL/(RL-kappa))*rL + (EXPECTATION(0)(Csih_r(+1)) - EXPECTATION(0)(z(+1)) - EXPECTATION(0)(pi(+1)) - kappa/(RL-kappa) * EXPECTATION(0)(rL(+1)));
    Xh_wnu = (1-zeta_w*beta_u) * (1+lambda_w) * (bh_u + varphi_u + (1+nu)*Lh + (1+lambda_w)/lambda_w * (1+nu) * wh) + zeta_w*beta_u * ((1+lambda_w)/lambda_w * (1+nu) * ( EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(z(+1)) ) + EXPECTATION(0)(Xh_wnu(+1)));
    Xh_wnr = (1-zeta_w*beta_r) * (1+lambda_w) * (bh_r + varphi_r + (1+nu)*Lh + (1+lambda_w)/lambda_w * (1+nu) * wh) + zeta_w*beta_r * ((1+lambda_w)/lambda_w * (1+nu) * ( EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(z(+1)) ) + EXPECTATION(0)(Xh_wnr(+1)));
    Xh_wdu = (1-zeta_w*beta_u) * (Csih_u + Lh + (1+lambda_w)/lambda_w * wh) + zeta_w*beta_u * (1/lambda_w * ( EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(z(+1)) ) + EXPECTATION(0)(Xh_wdu(+1)));
    Xh_wdr = (1-zeta_w*beta_r) * (Csih_r + Lh + (1+lambda_w)/lambda_w * wh) + zeta_w*beta_r * (1/lambda_w * ( EXPECTATION(0)(pi(+1)) + EXPECTATION(0)(z(+1)) ) + EXPECTATION(0)(Xh_wdr(+1)));
    wh = (1-zeta_w)/(1+nu*(1+lambda_w)/lambda_w) * ( chi_wu * (Xh_wnu - Xh_wdu) + (1-chi_wu)*(Xh_wnr - Xh_wdr) ) +  zeta_w * ( wh(-1) - pi - z );
    Bh + coef1 * Bh_L = beta_u^(-1) * (Bh(-1) +r(-1)) + coef1 * beta_r^(-1) * Bh_L(-1) + coef2*coef1* rL + G/B * Gh - Y/B * Th - (beta_u^(-1)+coef1*beta_r^(-1)) * (z + pi); 
    Bh_L - RL/(RL-kappa)*rL = rho_B * (Bh_L(-1) - RL/(RL-kappa)*rL(-1)) + eps_B; 
    (Th-Gh*G)/(T-G) = phi_T * (Bh(-1) + coef1 * Bh_L(-1) - RL/(RL-kappa) * coef1 * rL(-1))/(1+coef1) + eps_T;
    r = rho_r * r(-1) + (1-rho_r)*(phi_pi * pi + phi_y *(Yh - Yh(-4) + z + z(-1) + z(-2) + z(-3))) + eps_m;
    zetah = zetaprime * Bh_L + eps_zeta;
    Yh = omega_u*C_u/Y * Ch_u + omega_r*C_r/Y * Ch_r + I/Y * Ih + G/Y * Gh + exp(-gamma) * rk * Kbar/Y * uh;

end;

resid(1);
steady;


shocks;
%   Syntax
%      var xx;
%      stderr xx;
    
    %%var B_L;
    %periods 1 2 3 4:10;
    %values -0.0625 -0.125 -0.1875 -0.25;
    %var eps_z;
    %stderr 1;
    %var eps_varphi_u;
    %stderr 1;
    %var eps_varphi_r;
    %stderr 1;
    %var eps_mu;
    %stderr 1;
    %var eps_b_u; 
    %stderr 1;
    %var eps_b_r;
    %stderr 1;
    %var eps_g;
    %stderr 1;
    %var nu_zeta;
    %stderr 1;
    var eps_lambda;
    stderr 1;
    %var eps_m;
    %stderr 1;
    %var eps_T;
    %stderr 1;
    %var eps_B;
    %stderr 1;

end;
check;


stoch_simul(periods=300, irf=100);