close all
clear all
clc
format long

addpath C:\dynare\4.4.2\mex\matlab\win64-7.8-8.2\
%% ===========================================================
% Preliminaries 
%=============================================================
order = 1;  %order of approximation  
Ts = 300001;  %simulation length
Per = 200;

% Calibrated parameters (Quarterly)
betta       = .995;
rho         =  0.3;
alfa        =  10;
del_y       = 0.5;
del_pi      = 1.968;  
g_ss        = 1.6696;
phi_g       = 0.9790947;
phi_y       = -0.0445;  
phi_x       = 0.978; 
c           = 0.0;
sig_x       = 0.0003432;
sig_y       = 0.0078;
sig_c       = 0.6*0.0078;
sig_g       = 0.01097;   
var_shock   =2; 

gam_0    = 0.1;
L        = 0.3;



%% ===========================================================
% Steady State and other parameters 
%=============================================================
b_ss        = log(0.8);       %The steady state of debt/gdp   %(3) 
annual_inf  = 2.5/100;
pi_ss       = log(1+(annual_inf/4));%(1)    %Arbitrary
eta_0       = -1.4;           %This is arbitrary
def_prob    = 0.01;
eps_bar_ss  = 2*def_prob-1;
del_0       =-log(betta)+(1-del_pi)*pi_ss-def_prob*L;
q_ss        = log(betta)-pi_ss-(1+eps_bar_ss)*L/2;
kappa       = exp(q_ss)-(1/exp(pi_ss));
eta_b       = log((exp(-g_ss)-kappa*exp(b_ss))/exp(eta_0))/log(exp(b_ss));

diff=1; 
crit = 0.0000001;
gu = 10000;
gd = 0;
gam_0_b=2;
while abs(diff)>crit;
    
    gam_g       =(-gam_0_b+eta_b+b_ss*eta_b)/(g_ss*(1-phi_g)+eps_bar_ss*sig_g+1);
    eps_check   = (1/(gam_g*sig_g))*(eta_0+eta_b*b_ss-gam_0_b-gam_g*(1-phi_g)*g_ss-gam_g*phi_g*g_ss);
    h_check     = 0.5+eps_check/2;
  

    diff = (h_check-def_prob);
    if diff>0 
         gd =gam_0_b;
    else
         gu = gam_0_b;
    end;
    gam_0_b=(1/2)*(gu+gd);
end;
gam_0       =gam_0_b
eps_bar_ss  = 2*def_prob-1;                  %(1)
h_ss        = def_prob ;                     %(2)
q_ss        = log(betta)-pi_ss-(1+eps_bar_ss)*L/2; %(3)
tau_ss      = eta_0+eta_b*b_ss;   %(4) 
dy_ss       =0;                   %(5) 
dc_ss       =0;                   %(6) 
g_ss        =g_ss;                %(7) 
x_ss        =0;                   %(8)
m_ss        =log(betta);          %(9)
r_ss        =-m_ss;               %(10)
v_ss        = -log(1-exp(m_ss));  %(11)
b_ss        = log(0.8);            %(12)
pi_ss       = log(1+(annual_inf/4)); %(13)

           
%% ===========================================================
% Write Dynare++ Source File: 'ENDO.mod' 
%=============================================================

% Create file
fid = fopen('CSS_def.mod','w+');
fprintf(fid, 'var tau dy dc x q pai m r g b v eps_bar h;\n');
fprintf(fid, 'varexo z;\n');
fprintf(fid, '\n');
fprintf(fid, 'parameters betta rho alfa eta_0 eta_b del_0 del_pi del_y phi_x phi_y phi_g sig_x sig_c sig_y sig_g gam_0 gam_g L\n');
fprintf(fid, 'pi_ss q_ss b_ss tau_ss dy_ss dc_ss g_ss x_ss m_ss r_ss v_ss eps_bar_ss h_ss;\n');
fprintf(fid, '\n');
fprintf(fid, 'betta = %3.6f;\n',betta);
fprintf(fid, 'rho  = %3.6f;\n',rho);
fprintf(fid, 'alfa = %3.6f;\n',alfa);
fprintf(fid, 'del_y = %3.6f;\n',del_y);
fprintf(fid, 'del_pi = %3.6f;\n',del_pi);
fprintf(fid, 'del_0 = %3.6f;\n',del_0);
fprintf(fid, 'eta_0 = %3.6f;\n',eta_0);
fprintf(fid, 'eta_b = %3.6f;\n',eta_b);
fprintf(fid, 'phi_x = %3.6f;\n',phi_x);
fprintf(fid, 'phi_y = %3.6f;\n',phi_y);
fprintf(fid, 'phi_g = %3.6f;\n',phi_g);
fprintf(fid, 'sig_x = %3.6f;\n',sig_x);
fprintf(fid, 'sig_c = %3.6f;\n',sig_c);
fprintf(fid, 'sig_g = %3.6f;\n',sig_g);
fprintf(fid, 'sig_y = %3.6f;\n',sig_y);
fprintf(fid, 'gam_0 = %3.6f;\n',gam_0_b);
fprintf(fid, 'gam_g = %3.6f;\n',gam_g);
fprintf(fid, 'L = %3.6f;\n',L);

