% Note: The chain provided with the code in mhall_RBC.mat has 25000 draws for file size considerations. 
% To produce the full posterior distribution the chain length should be set to S_.Options.EST.N = 750000 
% in line 20 of run_mcmc.m

addpath functions

%% Load chains

load('mhall_RBC.mat'); 

[nobs, T_sample] = size(Zobs);
v_ = var_dim ;
n  = v_.dim.n    ;
n_ = v_.dim.n_    ;
ns = v_.dim.ns    ;

%% Run smoother to get shocks and initial conditions (at the mode)

[ytT, n_T] = smoother(params_xx_mode,params_T_out_mode,S_);
Tz         = params_T_out_mode(1,:);

%% Run counterfactula with no break in trend growth

% Parameters

S_.param = param_func(params_xx_mode);
    
 % Initial guess for steady state

S_.initvals = [0.834755321 % ln c -- 1
    -0.095030773 % ln invest -- 2
    -0.973976034 % ln lambda -- 3
    0.005224883 % ln r_F -- 4
    4.631985754 % ln k -- 5
    -0.020702261 % ln l -- 6
    1.024943397 % ln w -- 7
    -4.461247468 % ln r_k -- 8
    1.220529165 % ln y --9
    0           % nx (can be negative) -- 10
    0           % ca (can be negative) -- 11
    0.003924883 % r_star -- 12
    0           % b_F (can be negative) -- 13
    0.014022286 % ln r -- 14
    0.0         % b  (can be negative) -- 15
    -0.2         % ln g -- 16
    0.18      % tr (can be negative) -- 17
    3.1845];    % p_b -- 18


% Structure 1 (z_ss low, gy_ss)

S_.options.gy = 1 ;  % sets gy_ss

st_1 = structural_matrices(S_);
%
[   S0_1, ...
    S1_1, ...
    S2_1, ...
    ~, ...
    ~, ...
    ~, ...
    ~, ...
    unique_1 ] = Smats(st_1.C_, st_1.GAM0, st_1.GAM1, st_1.PSI, st_1.PPI);

if unique_1 ~= 1
    out = penalty;
    return
end

if isempty(unique_1)
    out = penalty;
    return
end

y_bar_sims_1 = (eye(n)-S1_1) \ S0_1;

% BP Solution matrices
C_1 = S0_1(1:n_,:);
Q_1 = S1_1(1:n_,1:n_);
G_1 = S2_1(1:n_,:);

y_bar_1 = (eye(n_)-Q_1) \ C_1;

%% Structure 2 (z_ss high, gy_ss low , g_ss at structure 1)

S_.param.z_ss = S_.param.z_ss + S_.param.Delta_z ; 
S_.param.tau_w_ss = S_.param.tau_w_ss + S_.param.Delta_tauw;
S_.param.tau_K_ss = S_.param.tau_K_ss + S_.param.Delta_tauK;

S_.param.g_ss = st_1.y_steady(16);

S_.initvals = [0.834755321 % ln c -- 1
              -0.095030773 % ln invest -- 2
              -0.973976034 % ln lambda -- 3
               0.005224883 % ln r_F -- 4
               4.631985754 % ln k -- 5
              -0.020702261 % ln l -- 6
               1.024943397 % ln w -- 7
              -4.461247468 % ln r_k -- 8
               1.220529165 % ln y --9
               0           % nx (can be negative) -- 10
               0           % ca (can be negative) -- 11
               0.003924883 % r_star -- 12
               0           % b_F (can be negative) -- 13
               0.014022286 % ln r -- 14
               0.0         % b  (can be negative) -- 15
              -1.6         % ln gy -- 16
               0.18      % tr (can be negative) -- 17
               3.1845];    % p_b -- 18

S_.options.gy = 0 ;  % sets g_ss

st_2 = structural_matrices(S_);

