//This program codes up the model in the notes: "The Labor Market in the New Keynesian Model:
// Incorporating a Simple DMP Version of the Labor Market and Rediscovering the Shimer Puzzle".
//By running the code twice (first with DMP=0, then with DMP=1) you get graphs that compare
//the response of inflation and output to an expansionary monetary policy shock.
//You can control whether the DMP version of the model has hiring (hiring=1) or search costs (hiring=0)
//in this comparison.

//there are several exercises one can do with this code: (i) examine the response to 
//technology shocks, (ii) look at the responses of other variables; (iii) consider alternative 
//parameterizations (for example, raise the replacement ratio from 0.4 to 0.99 and see how
//the responses in the DMP model get much stronger...this is the Hagedorn and Manovskii result). 
//(iv) you can use the code to investigate (by setting order=2) the impulse responses when
//the model is soved by second order approximation. In addition, in this case, you can investigate 
//the role of including the 'pruning' parameter, which is the correct thing to do when order=2.
//For a detailed discussion, see
// http://faculty.wcas.northwestern.edu/~lchrist/course/Korea_2012/lecture_on_solving_rev.pdf 


//if DMP = 1, then do the DMP model, otherwise competitive labor market model
@# define DMP = 1

//if hiring = 1, work with the hiring cost specification, if hiring = 0 the search cost specification
//this parameter is ignored if DMP=0.
@# define hiring = 1


hhiring=1;
@#if hiring == 0
   hhiring=0;
@#endif


@# if DMP == 1
   var l x Q J V U mm f vartheta tight Y v wbar; //13 variables
@# endif
@# if DMP == 0
var l;
@#endif

var C R pibar F K s pstar a u;// 9 variables
varexo eps_a eps_u;
parameters beta epsil theta Rsteady alpha phi_pie phi_x nu rho kap D eta sigm sig delta phi; 

sig=0.5;//power on unemployment in matching function
l_ss=0.95;//aggregate employment
recshare=0.01;//recruitment costs, as a share of GDP
repl=0.4;//the replacement ratio, D/wbar
beta=0.99;//household discount factor
rho=0.9;//probability of a match surviving for one period
Q_ss=0.7;//the vacancy filling rate
vartheta_ss=1;//competitive real price of what sellers in the labor market offer
alpha=0.8;//smoothing parameter in Taylor rule
phi_pie=1.5;//coefficient on inflation in the Taylor rule
phi_x=0.05;//coefficient on deviation of GDP from steady state in Taylor rule
delta=0;//autocorrelation of monetary policy shock
pibar_ss=1;//inflation (gross)
pstar_ss=1;//Tack Yun distortion
epsil=6;//elasticity of demand for specialized intermediat good.
nu=1-(epsil-1)/epsil;//subsidy to monopolists...careful, the code has not been debugged for other values of nu.
s_ss=(epsil-1)/epsil;//marginal cost of monopoly intermediate good producers.
theta=0.75;//price stickiness parameter for intermediat good producers

phi=1;//utility on labor parameter in case DMP=0

[R_ss,mm_ss,x_ss,f_ss,sigmx,Y_ss,v_ss,kapx,J_ss,wbar_ss,Dx,C_ss,U_ss,V_ss,etax,tight_ss,F_ss,K_ss] = ...;
    sstate(sig,l_ss,recshare,repl,beta,rho,Q_ss,vartheta_ss,theta,hhiring); 

Rsteady=R_ss;
kap=kapx;
D=Dx;
eta=etax;
sigm=sigmx;

fprintf('input parameters for steady state\n');
fprintf(' sig = %5.3f, l = %5.3f, recruiting cost as share of gross output = %5.3f, replacement ratio = %5.3f\n', ...
    sig,l_ss,recshare,repl);
fprintf(' beta = %5.3f, rho = %5.3f, Q = %5.3f, vartheta = %5.3f\n',beta,rho,Q_ss,vartheta_ss);

fprintf('computed parameters for steady state\n');
fprintf(' R = %15.13f, m = %15.13f, x = %15.13f, f = %15.13f, sigm  = %15.13f\n',R_ss,mm_ss,x_ss,f_ss,sigm)
fprintf(' Y = %15.13f, v = %15.13f, kap = %15.13f\n',Y_ss,v_ss,kap);
fprintf(' J = %15.13f, wbar (real wage)  = %15.13f, D = %15.13f, C = %15.13f\n',J_ss,wbar_ss,D,C_ss);
fprintf(' U = %15.13f, V = %15.13f, eta = %15.13f, tight = %15.13f\n',U_ss,V_ss,eta,tight_ss);


//R gross nominal rate of interest
//mm pricing kernel
//x hiring rate
//f probability of searching worker finding a job
//sigm parameter in matching function
//Y gross output of homogeneous final good
//v vacancy rate
//kap cost for posting one vacancy
//J value of a worker to a firm
//wbar real wage
//D payment to unemployed workers
//C consumption
//U value of unemployment
//V value of employment to a worker
//eta share of surplus going to worker
//tight labor market tightness (i.e., total vacancies / total number unemployed)
//F denominator of optimal price for intermediate good producers that optimize price
//K numerator 

@#if DMP == 0
  l_ss=1;
  C_ss=1;
