%There are 3 variables (NeworOld, Graphs, KN2010cost)
%NeworOld=1 if we want New parameters and 0 if old.
%HN2010cost=1 if one wants investment adjustment cost like Kobayashi and Nutahara 2010.
%if HN2010cost=0 use of another cost function
%Graphs=1 then graphics are saved in pdf files.

@#define NeworOld = 0

@#define Graphs = 1

@#define KN2010cost = 1



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Endogenous Variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
var
C_1 $C$ % consumption.
N_1 $N$ % labour.
w_1 $w$ % real wage.
Pi_1 $\Pi$ % inflation CPI.
mc_1 $mc$ % real marginal cost.
Y_1 $Y$ % production.
R_1 $R$ % Nominal policy rate.
g_1 $g$ % growth techno shock.
A_1 $A$ %level techno shock
x1_1 $x1$ %first recurrence for price rigidity.
x2_1 $x2$ %second recurrence for price rigidity.
ptilde_1 $\tilde p$ %relative optimal price.
s_1 $s$ %price dispersion.
uc_1 ${u_C}$ %marginal utility of consumption.
ul_1 ${u_N}$ %marginal disutility of labour.
rk_1 $rk$ %real rental rate.
Phi_1 $\Phi$ %adjustment cost function.
DPhi_1 $D\Phi$ %derivative of the adjustment cost function.
I_1 $I$ %Investment.
K_1 $K$ %Capital stock.
q_1 $q$ %Tobin's Q.
nug_1 ${nu\_inter\_g}$ %intermediate variable growth techno shock.
nuA_1 ${nu\_inter\_A}$ %intermediate variable level techno shock.
;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Exogenous Variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
varexo
eA_1 $eA$% contemporary level techno shock
nuuA_1 ${\nu_A}$ % news level techno shock
eg_1 $eg$ %contemporary growth techno shock.
nuug_1 ${\nu_g}$ % news growth techno shock.
;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parameters 
beta_1 $\beta$
gamma_1 $\gamma$
sigmaN_1 ${\sigma_N}$
alpha_1 $\alpha$
sigmaphi_1 ${\sigma_{\phi}}$
delta_1 $\delta$
rhog_1 ${\rho_g}$
rhoA_1 $\rho_A$
gss_1 $g_{ss}$
Ass_1 $A_{ss}$
kappa_1 $\kappa$
eta_1 $\eta$
epsi_1 $\epsilon$
rhor_1 $\rho_r$
rhopi_1 $\rho_{\pi}$
rhoy_1 $\rho_y$
abar_1 $\bar a$
bbar_1 $\bar b$

%useful if choose alternative adjustment cost function
phik_1 $\phi_k$
   
;

beta_1=1.01358^(-0.25);%discount factor.
gamma_1=109.82;%weight of disutility of labour.
sigmaN_1 = 1;%Frisch elasticity.
alpha_1 = 0.4;%share of capital in production.
delta_1 =0.025;%depreciation rate.
rhog_1=0.83;%persistence of growth techno shock.
rhoA_1=0.83;%persistence of level techno shock.
gss_1=0;% steady state of growth techno.
Ass_1=0;%steady state of level techno.% log normal
kappa_1=1-0.36;%(1-kappa)=probability of price change.
eta_1 = 0.84;%price indexation.
epsi_1=6;% elasticity of substitution between goods (=6 because SS of markup =1.2)
rhor_1 = 0.81;%persistence of norminal interest rate.
rhopi_1=1.95;%weight of inflation in monetary policy.
rhoy_1 = 0.18;%weight of output growth in monetary policy.



%abar_1 and bbar_1 are function of delta_1 and sigmaphi_1 and such that PHI(delta)=delta and PHI(0)=0.
%if sigmaphi=0.02 then abar and bbar equal 0.01 and 0.092 respectively.
%if sigmaphi=1.01 then abar and bbar equal 0.9976 and  0.0025.
%The first value of sigmaphi is given in a code on Nutahara website and the second one is in the paper.
@#if NeworOld == 1
    sigmaphi_1 = 1.01;%investment adjustment cost parameter.
    abar_1 =  0.9976;%param in investment adjustment cost function.
    bbar_1 =  0.0025;%param in investment adjustment cost function.
@#else
    sigmaphi_1 = 0.02;%investment adjustment cost parameter.
    abar_1 =  0.0100;%param in investment adjustment cost function.
    bbar_1 =  0.0920;%param in investment adjustment cost function.
@#endif


phik_1 = 200;%useful if choose alternative adjustment cost function

  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;
%1 real wage:
    w_1 = -ul_1/uc_1; 

%2 Euler equation:
	uc_1 = beta_1*uc_1(+1)*(1+g_1(+1))^(-1)*R_1(+1)/Pi_1(+1); 