[   S0_2, ...
    S1_2, ...
    S2_2, ...
    ~, ...
    ~, ...
    ~, ...
    ~, ...
    unique_2 ] = Smats(st_2.C_, st_2.GAM0, st_2.GAM1, st_2.PSI, st_2.PPI);

if unique_2 ~= 1
   out = penalty;
   return
end

if isempty(unique_2)
   out = penalty;
   return
end

% BP Solution matrices
C_2 = S0_2(1:n_,:);
Q_2 = S1_2(1:n_,1:n_);
G_2 = S2_2(1:n_,:);

ytT_nobreak = zeros(n_,T_sample);
ytT_nobreak(:,1) = ytT(:,1);

y_non = y_bar_1;
y_non_nobreak = y_bar_1;

for t = 2:T_sample
    ytT_nobreak(:,t) = C_1 + Q_1*ytT_nobreak(:,t-1) + G_1*n_T(:,t);
    y_non_nobreak(:,t) = C_1 + Q_1*y_non_nobreak(:,t-1);
end

for t = 2:Tz-1
    y_non(:,t) = C_1 + Q_1*y_non(:,t-1);
end

for t = Tz:T_sample
    y_non(:,t) = C_2 + Q_2*y_non(:,t-1);
end

%% Construct Z(t)

Z = zeros(1, T_sample) ;
Z_nobreak = zeros(1,T_sample);

Z(1,1) = exp(log(GDP_pc(1)) - ytT(v_.y,1)) ;
Z_nobreak(1,1) = Z(1,1);

for t_ = 2:T_sample
    Z(1,t_) = Z(1,t_-1) * exp(ytT(v_.z,t_));
    Z_nobreak(1,t_) = Z_nobreak(1,t_-1) * exp(ytT_nobreak(v_.z,t_));
end


Z_non = zeros(1, T_sample);
Z_non_nobreak = zeros(1,T_sample);

Z_non(1,1) = exp(log(GDP_pc(1)) - y_non(v_.y,1)) ;
Z_non_nobreak(1,1) = Z_non(1,1);

for t_ = 2:T_sample
    Z_non(1,t_) = Z_non(1,t_-1) * exp(y_non(v_.z,t_));
    Z_non_nobreak(1,t_) = Z_non_nobreak(1,t_-1) * exp(y_non_nobreak(v_.z,t_));
end

%% Labour income tax revenue

labIncome_tax = ytT(v_.tau_w,:) .* Z.* exp(ytT(v_.w,:)) .* exp(ytT(v_.l,:)) ;

labIncome_tax_nobreak = ytT_nobreak(v_.tau_w,:) .* Z_nobreak.* exp(ytT_nobreak(v_.w,:)) .* exp(ytT_nobreak(v_.l,:)) ;

labIncome_tax_tau = ytT_nobreak(v_.tau_w,:) .* Z.* exp(ytT(v_.w,:)) .* exp(ytT(v_.l,:)) ;
labIncome_tax_w = ytT(v_.tau_w,:) .* Z.* exp(ytT_nobreak(v_.w,:)) .* exp(ytT(v_.l,:)) ;
labIncome_tax_l = ytT(v_.tau_w,:) .* Z.* exp(ytT(v_.w,:)) .* exp(ytT_nobreak(v_.l,:)) ;
labIncome_tax_Z = ytT(v_.tau_w,:) .* Z_nobreak.* exp(ytT(v_.w,:)) .* exp(ytT(v_.l,:)) ;

labIncome_tax_non = y_non(v_.tau_w,:) .* Z_non.* exp(y_non(v_.w,:)) .* exp(y_non(v_.l,:)) ;

labIncome_tax_non_nobreak = y_non_nobreak(v_.tau_w,:) .* Z_non_nobreak.* exp(y_non_nobreak(v_.w,:)) .* exp(y_non_nobreak(v_.l,:)) ;

