%There are 3 variables (Graphs, KN2010cost and WPparam)
%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.
% If WPparam=1, Working paper's parameters are used.
%If WPparam=0, published paper's parameters are used.

@#define Graphs = 1

@#define KN2010cost = 1

@#define WPparam = 0


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Endogenous Variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
var
C_1 $C$ % consumption.
N_1 $N$ % labour.
w_1 $w$ % real wage.
Pi_1 $\Pi$ % inflation CPI.
X_1 $X$ % markup.
Y_1 $Y$ % production.
R_1 $R$ % Nominal policy rate.
g_1 $g$ % growth techno shock.
A_1 $A$ %level techno shock
lambda_1 ${\lambda}$ %marginal utility of consumption.
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.
log_C1
;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%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$
qbar $\bar q$

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



@#if WPparam == 1
    beta_1=0.99;%discount factor.
    gamma_1=109.82;%weight of disutility of labour.
    sigmaN_1 = 1;%Frisch elasticity.
    alpha_1 = 0.36;%share of capital in production.
    delta_1 =0.02;%depreciation rate.
    rhog_1=0.95;%persistence of growth techno shock.
    rhoA_1=0.95;%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.25;%(1-kappa)=probability of price change.
    eta_1 = 0;%price indexation.
    epsi_1=21;% elasticity of substitution between goods (=21 because SS of markup =1.05)
    rhor_1 = 0.73;%persistence of nominal interest rate.
    rhopi_1=1.5;%weight of inflation in monetary policy.
    rhoy_1 = 0.2;%weight of output growth in monetary policy.
@#else
    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 nominal interest rate.
    rhopi_1=1.95;%weight of inflation in monetary policy.
    rhoy_1 = 0.18;%weight of output growth in monetary policy.
@#endif


%abar_1 and bbar_1 are function of delta_1 and sigmaphi_1 and such that PHI(delta)=delta and PHI(0)=0.
%The first value of sigmaphi is given in a code on Nutahara website and the second one is in the paper.
@#if WPparam == 1
    sigmaphi_1 = 1.01;%investment adjustment cost parameter.
    abar_1 =   0.0118;%0.0066;%param in investment adjustment cost function.
    bbar_1 =   0.0896;%0.0803;%param in investment adjustment cost function.
    qbar = 1;%
@#else
    sigmaphi_1 = 1.01;%investment adjustment cost parameter.
    abar_1 =   0.0148;%0.0066;%param in investment adjustment cost function.
    bbar_1 =   0.1064;%0.0803;%param in investment adjustment cost function.
    qbar = 1;%
@#endif


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

  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;
%1 Marginal utility of consumption:
    lambda_1 = 1/C_1;

%2 real wage:
    gamma_1*C_1*N_1^(sigmaN_1) = w_1;

%3 Investissement:
    %q_1  = beta_1*lambda_1(+1)/lambda_1*(1+g_1(+1))^(-1)*(q_1(+1)*(1-delta_1) + rk_1(+1) 
    %+ q_1(+1)*(Phi_1(+1)-DPhi_1(+1)*I_1(+1)/K_1*(1+g_1(+1))));
    q_1  = beta_1*lambda_1(+1)/lambda_1*(exp(g_1(+1)))^(-1)*(q_1(+1)*(1-delta_1) + rk_1(+1) 
    + q_1(+1)*(Phi_1(+1)-DPhi_1(+1)*I_1(+1)/K_1*exp(g_1(+1))));

%4 Euler equation:
	1 = beta_1*lambda_1(+1)/lambda_1*(exp(g_1(+1)))^(-1)*R_1(+1)/Pi_1(+1);

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

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

%7 marginal productivity of capital:
    rk_1 = alpha_1/X_1*Y_1/K_1(-1)*(exp(g_1));

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

%9 Law of motion of Capital:
    K_1 = ((1-delta_1)*K_1(-1) + Phi_1*K_1(-1))/(exp(g_1));

%10 New philips curve:
    log(Pi_1) = beta_1/(1+eta_1*beta_1)*log(Pi_1(+1)) + eta_1/(1+eta_1*beta_1)*log(Pi_1(-1)) - (1-kappa_1)*(1-kappa_1*beta_1)/(kappa_1*(1+eta_1*beta_1))*log(X_1/(epsi_1/(epsi_1-1)));

%11 market clearing:
    C_1 + I_1 = Y_1;

%12 Monetary policy:
    log(R_1*beta_1) =log(R_1(-1)*beta_1)*(rhor_1) + (log(Pi_1(+1))*rhopi_1 + log(Y_1/Y_1(-1)*(exp(g_1)))*rhoy_1)*(1-rhor_1);
    %log(R_1*beta_1) =log(R_1(-1)*beta_1)*(rhor_1) + (log(Pi_1(+1))*rhopi_1 + log(Y_1/Y_1(-1))*rhoy_1)*(1-rhor_1);
%13 14 15 16 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;

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

