%% Compute steady state at benchmark case
clear all;

parameter=calibration;
BETA=parameter(1); DELTA=parameter(2); THETA_S=parameter(3); THETA_P=parameter(4);
ALPHA=parameter(5); As=parameter(6); Ap=parameter(7); ETA=parameter(8); GAMMA=parameter(9);
TAU_S=parameter(10); TAU_P=parameter(11); SIGMA=parameter(12); OMEGA_I=parameter(13); MU=parameter(14);
ZETA =parameter(15);EPSILON=parameter(16);PI=parameter(17); RL=parameter(18); 
KM=parameter(19); BETA_B=parameter(20); BETA_E=parameter(21);PHI=parameter(22); OMEGAM=parameter(32);
ALPHA_S=parameter(23);ALPHA_P =parameter(24);OMEGA_P=parameter(25); ALPHA_SB=parameter(30);
KAPPA=parameter(26); XI=parameter(27);KAPPA_M=parameter(28);BETA_SB=parameter(29); ZETAL=parameter(31);


[lcr,fval,cflagr]=fsolve(@(x) solve_steady(x(1),x(2),x(3),x(4),x(5),PI,DELTA,THETA_S,THETA_P,ALPHA,As,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA_E,BETA,EPSILON,BETA_SB,ALPHA_SB,ZETAL),[7.2,6.12,0.924,0.22,0.39],optimset('TolFun',1e-10));
lp=lcr(1);ls=lcr(2);omegas=lcr(3);hp=lcr(4);pp=lcr(5);
[y,yp,ys,ch,h,d,m,ks,wage,DIV_S,hs,kp,lp_bar,DIV_P,ps,Rs,Rp,Rp_bar,rb,Xb,F_xi,PI_CB,PI_SB,b,k,i,DWL,totutility,bg]=solve_star(lp,ls,omegas,hp,pp,PI,DELTA,THETA_S,THETA_P,ALPHA,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA,EPSILON,KAPPA_M,BETA_SB,ALPHA_SB,ZETAL);
pr_IN1=F_xi; int_xi_IN1=KM/(KM-1)*(OMEGAM-OMEGAM^KM*omegas^(1-KM));omegas_IN1=omegas;
save INITIAL_INFO pr_IN1 int_xi_IN1 omegas_IN1 


%% generate table III1 in appendix 
vTAU_P=[0.4 0.7 ];   
vCAR=[0.1 0.15];
for t=1:1:length(vCAR)
        GAMMA=vCAR(t);
for j=1:1:length(vTAU_P)
    TAU_P=vTAU_P(j);
    if j==1
  [lcr,fval,cflagr(t,j)]=fsolve(@(x) solve_steady(x(1),x(2),x(3),x(4),x(5),PI,DELTA,THETA_S,THETA_P,ALPHA,As,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA_E,BETA,EPSILON,BETA_SB,ALPHA_SB,ZETAL),[7.2,6.12,0.924,0.22,0.39],optimset('TolFun',1e-10));
    else 
  [lcr,fval,cflagr(t,j)]=fsolve(@(x) solve_steady(x(1),x(2),x(3),x(4),x(5),PI,DELTA,THETA_S,THETA_P,ALPHA,As,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA_E,BETA,EPSILON,BETA_SB,ALPHA_SB,ZETAL),[lp_ss(t,j-1),ls_ss(t,j-1),omegas_ss(t,j-1),hp_ss(t,j-1),pp_ss(t,j-1)],optimset('TolFun',1e-10));
    end
   lp=lcr(1);ls=lcr(2);omegas=lcr(3);hp=lcr(4);pp=lcr(5);
[y,yp,ys,ch,h,d,m,ks,wage,DIV_S,hs,kp,lp_bar,DIV_P,ps,Rs,Rp,Rp_bar,rb,Xb,F_xi,PI_CB,PI_SB,b,k,i,DWL,utility,bg]=solve_star(lp,ls,omegas,hp,pp,PI,DELTA,THETA_S,THETA_P,ALPHA,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA,EPSILON,KAPPA_M,BETA_SB,ALPHA_SB,ZETAL);
y_ss(t,j)=y;ys_ss(t,j)=ys;yp_ss(t,j)=yp;F_xi_ss(t,j)=F_xi;omegas_ss(t,j)=omegas;
k_ss(t,j)=k;i_ss(t,j)=i; lls_ss(t,j)=lp/ls; yys_ss(t,j)=yp/ys;hp_ss(t,j)=hp;pp_ss(t,j)=pp;
Rs_ss(t,j)=Rs;Rp_ss(t,j)=Rp;ch_ss(t,j)=ch;m_ss(t,j)=m;
ls_ss(t,j)=ls;lp_ss(t,j)=lp;lp_bar_ss(t,j)=lp_bar;
utility_ss(t,j)=utility;rb_ss(t,j)=rb;b_ss(t,j)=b;Rp_bar_ss(t,j)=Rp_bar;
con_u(t,j)=(log(ch)-KAPPA*h^(1+ETA)/(1+ETA))/(1-BETA);gdp_ss(t,j)=ch+i;
end 
end
TABLE1_appendix=[utility_ss(1,1) gdp_ss(1,1) lls_ss(1,1) yys_ss(1,1);              
                 utility_ss(2,1) gdp_ss(2,1) lls_ss(2,1) yys_ss(2,1);
                 utility_ss(1,2) gdp_ss(1,2) lls_ss(1,2) yys_ss(1,2); 
                 utility_ss(2,2) gdp_ss(2,2) lls_ss(2,2) yys_ss(2,2)]