labIncome_tax_non_tau = y_non_nobreak(v_.tau_w,:) .* Z_non.* exp(y_non(v_.w,:)) .* exp(y_non(v_.l,:)) ;
labIncome_tax_non_w = y_non(v_.tau_w,:) .* Z_non.* exp(y_non_nobreak(v_.w,:)) .* exp(y_non(v_.l,:)) ;
labIncome_tax_non_l = y_non(v_.tau_w,:) .* Z_non.* exp(y_non(v_.w,:)) .* exp(y_non_nobreak(v_.l,:)) ;
labIncome_tax_non_Z = y_non(v_.tau_w,:) .* Z_non_nobreak.* exp(y_non(v_.w,:)) .* exp(y_non(v_.l,:)) ;

%% Capital income tax revenue

capIncome_tax = ytT(v_.tau_K,:) .* Z.* exp(ytT(v_.r_K,:)) .* exp(ytT(v_.kk,:)) ;

capIncome_tax_nobreak = ytT_nobreak(v_.tau_K,:) .* Z_nobreak.* exp(ytT_nobreak(v_.r_K,:)) .* exp(ytT_nobreak(v_.kk,:)) ;

capIncome_tax_tau = ytT_nobreak(v_.tau_K,:) .* Z.* exp(ytT(v_.r_K,:)) .* exp(ytT(v_.kk,:)) ;
capIncome_tax_r = ytT(v_.tau_K,:) .* Z.* exp(ytT_nobreak(v_.r_K,:)) .* exp(ytT(v_.kk,:)) ;
capIncome_tax_k = ytT(v_.tau_K,:) .* Z.* exp(ytT(v_.r_K,:)) .* exp(ytT_nobreak(v_.kk,:)) ;
capIncome_tax_Z = ytT(v_.tau_K,:) .* Z_nobreak.* exp(ytT(v_.r_K,:)) .* exp(ytT(v_.kk,:)) ;


capIncome_tax_non = y_non(v_.tau_K,:) .* Z_non.* exp(y_non(v_.r_K,:)) .* exp(y_non(v_.kk,:)) ;

capIncome_tax_non_nobreak = y_non_nobreak(v_.tau_K,:) .* Z_non_nobreak.* exp(y_non_nobreak(v_.r_K,:)) .* exp(y_non_nobreak(v_.kk,:)) ;

capIncome_tax_non_tau = y_non_nobreak(v_.tau_K,:) .* Z_non.* exp(y_non(v_.r_K,:)) .* exp(y_non(v_.kk,:)) ;
capIncome_tax_non_r = y_non(v_.tau_K,:) .* Z_non.* exp(y_non_nobreak(v_.r_K,:)) .* exp(y_non(v_.kk,:)) ;
capIncome_tax_non_k = y_non(v_.tau_K,:) .* Z_non.* exp(y_non(v_.r_K,:)) .* exp(y_non_nobreak(v_.kk,:)) ;
capIncome_tax_non_Z = y_non(v_.tau_K,:) .* Z_non_nobreak.* exp(y_non(v_.r_K,:)) .* exp(y_non(v_.kk,:)) ;

%% Consumption tax revenue

consumption_tax = ytT(v_.tau_c,:) .* Z.* exp(ytT(v_.c,:)) ;

consumption_tax_nobreak = ytT_nobreak(v_.tau_c,:) .* Z_nobreak.* exp(ytT_nobreak(v_.c,:)) ;

consumption_tax_tau = ytT_nobreak(v_.tau_c,:) .* Z.* exp(ytT(v_.c,:)) ;
consumption_tax_c = ytT(v_.tau_c,:) .* Z.* exp(ytT_nobreak(v_.c,:))  ;
consumption_tax_Z = ytT(v_.tau_c,:) .* Z_nobreak.* exp(ytT(v_.c,:))  ;


consumption_tax_non = y_non(v_.tau_c,:) .* Z_non.* exp(y_non(v_.c,:)) ;

consumption_tax_non_nobreak = y_non_nobreak(v_.tau_c,:) .* Z_non_nobreak.* exp(y_non_nobreak(v_.c,:)) ;