@#endif

model;

//equations pertaining to the household and goods producing firms (7 equations)
//these 7 equations involve 8 variables: C,K,F,pibar,R,pstar,s,vartheta

    1/C - (beta/C(+1))*(R/pibar(+1));//intertemporal household equation
    1 + (pibar(+1)^(epsil-1))*beta*theta*F(+1) - F;//recursive representation for F
    F*(((1-theta*pibar^(epsil-1))/(1-theta))^(1/(1-epsil))) - K;//recursive representation for K
    (epsil/(epsil-1))*s + beta*theta*(pibar(+1)^epsil)*K(+1) - K;//static connection between F and K
// dynamic equation for pstar
    1/pstar - ((1-theta)*(((1-theta*(pibar)^(epsil-1))/(1-theta))^(epsil/(epsil-1))) + theta*(pibar^epsil)/pstar(-1));
// policy rule
    R/Rsteady = ((R(-1)/Rsteady)^alpha) * exp( (1-alpha) * ( phi_pie*(pibar-1) + phi_x*( log(C/steady_state(C)) ) ) + u );

@#if DMP == 0
    s=(1-nu)*C*l^phi*exp(-a);
    C=pstar*exp(a)*l;
@#endif
    
@#if DMP == 1

    s=(1-nu)*vartheta*exp(-a);

//equations that relate to the labor market and final goods market clearing 
//these 13 equations involve 12 new (i.e., not seen in the previous equations)
//variables: l, x, Q, J, wbar, mm, V, f, U, tight, Y, v
    l=(rho+x)*l(-1);//law of motion of employment (1)

@#if hiring == 0
    kap=Q*J;//free entry condition (2)
    C + v*l(-1)*kap - Y;//aggregate resources (10) 
@#endif

@#if hiring == 1
    kap=J;//free entry condition (2)'
    C + x*l(-1)*kap - Y;//aggregate resources (10)'
@#endif

    J=vartheta-wbar+rho*mm(+1)*J(+1);// value function of firm (3)
    V=wbar+mm(+1)*(rho*V(+1)+(1-rho)*(f(+1)*V(+1)+(1-f(+1))*U(+1)));//value of being a worker (4)
    U=D+mm(+1)*(f(+1)*V(+1)+(1-f(+1))*U(+1));//value of unemployment (5) 
    J=((1-eta)/eta)*(V-U);// Nash sharing rule (6)
    Q=sigm*tight^(-sig); //(7)
    tight=v*l(-1)/(1-rho*l(-1)); // definition of labor market tightness (8)
    mm=beta*C(-1)/C;// pricing kernel (9)
    f = x*l(-1)/(1-rho*l(-1));//finding rate (11) 
    Y = pstar*exp(a)*l;// final homogeneous output (12)
    Q = x/v;// (13)

@#endif

// exogenous shocks    
    a = rho*a(-1) + eps_a ;
    u = delta*u(-1) - eps_u;

end;

initval;

C=C_ss;
l=l_ss;
R=R_ss;
pibar=pibar_ss;
F=F_ss;
K=K_ss;
s=s_ss;
pstar=pstar_ss;

@#if DMP == 1
x=x_ss;
Q=Q_ss;
J=J_ss;
V=V_ss;
U=U_ss;
mm=mm_ss;
f=f_ss;
vartheta=vartheta_ss;
tight=tight_ss;
Y=Y_ss;
v=v_ss;
wbar=wbar_ss;
@#endif
a=0;
u=0;

end;

shocks;
var eps_a;
stderr .01;
var eps_u;
stderr .01;
end;

steady;

stoch_simul(order=1,irf=6,nograph);

@#if DMP == 0
infl=100*pibar_eps_u/pibar_ss;
gdp=100*C_eps_u/C_ss;
infl_tech=100*pibar_eps_a/pibar_ss;
gdp_tech=100*C_eps_a/C_ss;
save dmp infl gdp infl_tech gdp_tech
@#endif

@#if DMP == 1
load dmp infl gdp infl_tech gdp_tech
tt=1:length(gdp);

tech=0;// tech = 0 if you want to just look at responses to inflation, tech = 1 if also technology


ia=2;
ib=1;
if tech == 1
  ib=2;
end
 
subplot(ia,ib,1)
plot(tt,gdp,tt,100*C_eps_u/C_ss,'*-')
legend('competitive markets','DMP labor market')
title('consumption resp to policy')
ylabel('% deviation from ss')
axis tight

il=2;
if tech == 1
  il=3;
end

subplot(ia,ib,2)
plot(tt,infl,tt,100*pibar_eps_u/pibar_ss,'*-')
title('inflation resp to policy')
ylabel('% deviation from ss')
axis tight

hhiring
[100*C_eps_u/C_ss 100*pibar_eps_u/pibar_ss]

if tech == 1
subplot(ia,ib,4)
plot(tt,infl_tech,tt,100*pibar_eps_a/pibar_ss,'*-')
title('inflation resp to tech')
ylabel('% deviation from ss')
axis tight

subplot(ia,ib,2)
plot(tt,gdp_tech,tt,100*C_eps_a/C_ss,'*-')
title('con resp to tech')
ylabel('% deviation from ss')
axis tight

end

@#endif