fprintf(fid, 'pi_ss = %3.6f;\n',pi_ss);
fprintf(fid, 'q_ss  = %3.6f;\n',q_ss);
fprintf(fid, 'b_ss = %3.6f;\n', b_ss);
fprintf(fid, 'tau_ss = %3.6f;\n',tau_ss);
fprintf(fid, 'dy_ss = %3.6f;\n',dy_ss);
fprintf(fid, 'dc_ss = %3.6f;\n',dc_ss);
fprintf(fid, 'g_ss = %3.6f;\n',g_ss);
fprintf(fid, 'x_ss = %3.6f;\n',x_ss);
fprintf(fid, 'm_ss = %3.6f;\n',m_ss);
fprintf(fid, 'v_ss = %3.6f;\n',v_ss);
fprintf(fid, 'r_ss = %3.6f;\n',r_ss);
fprintf(fid, 'eps_bar_ss = %3.6f;\n',eps_bar_ss);
fprintf(fid, 'h_ss = %3.6f;\n',h_ss);

fprintf(fid, 'model;\n');
fprintf(fid, ' exp(tau(-1))+exp(q+b) =   exp(b(-1)-pai-dy)+exp(-g);\n');
fprintf(fid, ' tau                   =   eta_0+eta_b*b  ;\n');
fprintf(fid, ' dy                    =   phi_y*(tau(-1)-tau_ss)+ sig_y*z ;\n');
fprintf(fid, ' dc                    =   x(-1) + sig_c*z;\n');
fprintf(fid, ' x                     =   phi_x*x(-1)+sig_x*z ;\n');
fprintf(fid, ' g                     =   (1-phi_g)*g_ss+phi_g*g(-1)+sig_g*z;\n');
fprintf(fid, ' q                     =   -del_0-del_pi*pai-del_y*dy ;\n');
fprintf(fid, ' eps_bar               =  (1/(gam_g*sig_g))*(eta_0+eta_b*b - gam_0-gam_g*(1-phi_g)*g_ss- gam_g*phi_g*g ) ;\n');
fprintf(fid, '     h                 =   (1/2)+(eps_bar/2) ;\n');
fprintf(fid, ' exp(q)                =  exp(m(+1)-pai(+1))*exp((1+eps_bar)*(L/2));\n');
fprintf(fid, ' exp(m)                = (betta^(alfa/rho))*exp(dc*(alfa*(rho-1)/rho))*exp(r*(alfa-rho)/rho);\n');
fprintf(fid, ' exp(r)                = (exp(v)*exp(dc))/(exp(v)-1);\n');
fprintf(fid, ' exp(v)                = 1 + exp(m(+1)+v(+1)+dc(+1));\n');
fprintf(fid, 'end; \n');
fprintf(fid, '\n');
fprintf(fid, 'initval;\n');
fprintf(fid, 'tau = %3.12f;\n',tau_ss);
fprintf(fid, 'dy = %3.12f;\n', dy_ss);
fprintf(fid, 'dc = %3.12f;\n', dc_ss);
fprintf(fid, 'x = %3.12f;\n',  x_ss);
fprintf(fid, 'q = %3.12f;\n',  q_ss);
fprintf(fid, 'pai = %3.12f;\n',pi_ss);
fprintf(fid, 'm = %3.12f;\n',  m_ss);
fprintf(fid, 'r = %3.12f;\n',r_ss) ;
fprintf(fid, 'g = %3.12f;\n',g_ss);
fprintf(fid, 'b = %3.12f;\n',b_ss);
fprintf(fid, 'v = %3.12f;\n',v_ss);
fprintf(fid, 'eps_bar = %3.12f;\n',eps_bar_ss);
fprintf(fid, 'h = %3.12f;\n',(h_ss));
fprintf(fid, 'end;\n');
fprintf(fid, '\n');
fprintf(fid, 'vcov = [%3.12f];\n',var_shock);
fprintf(fid, '\n');
fprintf(fid, 'order = %i;\n',order);
fclose(fid);
%% ===========================================================
% Solve Model using Dynare++
%=============================================================

% Housekeeping
%keep Ts
if exist('CSS_def.mat')
    delete CSS_def.mat  % prevent errantly using old .mat file when code doesn't solve
end

% Call executable
disp('Running...')
    !dynare++ CSS_def.mod
disp('Done!')