%3 First recurrence, price rigidity:
    x1_1 = Y_1*mc_1 + kappa_1*beta_1*uc_1(+1)/uc_1*Pi_1(+1)^(epsi_1+1)*Pi_1^(-eta_1*epsi_1)*x1_1(+1); 

%4 Second recurrence, price rigidity:
    x2_1 = Y_1 + kappa_1*beta_1*uc_1(+1)/uc_1*Pi_1(+1)^(epsi_1)*Pi_1^(-eta_1*epsi_1)*x2_1(+1);

%5 Relative price for non-optimizing firms:
    ptilde_1=epsi_1/(epsi_1-1)*x1_1/x2_1;

%6 Inflation:
    1= kappa_1*Pi_1^(epsi_1-1)*Pi_1(-1)^(eta_1*(1-epsi_1)) + (1-kappa_1)*(ptilde_1)^(1-epsi_1);

%7 Price dispersion:
    s_1 = (1-kappa_1)*(ptilde_1)^(-epsi_1) + kappa_1*(Pi_1/Pi_1(-1)^(eta_1))^(epsi_1)*s_1(-1);
    
%8 Monetary policy:
    R_1*beta_1 =(R_1(-1)*beta_1)^(rhor_1)*(Pi_1(+1)^rhopi_1*(Y_1/Y_1(-1)*(1+g_1))^rhoy_1)^(1-rhor_1);

%9 Production:
	Y_1 = exp(A_1)*(K_1(-1)/(1+g_1))^(alpha_1)*N_1^(1-alpha_1)/s_1;

%10 11 12 13 Productivity:
    %level
    A_1 = (1-rhoA_1)*Ass_1 + rhoA_1*A_1(-1) + 0.01*(eA_1 + nuA_1(-4));
    nuA_1 = nuuA_1;
    %growth
    g_1 = (1-rhog_1)*gss_1 + rhog_1*g_1(-1) + 0.01*(eg_1 + nug_1(-4));
    nug_1 = nuug_1;

%14 Marginal cost:
    w_1 = mc_1*(1-alpha_1)*Y_1/N_1;

%15 marginal productivity of capital:
    rk_1 = alpha_1*mc_1*Y_1/K_1(-1)*(1+g_1);

%16 Capital Adjustment cost:
    @#if KN2010cost == 1
        Phi_1 = (1-delta_1) + sigmaphi_1*log(I_1/K_1(-1)*(1+g_1)+abar_1) + bbar_1;
    @#else
        Phi_1 = I_1/K_1(-1)*(1+g_1) + (1-delta_1) - phik_1/2*(I_1/K_1(-1)*(1+g_1)-delta_1)^2;
    @#endif
%17 Derivative of Capital Adjustment cost:
    @#if KN2010cost == 1
        DPhi_1 = sigmaphi_1/(I_1/K_1(-1)*(1+g_1)+abar_1);
    @#else
        DPhi_1 = 1 - phik_1*(I_1/K_1(-1)*(1+g_1)-delta_1);
    @#endif
%18 Law of motion of Capital:
    K_1 = Phi_1*K_1(-1)/(1+g_1);

%19 Tobin's Q:
    q_1 = 1/DPhi_1;

%20 Investissement:
    q_1  = beta_1*uc_1(+1)/uc_1*(1+g_1(+1))^(-1)*(rk_1(+1) 
    + q_1(+1)*(Phi_1(+1)-DPhi_1(+1)*I_1(+1)/K_1*(1+g_1(+1))));

%21 Marginal utility of consumption:
    uc_1 = 1/C_1;

%22 Marginal utility of labour:
    ul_1 = -gamma_1*N_1^(sigmaN_1);

%23 Goods market clearing:
    C_1 + I_1 = Y_1;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initval;
N_1=0.1;%guess
eA_1= 0;
nuA_1= 0;
nuuA_1= 0;
eg_1= 0;
nug_1= 0;
nuug_1= 0;
A_1 = Ass_1;
g_1 = gss_1;
ptilde_1=1;
mc_1 = (epsi_1-1)/epsi_1;
s_1=1;
Pi_1 = 1;
R_1 = 1/beta_1;
Phi_1 = 1;
DPhi_1 = 1;
q_1 = 1;
rk_1 = q_1/beta_1 - q_1*(1-delta_1);% need q
ul_1=-gamma_1*N_1^(sigmaN_1);
K_1 = (rk_1/(mc_1*alpha_1*exp(A_1)*N_1^(1-alpha_1)))^(1/(alpha_1-1));
Y_1 = exp(A_1)*K_1^(alpha_1)*N_1^(1-alpha_1)/s_1;
I_1 = K_1*delta_1;% need K
C_1 = Y_1-I_1;%needd y and I
uc_1 = C_1^(-1);% need c
w_1 = -ul_1/uc_1;% need c
mc_1 = w_1*K_1^(-alpha_1)*N_1^(alpha_1)/((1-alpha_1)*exp(A_1));
x1_1 = Y_1*mc_1/(1-kappa_1*beta_1); %need y and mc
x2_1 = Y_1/(1-kappa_1*beta_1);%need y 
end;

