clear all

%%Manipulation of the output
%choose 1 for seeing graphics, 0 otherwise 
graphics=0;
%choose which part of the output shall be depicted (number of periods)
visible_part=30;

%% Calibration
%Parameters
c_alpha_value=0.00797;
c_beta_value=0.99; 
c_pi_value=0.95;
c_rho_value=0.95;
c_W_B_value=0.0011;
c_W_H_value=0.045;

c_shock_epsilon=-0.1;

%{ 
%Taken over, to check
pbar=1;
delta_p=1;
phi=10;
spread=0.0025;
c_inno_p=0.01;
c_gamma=.99;
totmc=.1*.096;
%}

%Predetermined variables and/or end values
K_T_value=1;

%Steady state values
phi_ss_value=10; 
q_ss_value=1; 
spread_ss_value=0.0025;

%To allow for 'parameter_value' structure in dynare 
% while changing only values here in MatLab
c_alpha=c_alpha_value;
c_beta=c_beta_value; 
c_pi=c_pi_value;
c_rho=c_rho_value;
c_W_B=c_W_B_value;
c_W_H=c_W_H_value;
%c_inno_z=c_inno_z_value;
%c_sigma_epsilon=c_sigma_epsilon_value;
K_T=K_T_value;
phi_ss=phi_ss_value; 
q_ss=q_ss_value; 
spread_ss=spread_ss_value;

%Definition of variables in the steady state
r_bar_ss=1/c_beta; %eq.7 
r_B_ss=r_bar_ss+spread_ss;

Z_ss=(r_B_ss-1)*q_ss; %eq.14
K_H_ss=(1/c_alpha)*((Z_ss+q_ss)*c_beta-q_ss); %eq.9
K_B_ss=K_T-K_H_ss;

N_ss=q_ss*K_B_ss/phi_ss;
c_W_B=N_ss*(1-c_pi*((r_B_ss-r_bar_ss)*phi_ss+r_bar_ss)); %at the steady state, the endowment must be equal to the net worth of the exited bank; eq.17

Y_ss=Z_ss+Z_ss*c_W_H+c_W_B;
C_B_ss=(1-c_pi)*N_ss*(phi_ss*spread_ss+r_bar_ss); %eq. 26 + 10 + 15
C_H_ss=Y_ss-(c_alpha/2)*(K_H_ss)^2-C_B_ss;
c_theta=((1-c_pi)*c_beta*(r_bar_ss+phi_ss*spread_ss))/(phi_ss*(1-c_pi*c_beta*(r_bar_ss+phi_ss*spread_ss)));

mu_ss=c_beta*(1-c_pi+c_pi*phi_ss*c_theta)*spread_ss;
nu_ss=c_beta*(1-c_pi+c_pi*phi_ss*c_theta)*r_bar_ss;

D_ss=(1/r_bar_ss)*(K_B_ss-N_ss); %eq.10

%Structured variables
used_parameters=struct('theta',c_theta,'pi',c_pi,'alpha',c_alpha,'beta',c_beta,'Z_ss',Z_ss,'W_B',c_W_B,'W_H',c_W_H,'K_T',K_T,'rho',c_rho);


%% Recession
c_shock_epsilon_value=c_shock_epsilon;
dynare Paper_basic_DSGE noclearall

q_z=q;
K_H_z=K_H;
D_z=D;
phi_z=phi;
Y_z=Y;
N_z=N;
r_bar_z=r_bar;
C_H_z=C_H;
C_B_z=C_B;

%%Find the threshold price at each period of time 
%no bank run as long as the bank's capital is sufficient to cover deposits
q_threshold_no_run=D_z(1:end-1)./(1-K_H_z(1:end-1))-Z(2:end);

%% Graphical representation of the output
if graphics==1

    Z_dev=(Z(1:end)-Z_ss)*100/Z_ss;
    Y_dev=(Y(1:end)-Y_ss)*100/Y_ss;
    q_dev=(q(1:end)-q_ss)*100/q_ss;
    K_H_dev=(K_H(1:end)-K_H_ss)*100/K_H_ss;
    C_H_dev=(C_H(1:end)-C_H_ss)*100/C_H_ss;
    N_dev=(N(1:end)-N_ss)*100/N_ss;
    r_bar_dev=(r_bar(1:end)-r_bar_ss)*100/r_bar_ss;
    phi_dev=(phi(1:end)-phi_ss)*100/phi_ss;
    C_B_dev=(C_B(1:end)-C_B_ss)*100/C_B_ss;


    figure
    set(gcf,'numbertitle','off','name','Recession in the baseline model without bank runs')

    subplot(3,3,1)
    plot(Z_dev,'LineWidth',2)     
    title('Z')
    ylabel('%\Delta from steady state') 
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,2)
    plot(Y_dev,'LineWidth',2)     
    title('Y')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,3)
    plot(q_dev,'LineWidth',2)     
    title('q')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,4)
    plot(K_H_dev,'LineWidth',2)     
    title('K^H')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,5)
    plot(C_H_dev,'LineWidth',2)     
    title('C^H')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,6)
    plot(N_dev,'LineWidth',2)     
    title('N')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,7)
    plot(r_bar_dev,'LineWidth',2)     
    title('$\bar{r}$','interpreter','latex')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,8)
    plot(phi_dev,'LineWidth',2)     
    title('\phi')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,9)
    plot(C_B_dev,'LineWidth',2)     
    title('C^B')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])
    