consumption_tax_non_tau = y_non_nobreak(v_.tau_c,:) .* Z_non.* exp(y_non(v_.c,:)) ;
consumption_tax_non_c = y_non(v_.tau_c,:) .* Z_non.* exp(y_non_nobreak(v_.c,:))  ;
consumption_tax_non_Z = y_non(v_.tau_c,:) .* Z_non_nobreak.* exp(y_non(v_.c,:))  ;

%% Government spending

spending =  Z.* exp(ytT(v_.g,:)) ;

spending_nobreak = exp(ytT_nobreak(v_.g,:)) .* Z_nobreak ;

spending_g = exp(ytT_nobreak(v_.g,:)) .* Z ;
spending_Z = exp(ytT(v_.g,:)) .* Z_nobreak  ;

spending_non =  Z_non.* exp(y_non(v_.g,:)) ;

spending_non_nobreak = exp(y_non_nobreak(v_.g,:)) .* Z_non_nobreak ;

spending_non_g = exp(y_non_nobreak(v_.g,:)) .* Z_non ;
spending_non_Z = exp(y_non(v_.g,:)) .* Z_non_nobreak  ;


%% Lump sum Taxes

lumpsum_tax =  Z.* ytT(v_.tr,:) ;

lumpsum_tax_nobreak = ytT_nobreak(v_.tr,:) .* Z_nobreak ;

lumpsum_tax_tr = ytT_nobreak(v_.tr,:) .* Z ;
lumpsum_tax_Z = ytT(v_.tr,:) .* Z_nobreak  ;


lumpsum_tax_non =  Z_non.* y_non(v_.tr,:) ;

lumpsum_tax_non_nobreak = y_non_nobreak(v_.tr,:) .* Z_non_nobreak ;

lumpsum_tax_non_tr = y_non_nobreak(v_.tr,:) .* Z_non ;
lumpsum_tax_non_Z = y_non(v_.tr,:) .* Z_non_nobreak  ;

%% Surplus

surplus =  Z.* ytT(v_.s,:) ;

surplus_nobreak = ytT_nobreak(v_.s,:) .* Z_nobreak ;

surplus_s = ytT_nobreak(v_.s,:) .* Z ;
surplus_Z = ytT(v_.s,:) .* Z_nobreak  ;


surplus_non =  Z_non.* y_non(v_.s,:) ;

surplus_non_nobreak = y_non_nobreak(v_.s,:) .* Z_non_nobreak ;

surplus_non_s = y_non_nobreak(v_.s,:) .* Z_non ;
surplus_non_Z = y_non(v_.s,:) .* Z_non_nobreak  ;


%% Plot nonstochastic responses

t_ = 1983:.25:2018;
t_start = find(t_==2000);

figure('Name', "Fiscal Decomposition", 'Position', [1 1 800 500])
set(0,'defaulttextinterpreter','latex')