steady;
resid(1);
check;


shocks;
var nuug_1=1;
var nuuA_1=1;
var eg_1 = 1;
var eA_1 = 1;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stoch_simul(order=1,irf=20,nograph);


write_latex_static_model;
write_latex_dynamic_model;

@#if Graphs == 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Index for graph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
indC = strmatch('C_1',M_.endo_names,'exact');
indK = strmatch('K_1',M_.endo_names,'exact');
indY = strmatch('Y_1',M_.endo_names,'exact');
indPi = strmatch('Pi_1',M_.endo_names,'exact');
indq = strmatch('q_1',M_.endo_names,'exact');
indR = strmatch('R_1',M_.endo_names,'exact');
induc = strmatch('uc_1',M_.endo_names,'exact');
indrk = strmatch('rk_1',M_.endo_names,'exact');
indI = strmatch('I_1',M_.endo_names,'exact');
indN = strmatch('N_1',M_.endo_names,'exact');
indmc = strmatch('mc_1',M_.endo_names,'exact');
indw = strmatch('w_1',M_.endo_names,'exact');




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Figure 1 of Kobayashi et Nutahara
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(4,3,1)
plot([0 oo_.irfs.C_1_nuug_1(1:19)]/oo_.steady_state(indC)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Consumption')
subplot(4,3,2)
plot([0 oo_.irfs.K_1_nuug_1(1:19)]/oo_.steady_state(indK)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Capital')
subplot(4,3,3)
plot([0 oo_.irfs.Y_1_nuug_1(1:19)]/oo_.steady_state(indY)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Output')

subplot(4,3,4)
plot([0 oo_.irfs.Pi_1_nuug_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Inflation')
subplot(4,3,5)
plot([0 oo_.irfs.q_1_nuug_1(1:19)]/oo_.steady_state(indq)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Tobins Q')
subplot(4,3,6)
plot([0 oo_.irfs.R_1_nuug_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Nominal Interest Rate')

subplot(4,3,7)
plot([0 oo_.irfs.uc_1_nuug_1(1:19)]/oo_.steady_state(induc)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Lambda')
subplot(4,3,8)
plot([0 oo_.irfs.rk_1_nuug_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real rental rate')
subplot(4,3,9)
plot([0 oo_.irfs.I_1_nuug_1(1:19)]/oo_.steady_state(indI)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Investment')

subplot(4,3,10)
plot([0 oo_.irfs.N_1_nuug_1(1:19)]/oo_.steady_state(indN)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Labour')
subplot(4,3,11)
plot([0 (1./(oo_.irfs.mc_1_nuug_1(1:19)+oo_.steady_state(indmc))-1/oo_.steady_state(indmc))]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Markup')
subplot(4,3,12)
plot([0 oo_.irfs.w_1_nuug_1(1:19)]/oo_.steady_state(indw)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real wage')

@#if KN2010cost == 1
    @#if NeworOld == 1
        print -dpdf Figure1_NewSigmaphi;
    @#else
        print -dpdf Figure1_OldSigmaphi;
    @#endif
@#else
    print -dpdf Figure1_myCost;
@#endif

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Figure 3 of Kobayashi et Nutahara
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(4,3,1)
plot([0 oo_.irfs.C_1_nuuA_1(1:19)]/oo_.steady_state(indC)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Consumption')
subplot(4,3,2)
plot([0 oo_.irfs.K_1_nuuA_1(1:19)]/oo_.steady_state(indK)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Capital')
subplot(4,3,3)
plot([0 oo_.irfs.Y_1_nuuA_1(1:19)]/oo_.steady_state(indY)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Output')

subplot(4,3,4)
plot([0 oo_.irfs.Pi_1_nuuA_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Inflation')
subplot(4,3,5)
plot([0 oo_.irfs.q_1_nuuA_1(1:19)]/oo_.steady_state(indq)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Tobins Q')
subplot(4,3,6)
plot([0 oo_.irfs.R_1_nuuA_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Nominal Interest Rate')

subplot(4,3,7)
plot([0 oo_.irfs.uc_1_nuuA_1(1:19)]/oo_.steady_state(induc)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Lambda')
subplot(4,3,8)
plot([0 oo_.irfs.rk_1_nuuA_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real rental rate')
subplot(4,3,9)
plot([0 oo_.irfs.I_1_nuuA_1(1:19)]/oo_.steady_state(indI)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Investment')

subplot(4,3,10)
plot([0 oo_.irfs.N_1_nuuA_1(1:19)]/oo_.steady_state(indN)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Labour')
subplot(4,3,11)
plot([0 (1./(oo_.irfs.mc_1_nuuA_1(1:19)+oo_.steady_state(indmc))-1/oo_.steady_state(indmc))]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Markup')
subplot(4,3,12)
plot([0 oo_.irfs.w_1_nuuA_1(1:19)]/oo_.steady_state(indw)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real wage')

@#if KN2010cost == 1
    @#if NeworOld == 1
        print -dpdf Figure3_NewSigmaphi;
    @#else
        print -dpdf Figure3_OldSigmaphi;
    @#endif
@#else
    print -dpdf Figure3_myCost;
@#endif


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Kobayashi et Nutahara contemporary growth shock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(4,3,1)
plot([0 oo_.irfs.C_1_eg_1(1:19)]/oo_.steady_state(indC)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Consumption')
subplot(4,3,2)
plot([0 oo_.irfs.K_1_eg_1(1:19)]/oo_.steady_state(indK)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Capital')
subplot(4,3,3)
plot([0 oo_.irfs.Y_1_eg_1(1:19)]/oo_.steady_state(indY)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Output')

subplot(4,3,4)
plot([0 oo_.irfs.Pi_1_eg_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Inflation')
subplot(4,3,5)
plot([0 oo_.irfs.q_1_eg_1(1:19)]/oo_.steady_state(indq)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Tobins Q')
subplot(4,3,6)
plot([0 oo_.irfs.R_1_eg_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Nominal Interest Rate')

subplot(4,3,7)
plot([0 oo_.irfs.uc_1_eg_1(1:19)]/oo_.steady_state(induc)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Lambda')
subplot(4,3,8)
plot([0 oo_.irfs.rk_1_eg_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real rental rate')
subplot(4,3,9)
plot([0 oo_.irfs.I_1_eg_1(1:19)]/oo_.steady_state(indI)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Investment')

subplot(4,3,10)
plot([0 oo_.irfs.N_1_eg_1(1:19)]/oo_.steady_state(indN)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Labour')
subplot(4,3,11)
plot([0 (1./(oo_.irfs.mc_1_eg_1(1:19)+oo_.steady_state(indmc))-1/oo_.steady_state(indmc))]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Markup')
subplot(4,3,12)
plot([0 oo_.irfs.w_1_eg_1(1:19)]/oo_.steady_state(indw)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real wage')



@#if KN2010cost == 1
    @#if NeworOld == 1
        print -dpdf ContempGrowthShock_NewSigmaphi;
    @#else
        print -dpdf ContempGrowthShock_OldSigmaphi;
    @#endif
@#else
    print -dpdf ContempGrowthShock_myCost;
@#endif

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Kobayashi et Nutahara contemporary growth shock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(4,3,1)
plot([0 oo_.irfs.C_1_eA_1(1:19)]/oo_.steady_state(indC)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Consumption')
subplot(4,3,2)
plot([0 oo_.irfs.K_1_eA_1(1:19)]/oo_.steady_state(indK)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Capital')
subplot(4,3,3)
plot([0 oo_.irfs.Y_1_eA_1(1:19)]/oo_.steady_state(indY)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Output')

subplot(4,3,4)
plot([0 oo_.irfs.Pi_1_eA_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Inflation')
subplot(4,3,5)
plot([0 oo_.irfs.q_1_eA_1(1:19)]/oo_.steady_state(indq)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Tobins Q')
subplot(4,3,6)
plot([0 oo_.irfs.R_1_eA_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Nominal Interest Rate')

subplot(4,3,7)
plot([0 oo_.irfs.uc_1_eA_1(1:19)]/oo_.steady_state(induc)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Lambda')
subplot(4,3,8)
plot([0 oo_.irfs.rk_1_eA_1(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real rental rate')
subplot(4,3,9)
plot([0 oo_.irfs.I_1_eA_1(1:19)]/oo_.steady_state(indI)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Investment')

subplot(4,3,10)
plot([0 oo_.irfs.N_1_eA_1(1:19)]/oo_.steady_state(indN)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Labour')
subplot(4,3,11)
plot([0 (1./(oo_.irfs.mc_1_eA_1(1:19)+oo_.steady_state(indmc))-1/oo_.steady_state(indmc))]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Markup')
subplot(4,3,12)
plot([0 oo_.irfs.w_1_eA_1(1:19)]/oo_.steady_state(indw)*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real wage')


@#if KN2010cost == 1
    @#if NeworOld == 1
        print -dpdf ContempGrowthlevel_NewSigmaphi;
    @#else
        print -dpdf ContempGrowthlevel_OldSigmaphi;
    @#endif
@#else
    print -dpdf ContempGrowthlevel_myCost;
@#endif

@#endif