save table1_appendix TABLE1_appendix
%% generate transition path from tau_p=0.4 to tau_p=0.7

clear all

    GAMMA_P=0.1;
    CHANGE=0.7;
    save PARAMETER_CHOICE  CHANGE
    dynare base_transition
    mean_gdp =mean(gdp(1:200));
    vol_gdp =var(gdp(1:200));
    mean_utility =mean(totutility(1:200));
    mean_ch =mean(ch(1:200));
    mean_lls =mean(lls(1:200));
    mean_yys =mean(yys(1:200));
    vol_ch =var(ch(1:200));
    vol_lls =var(lls(1:200));
    vol_yys =var(yys(1:200));
    mean_b =mean(b(1:200));
    mean_lpbar =mean(lp_bar(1:200));
    vol_b =var(b(1:200));
    vol_lpbar =var(lp_bar(1:200));
    
    Welfare0(1)=utility(1);
    for t=1:1:200
Welfare0(t+1)=Welfare0(t)+BETA^(t)*utility(t+1);
    end
Welfare =Welfare0(201);

save Transition_welfare 

%% generate fig III1 in appendix
clear
load Transition_welfare
figure;

subplot(2,3,1)
plot(gdp(1:200),'k-','LineWidth',1.5)
title('GDP','FontSize',14,'interpreter','latex')
subplot(2,3,2)
plot(ch(1:200),'k-','LineWidth',1.5)
title('消费','FontSize',14,'interpreter','latex')
subplot(2,3,3)
plot(i(1:200),'k-','LineWidth',1.5)
title('投资','FontSize',14,'interpreter','latex')
subplot(2,3,4)
plot(lp_bar(1:200),'k-','LineWidth',1.5)
title('影子银行贷款规模','FontSize',14,'interpreter','latex')
subplot(2,3,5)
plot(lls(1:200),'k-','LineWidth',1.5)
title('信贷结构','FontSize',14,'interpreter','latex')
subplot(2,3,6)
plot(yys(1:200),'k-','LineWidth',1.5)
title('产业结构','FontSize',14,'interpreter','latex')
xlabel('时期','FontSize',14,'interpreter','latex')


%% generate table III2 in appendix
clear all

    GAMMA_P=0.15;
    CHANGE=0.7;
    save PARAMETER_CHOICE  CHANGE
    dynare base_transition
    mean_gdp1 =mean(gdp(1:200));
    vol_gdp1 =var(gdp(1:200));
    mean_utility1 =mean(totutility(1:200));
    mean_ch1 =mean(ch(1:200));
    mean_lls1 =mean(lls(1:200));
    mean_yys1 =mean(yys(1:200));
    vol_ch1 =var(ch(1:200));
    vol_lls1 =var(lls(1:200));
    vol_yys1 =var(yys(1:200));
    mean_b1 =mean(b(1:200));
    mean_lpbar1 =mean(lp_bar(1:200));
    vol_b1 =var(b(1:200));
    vol_lpbar1 =var(lp_bar(1:200));
    
    Welfare0(1)=utility(1);
    for t=1:1:200
Welfare0(t+1)=Welfare0(t)+BETA^(t)*utility(t+1);
    end
Welfare1 =Welfare0(201);

load Transition_welfare 
TABLEIII2=[Welfare Welfare1;
           mean_gdp mean_gdp1;
           mean_lls  mean_lls1;
           mean_yys mean_yys1;
           mean_lpbar mean_lpbar1;
           vol_gdp.*10000 vol_gdp1.*10000;
           vol_lls.*10000 vol_lls1.*10000;
           vol_yys.*10000 vol_yys1.*10000]
save tableIII2 TABLEIII2