end

%% Unanticipated run
%Simulation parameters
t_run=3; %period of the run
T=27; %periods of simulation

% Initialization of output
C_H_run=zeros(T,1);
epsilon_t=zeros(T,1);
q_z_run=ones(T,1);
path_Z_shocked=zeros(T,1);
q_run=zeros(T,1);
%% Simulation

for j=1:T
    %calculate the shock in the (T-j+1)th period 
    %(as a shock to the steady state Z)
    %i.e. find values from the steady state backwards to the run
    c_shock_epsilon_value=c_rho^(T-j)*c_shock_epsilon;
    %evaluate the model with this shock
    dynare Paper_bank_run_DSGE noclearall
    %in order to derive consumption and asset prices when the bank is not active
    %first the shocks in the preceding two periods are calculated 
    Z_shocked=exp(-c_rho.^(T-j:T-j+1)'*(1-c_rho))*Z_ss
    %the path of Z is written (most recent on top)
    %then, the HH consumption in the previous period is calculated
    %given that no banks operated  (W_B=C_B and K_H=K_T)
    C_H_onlyHH=exp(-c_rho.^(T-j)*(1-c_rho))*Z_ss-(c_alpha/2)*(K_T^2);   
    %asset prices
    q_z_run(2:end)=q(2:T)
    q_z_run(1)=c_beta*C_H_onlyHH/C_H(2)*(Z_shocked(2)+q_z_run(2))-c_alpha*K_T
    q_run(T-j+1)=q_z_run(1)
    %reaching the actual run
    %all the variables have the same structure
    %starting with a steady state they experience the shock 
    %then they have the normal path until the run
    %then there is the first period after the run (special for some)
    %finally, the path convergence back to the steady state
    if j==T-t_run+1
        q_run=[q_ss;q_z(2:t_run);q_z_run;q_ss];
        K_H_run=[K_H_ss;K_H_z(2:t_run);K_T;K_H(2:T+1)];
        D_run=[D_ss; D_z(2:t_run);0;D(2:T+1)];
        N_run=[N_ss;N_z(2:t_run);c_W_B;N(2:T+1)];
        C_H_run=[C_H_ss;C_H_z(2:t_run);C_H_onlyHH;C_H(t_run+2:T+1)];
        C_B_run=[C_B_ss;C_B_z(2:t_run);0;C_B(2:T+1)];
        Y_run=[Y_ss;Y_z(2:t_run);C_H_onlyHH+c_W_B+(c_alpha/2)*(K_T^2);Y(2:T-1)];
        phi_run=[phi_ss;phi_z(2:t_run);0;phi(2:T+1)];
        r_bar_run=[r_bar_ss;r_bar_z(2:t_run);1/c_beta*(C_H_run(t_run+2:end)./C_H_run(t_run+1:end-1));r_bar_ss];
    end

end
%% Graphical representation of the output
if graphics==1

    Z_dev_run=(path_Z_shocked(1:end)-Z_ss)/Z_ss;
    Y_dev_run=(Y_run(1:end)-Y_ss)/Y_ss;
    q_dev_run=(q_run(1:end)-q_ss)/q_ss;
    K_H_dev_run=(K_H_run(1:end)-K_H_ss)/K_H_ss;
    C_H_dev_run=(C_H_run(1:end)-C_H_ss)/C_H_ss;
    N_dev_run=(N_run(1:end)-N_ss)/N_ss;
    r_bar_dev_run=(r_bar_run(1:end)-r_bar_ss)/r_bar_ss;
    phi_dev_run=(phi_run(1:end)-phi_ss)/phi_ss;
    C_B_dev_run=(C_B_run(1:end)-C_B_ss)/C_B_ss;


   
     figure
    set(gcf,'numbertitle','off','name','Recession in the baseline model without bank runs')

    subplot(3,3,1)
    plot(Z_dev_run,'LineWidth',2)   
    hold on
    plot(Z_dev,'r');
    title('Z')
    ylabel('%\Delta from steady state') 
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,2)
    plot(Y_dev_run,'LineWidth',2) 
    hold on
    plot(Y_dev,'r');
    title('Y')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,3)
    plot(q_dev_run,'LineWidth',2)
    hold on
    plot(q_dev,'r');
    title('q')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,4)
    plot(K_H_dev_run,'LineWidth',2) 
    hold on
    plot(K_H_dev,'r');
    title('K^H')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,5)
    plot(C_H_dev_run,'LineWidth',2)
    hold on
    plot(C_H_dev,'r');
    title('C^H')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,6)
    plot(N_dev_run,'LineWidth',2)  
    hold on
    plot(N_dev,'r');
    title('N')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,7)
    plot(r_bar_dev_run,'LineWidth',2)    
    hold on
    plot(r_bar_dev,'r');
    title('$\bar{r}$','interpreter','latex')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,8)
    plot(phi_dev_run,'LineWidth',2) 
    hold on
    plot(phi_dev,'r');
    title('\phi')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])

    subplot(3,3,9)
    plot(C_B_dev_run,'LineWidth',2) 
    hold on
    plot(C_B_dev,'r');
    title('C^B')
    ylabel('%\Delta from steady state')
    axis([-inf,visible_part,-inf,inf])
end