log_C1= log(C_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;
X_1 = epsi_1/(epsi_1-1);
Pi_1 = 1;
R_1 = 1/beta_1;
Phi_1 = delta_1;
DPhi_1 = sigmaphi_1*delta_1/qbar/(delta_1+abar_1);
q_1 = 1/DPhi_1;
rk_1 = q_1/beta_1 - q_1*(1-delta_1) - q_1*delta_1*(1-DPhi_1);% need q
K_1 = (rk_1*X_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);
I_1 = K_1*delta_1;% need K
C_1 = Y_1-I_1;%need y and I
lambda_1 = C_1^(-1);% need c
w_1 = gamma_1*C_1*N_1^(sigmaN_1);% need c
log_C1= log(C_1);
end;



resid(1);
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=40,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');
indlambda = strmatch('lambda_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');
indX = strmatch('X_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)
%accumulate trend shock.
C_1_cum = oo_.irfs.C_1_nuug_1/oo_.steady_state(indC) + cumsum(oo_.irfs.g_1_nuug_1);
plot([0 C_1_cum(1:19)]*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)
%accumulate trend shock.
K_1_cum = oo_.irfs.K_1_nuug_1/oo_.steady_state(indK) + cumsum(oo_.irfs.g_1_nuug_1);
plot([0 K_1_cum(1:19)]*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)
%accumulate trend shock.
Y_1_cum = oo_.irfs.Y_1_nuug_1/oo_.steady_state(indY) + cumsum(oo_.irfs.g_1_nuug_1);
plot([0 Y_1_cum(1:19)]*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.lambda_1_nuug_1(1:19)]/oo_.steady_state(indlambda)*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)
%accumulate trend shock.
I_1_cum = oo_.irfs.I_1_nuug_1/oo_.steady_state(indI) + cumsum(oo_.irfs.g_1_nuug_1);
plot([0 I_1_cum(1:19)]*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 oo_.irfs.X_1_nuug_1(1:19)]/oo_.steady_state(indX)*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)
%accumulate trend shock.
w_1_cum = oo_.irfs.w_1_nuug_1/oo_.steady_state(indw) + cumsum(oo_.irfs.g_1_nuug_1);
plot([0 w_1_cum(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real wage')

@#if KN2010cost == 1
    @#if WPparam==1
        print -dpdf Figure1_WPparam;
    @#else
        print -dpdf Figure1_Paperparam;
    @#endif
@#else
    @#if WPparam==1
        print -dpdf Figure1_WPparam_mycost;
    @#else
        print -dpdf Figure1_Paperparam_mycost;
    @#endif
@#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.lambda_1_nuuA_1(1:19)]/oo_.steady_state(indlambda)*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 oo_.irfs.X_1_nuug_1(1:19)]/oo_.steady_state(indX)*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 WPparam==1
        print -dpdf Figure3_WPparam;
    @#else
        print -dpdf Figure3_Paperparam;
    @#endif
@#else
    @#if WPparam==1
        print -dpdf Figure3_WPparam_mycost;
    @#else
        print -dpdf Figure3_Paperparam_mycost;
    @#endif
@#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)
%accumulate trend shock.
C_1_cum = oo_.irfs.C_1_eg_1/oo_.steady_state(indC) + cumsum(oo_.irfs.g_1_eg_1);
plot([0 C_1_cum(1:19)]*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)
%accumulate trend shock.
K_1_cum = oo_.irfs.K_1_eg_1/oo_.steady_state(indK) + cumsum(oo_.irfs.g_1_eg_1);
plot([0 K_1_cum(1:19)]*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)
%accumulate trend shock.
Y_1_cum = oo_.irfs.Y_1_eg_1/oo_.steady_state(indY) + cumsum(oo_.irfs.g_1_eg_1);
plot([0 Y_1_cum(1:19)]*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.lambda_1_eg_1(1:19)]/oo_.steady_state(indlambda)*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)
%accumulate trend shock.
I_1_cum = oo_.irfs.I_1_eg_1/oo_.steady_state(indI) + cumsum(oo_.irfs.g_1_eg_1);
plot([0 I_1_cum(1:19)]*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 oo_.irfs.X_1_eg_1(1:19)]/oo_.steady_state(indX)*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)
%accumulate trend shock.
w_1_cum = oo_.irfs.w_1_eg_1/oo_.steady_state(indw) + cumsum(oo_.irfs.g_1_eg_1);
plot([0 w_1_cum(1:19)]*100)
hold on
line('XData', [0 20], 'YData', [0 0])
hold off
title('Real wage')


@#if KN2010cost == 1
    @#if WPparam==1
        print -dpdf ContempGrowthShock_WPparam;
    @#else
        print -dpdf ContempGrowthShock_Paperparam;
    @#endif
@#else
    @#if WPparam==1
        print -dpdf ContempGrowthShock_WPparam_mycost;
    @#else
        print -dpdf ContempGrowthShock_Paperparam_mycost;
    @#endif
@#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.lambda_1_eA_1(1:19)]/oo_.steady_state(indlambda)*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 oo_.irfs.X_1_eg_1(1:19)]/oo_.steady_state(indX)*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 WPparam==1
        print -dpdf ContempLevelShock_WPparam;
    @#else
        print -dpdf ContempLevelShock_Paperparam;
    @#endif
@#else
    @#if WPparam==1
        print -dpdf ContempLevelShock_WPparam_mycost;
    @#else
        print -dpdf ContempLevelShock_Paperparam_mycost;
    @#endif
@#endif


@#endif