subplot(2,2,1), plot(t_, (labIncome_tax_non./labIncome_tax_non_nobreak -1)*100,'Color', [0.6 0.6 0.6], 'linewidth',2,'Marker','o','MarkerSize',2), title('Labour Income Tax Revenue') 
    hold on
    plot(t_,(labIncome_tax_non./labIncome_tax_non_tau -1)*100, 'Color',[0.44, 0.16, 0.39], 'LineWidth', 1.2)
    plot(t_,(labIncome_tax_non./labIncome_tax_non_w-1)*100, 'Color',[0.65, 0.16, 0.16], 'LineWidth', 1.2)
    plot(t_,(labIncome_tax_non./labIncome_tax_non_l-1)*100, 'Color',[1.0, 0.75, 0.0], 'LineWidth', 1.2)
    plot(t_,(labIncome_tax_non./labIncome_tax_non_Z-1)*100, 'Color',[0.0, 0.4, 0.65], 'LineWidth', 1.2)
    xlim([2000 2018]);
     xticks([2000, 2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    xticklabels([2000,2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    ylim([-25,10]);
    ylabel('\%' )
    set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')
    lgnd = legend( '$\tau_w w l Z$' ,'$\tau_w$', '$w$', '$l$', '$Z$');
    set(lgnd,'Box','off', 'Interpreter','latex', 'NumColumns',5, 'location','southoutside')

    
subplot(2,2,2), plot(t_, (capIncome_tax_non./capIncome_tax_non_nobreak-1)*100, 'Color',[0.6 0.6 0.6], 'linewidth',2,'Marker','o','MarkerSize',2), title('Capital Income Tax Revenue') 
    hold on
    plot(t_,(capIncome_tax_non./capIncome_tax_non_tau-1)*100, 'Color',[0.44, 0.16, 0.39], 'LineWidth', 1.2)
    plot(t_,(capIncome_tax_non./capIncome_tax_non_r-1)*100, 'Color',[0.65, 0.16, 0.16], 'LineWidth', 1.2)
    plot(t_,(capIncome_tax_non./capIncome_tax_non_k-1)*100, 'Color',[1.0, 0.75, 0.0], 'LineWidth', 1.2)
    plot(t_,(capIncome_tax_non./capIncome_tax_non_Z-1)*100, 'Color',[0.0, 0.4, 0.65], 'LineWidth', 1.2)
    xlim([2000 2018]);
     xticks([2000, 2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    xticklabels([2000,2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    ylabel('\%' )
    set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')
    lgnd = legend( '$\tau_K r^K k Z$' ,'$\tau_K$', '$r^K$', '$k$', '$Z$');
    set(lgnd,'Box','off', 'Interpreter','latex', 'NumColumns',5, 'location','southoutside')

    
subplot(2,2,3), plot(t_, (consumption_tax_non./consumption_tax_non_nobreak-1)*100, 'Color',[0.6 0.6 0.6], 'linewidth',2,'Marker','o','MarkerSize',2), title('Consumption Tax Revenue') 
    hold on
    plot(t_,(consumption_tax_non./consumption_tax_non_tau-1)*100, 'Color',[0.44, 0.16, 0.39], 'LineWidth', 1.2)
    plot(t_,(consumption_tax_non./consumption_tax_non_c-1)*100, 'Color',[1.0, 0.75, 0.0], 'LineWidth', 1.2)
    plot(t_,(consumption_tax_non./consumption_tax_non_Z-1)*100, 'Color',[0.0, 0.4, 0.65], 'LineWidth', 1.2)
    xlim([2000 2018]);
     xticks([2000, 2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    xticklabels([2000,2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])    
    ylim([-25,10]);
    ylabel('\%' )
    set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')
    lgnd = legend( '$\tau_c c Z$' ,'$\tau_c$', '$c$', '$Z$');
    set(lgnd,'Box','off', 'Interpreter','latex', 'NumColumns',5, 'location','southoutside')
    
subplot(2,2,4), plot(t_, (spending_non./spending_non_nobreak-1)*100, 'Color',[0.6 0.6 0.6], 'linewidth',2,'Marker','o','MarkerSize',2), title('Government Spending') 
    hold on
    plot(t_,(spending_non./spending_non_g-1)*100, 'Color',[1.0, 0.75, 0.0], 'LineWidth', 1.2)
    plot(t_,(spending_non./spending_non_Z-1)*100, 'Color',[0.0, 0.4, 0.65], 'LineWidth', 1.2)
    xlim([2000 2018]);
    xticks([2000, 2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    xticklabels([2000,2002, 2004,2006, 2008, 2010,2012,2014,2016,2018])
    ylabel('\%' )
    set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')
    lgnd = legend( '$g Z$', '$g$', '$Z$');
    set(lgnd,'Box','off', 'Interpreter','latex', 'NumColumns',5, 'location','southoutside')
    
set(gcf,'Color','w')
set(findall(gcf,'-property','FontSize'),'FontSize',12);
print -depsc '-r600' Figure_11.eps