% Load workfile 
load CSS_def.mat


%%============================================================
% Steady States
%=============================================================
dyn_dy_ss= dyn_steady_states(1,1);
dyn_q_ss= dyn_steady_states(2,1);
dyn_r_ss= dyn_steady_states(3,1);
dyn_tau_ss= dyn_steady_states(4,1);
dyn_x_ss= dyn_steady_states(5,1);
dyn_g_ss= dyn_steady_states(6,1);
dyn_b_ss= dyn_steady_states(7,1);
dyn_dc_ss= dyn_steady_states(8,1);
dyn_pi_ss= dyn_steady_states(9,1);
dyn_m_ss= dyn_steady_states(10,1);
dyn_v_ss= dyn_steady_states(11,1);
%% ===========================================================
% Simulation
%=============================================================

% Simulate variables in 'CSS.mod'

shocks = zeros(1,Ts)./0;
Vdyn = dynare_simul('CSS_def.mat',shocks);


%% ===========================================================
% Compute & Print Summary Statistics
%=============================================================

fprintf('\n');
fprintf('Moments (QUARTERLY)');
fprintf('\n');

% OUTPUT
dy_mean = mean(exp(Vdyn(dyn_i_dy,:)))-1;
dy_std = std(exp(Vdyn(dyn_i_dy,:)));
dy_AC1 = autocorr(Vdyn(dyn_i_dy,:),1);

fprintf('\n')
fprintf('E(dy):                 %3.6f\n', dy_mean);
fprintf('Std(dy):               %3.6f\n', dy_std);
fprintf('AC1(dy):               %3.6f\n', dy_AC1);

% TAX
tau_mean = mean(exp(Vdyn(dyn_i_tau,:)));
tau_std = std(exp(Vdyn(dyn_i_tau,:)));
tau_AC1 = autocorr(Vdyn(dyn_i_tau,:),1);

fprintf('\n')
fprintf('E(Tau)):                %3.6f\n', tau_mean);
fprintf('Std(tau):               %3.6f\n', tau_std);
fprintf('AC1(tau):               %3.6f\n', tau_AC1);


% DEBT
b_mean  = mean(exp(Vdyn(dyn_i_b,:)));
b_std = std(exp(Vdyn(dyn_i_b,:)));
b_AC1 = autocorr(Vdyn(dyn_i_b,:),1);

fprintf('\n')
fprintf('E(B)):                %3.6f\n', b_mean);
fprintf('Std(B):               %3.6f\n', b_std);
fprintf('AC1(B):               %3.6f\n', b_AC1);

% BOND PRICE 
q_mean  = mean(exp(Vdyn(dyn_i_q,:)));
q_std = std(exp(Vdyn(dyn_i_q,:)));
q_AC1 = autocorr(Vdyn(dyn_i_q,:),1);

fprintf('\n')
fprintf('E(Q)):                %3.6f\n', q_mean);
fprintf('Std(Q):               %3.6f\n', q_std);
fprintf('AC1(Q):               %3.6f\n', q_AC1);


% INFLATION
pi_mean  = mean(exp(Vdyn(dyn_i_pai,:)));
pi_std = std(exp(Vdyn(dyn_i_pai,:)));
pi_AC1 = autocorr(Vdyn(dyn_i_pai,:),1);

fprintf('\n')
fprintf('E(PI)):                %3.6f\n', pi_mean);
fprintf('Std(PI):               %3.6f\n', pi_std);
fprintf('AC1(PI):               %3.6f\n', pi_AC1);


%% CORRELATIONS
cov_dy_tau=cov(Vdyn(1,:), Vdyn(4,:));
cov_dy_pi=cov(Vdyn(1,:),  Vdyn(9,:));
cov_dy_q=cov(Vdyn(1,:),  Vdyn(2,:));
cov_dy_g=cov(Vdyn(1,:),  Vdyn(6,:));
cov_pi_q=cov(Vdyn(9,:), Vdyn(2,:));
cov_pi_b=cov(Vdyn(9,:), Vdyn(7,:));


corr_dy_tau =  cov_dy_tau(2,1) /(std(Vdyn(1,:))*std(Vdyn(4,:)))
corr_dy_pi  =  cov_dy_pi(2,1) / (std(Vdyn(1,:))*std(Vdyn(9,:)))
corr_dy_q   =  cov_dy_q(2,1) / (std(Vdyn(1,:))*std(Vdyn(2,:)))
corr_dy_g   =  cov_dy_g(2,1) / (std(Vdyn(1,:))*std(Vdyn(6,:)))
corr_pi_q   =  cov_pi_q(2,1) / (std(Vdyn(9,:))*std(Vdyn(2,:)))
corr_pi_b   =  cov_pi_b(2,1) / (std(Vdyn(9,:))*std(Vdyn(7,:)))