% 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 estimation results
load mhall_RBC.mat

%% Data
aus_data;
T = size(Zobs,2);

%% Pick draws

rng('default')
rng(20);
RandStream.getGlobalStream;

% Cut out the burnin draws
xx_J = params_xx_out(:,:);
N_sort = length(xx_J);

% Select the draws to sample over
JJ = 1000;        % Number of draws
J = round(N_sort*rand(JJ,1));
J = sort(J);
Reps = JJ;

% Posterior distribution of the parameters
xx_J = params_xx_out(:,:);
xx_J = xx_J(J,:);

% Posterior distribution of the date breaks
T_J =  params_T_out(:,:);
T_J = T_J(J,:);


%% Variable and dimensions

v_ = var_dim ;

n     = v_.dim.n    ;
n_    = v_.dim.n_    ;
ns    = v_.dim.ns    ;


T_long = T + 180; % Longer-run projections

Y_non_dy = zeros(Reps,T);
Y_non_gy = zeros(Reps,T);
Y_non_dinvest = zeros(Reps,T);
Y_non_dc = zeros(Reps,T);
Y_non_nx = zeros(Reps,T);
Y_non_dw = zeros(Reps,T);
Y_non_r = zeros(Reps,T);
Y_non_r_star = zeros(Reps,T);
Y_non_by = zeros(Reps,T);
Y_non_rev_cy = zeros(Reps,T);
Y_non_rev_wy = zeros(Reps,T);
Y_non_rev_Ky = zeros(Reps,T);
Y_non_tr = zeros(Reps,T);
Y_non_s = zeros(Reps,T);
Y_non_y = zeros(Reps,T);
Y_non_w = zeros(Reps,T);
Y_non_l = zeros(Reps,T);
Y_non_rK = zeros(Reps,T);
Y_non_k = zeros(Reps,T);

Y_non_dy_long = zeros(Reps,T_long);
Y_non_gy_long = zeros(Reps,T_long);
Y_non_dc_long = zeros(Reps,T_long);
Y_non_nx_long = zeros(Reps,T_long);
Y_non_dw_long = zeros(Reps,T_long);
Y_non_r_long = zeros(Reps,T_long);
Y_non_rstar_long = zeros(Reps,T_long);
Y_non_by_long = zeros(Reps,T_long);
Y_non_rev_cy_long = zeros(Reps,T_long);
Y_non_rev_wy_long = zeros(Reps,T_long);
Y_non_rev_Ky_long = zeros(Reps,T_long);
Y_non_k_long =  zeros(Reps,T_long);


y_non_all = zeros(n,T_long,Reps);


for gg = 1:Reps

    x_ = xx_J(gg,:);

    %% Parameters

    S_.param = param_func(x_);

    %% 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_,:);


    %% Plot simulated data against data

    y_non = zeros(n,T);
    y_non_long = zeros(n,T_long);

    Tz   = T_J(gg,1);
    Tsig = T_J(gg,2);

    y_non(:,1) = y_bar_sims_1;
    y_non_long(:,1) = y_bar_sims_1;

    % with the break in z_ss
    for t = 2:Tz-1
        y_non(:,t) = S0_1 + S1_1*y_non(:,t-1);
        y_non_long(:,t) = S0_1 + S1_1*y_non_long(:,t-1);
    end
    for t = Tz:T
        y_non(:,t) = S0_2 + S1_2*y_non(:,t-1);
        y_non_long(:,t) = S0_2 + S1_2*y_non_long(:,t-1);
    end
    for t = T+1:T_long
        y_non_long(:,t) = S0_2 + S1_2*y_non_long(:,t-1);
    end


    % Store all variables
    y_non_all(:,:,gg) = y_non_long;


    Y_non_dy(gg,:)      = y_non(v_.dy_obs,:);
    Y_non_dinvest(gg,:) = y_non(v_.dinvest_obs,:);
    Y_non_dc(gg,:)      = y_non(v_.dc_obs,:);
    Y_non_nx(gg,:)      = y_non(v_.nx_obs,:);
    Y_non_r(gg,:)       = y_non(v_.r_obs,:);
    Y_non_r_star(gg,:)  = y_non(v_.r_star_obs,:);
    Y_non_gy(gg,:)      = y_non(v_.gy_obs,:);
    Y_non_by(gg,:)      = y_non(v_.by_obs,:);
    Y_non_dw(gg,:)      = y_non(v_.dw_obs,:);
    Y_non_rev_cy(gg,:)  = y_non(v_.rev_c_y_obs,:);
    Y_non_rev_wy(gg,:)  = y_non(v_.rev_w_y_obs,:);
    Y_non_rev_Ky(gg,:)  = y_non(v_.rev_K_y_obs,:);
    Y_non_tr(gg,:)  = y_non(v_.tr,:);
    Y_non_s(gg,:)  = y_non(v_.s,:);
    Y_non_y(gg,:)  = y_non(v_.y,:);
    Y_non_w(gg,:)  = y_non(v_.w,:);
    Y_non_l(gg,:)  = y_non(v_.l,:);
    Y_non_rK(gg,:)  = y_non(v_.r_K,:);
    Y_non_k(gg,:)  = y_non(v_.kk,:);

    Y_non_dy_long(gg,:)      = y_non_long(v_.dy_obs,:);
    Y_non_dc_long(gg,:) = y_non_long(v_.dc_obs,:);
    Y_non_nx_long(gg,:)      = y_non_long(v_.nx_obs,:);
    Y_non_r_long(gg,:)       = y_non_long(v_.r_obs,:);
    Y_non_rstar_long(gg,:)       = y_non_long(v_.r_star_obs,:);
    Y_non_gy_long(gg,:)      = y_non_long(v_.gy_obs,:);
    Y_non_by_long(gg,:)      = y_non_long(v_.by_obs,:);
    Y_non_dw_long(gg,:)      = y_non_long(v_.dw_obs,:);
    Y_non_rev_cy_long(gg,:)  = y_non_long(v_.rev_c_y_obs,:);
    Y_non_rev_wy_long(gg,:)  = y_non_long(v_.rev_w_y_obs,:);
    Y_non_rev_Ky_long(gg,:)  = y_non_long(v_.rev_K_y_obs,:);
    Y_non_k_long(gg,:)  = y_non_long(v_.kk,:);


end



%% Figure

t_long = zeros(T_long,1);
t_long(1:T) = t_;
for i_ = (size(t_)+1):T_long, t_long(i_) = t_long(i_-1) + 0.25; end

[mat_quant_dy, matm_dy, num_quant_dy, ylimit_dy] = fanChart(100*Y_non_dy_long',zeros(1,1)', [-2 2]);
[mat_quant_dc, matm_dc, num_quant_dc, ylimit_dc] = fanChart(100*Y_non_dc_long',zeros(1,1)', [-2 2]);
[mat_quant_nx, matm_nx, num_quant_nx, ylimit_nx] = fanChart(100*Y_non_nx_long',zeros(1,1)', [-2 2]);
[mat_quant_dw, matm_dw, num_quant_dw, ylimit_dw] = fanChart(100*Y_non_dw_long',zeros(1,1)', [-2 2]);
[mat_quant_r, matm_r, num_quant_r, ylimit_r] = fanChart(400*Y_non_r_long',zeros(1,1)', [-2 2]);
[mat_quant_rstar, matm_rstar, num_quant_rstar, ylimit_rstar] = fanChart(400*Y_non_rstar_long',zeros(1,1)', [-2 2]);
[mat_quant_gy, matm_gy, num_quant_gy, ylimit_gy] = fanChart(100*exp(Y_non_gy_long)',zeros(1,1)', [-2 2]);
[mat_quant_by, matm_by, num_quant_by, ylimit_by] = fanChart(100*(Y_non_by_long)'/4,zeros(1,1)', [-2 2]);
[mat_quant_revcy, matm_revcy, num_quant_revcy, ylimit_revcy] = fanChart(100*(Y_non_rev_cy_long)',zeros(1,1)', [-2 2]);
[mat_quant_revwy, matm_revwy, num_quant_revwy, ylimit_revwy] = fanChart(100*(Y_non_rev_wy_long)',zeros(1,1)', [-2 2]);
[mat_quant_revKy, matm_revKy, num_quant_revKy, ylimit_revKy] = fanChart(100*(Y_non_rev_Ky_long)',zeros(1,1)', [-2 2]);
[mat_quant_k, matm_k, num_quant_k, ylimit_k] = fanChart(exp(Y_non_k_long)',zeros(1,1)', [-2 2]);

median_forecast = 10;

left_color = [0.8 0.8 0.8];
right_color = [0.5 0.5 0.5];
cmap = interp1([0, 1], [left_color; right_color], linspace(0, 1, 10));
cmap = [cmap; flipud(cmap)];

figure('Name', "Non-Stochastic Transition Fan Chart", 'Position', [1 1 1050 1000])
%sgtitle('NonStochastic Transition Fan Chart')
set(0,'defaulttextinterpreter','latex')

subplot(4,3,1)
h=area(matm_dy(:,:),-2);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_dy(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(1,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Output Growth')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,2)
h=area(matm_dc(:,:),-2);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_dc(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(2,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Consumption Growth')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,3)
h=area(matm_nx(:,:),0);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_nx(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(3,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
title('Net Exports-to-Output Ratio')
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,4)
h=area(matm_dw(:,:),-4);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_dw(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(4,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Wage Growth')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,5)
h=area(matm_r(:,:),-0.8);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_r(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(400*(Zobs(5,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Domestic Real Interest Rate')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,6)
h=area(matm_rstar(:,:),-2.6);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_rstar(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(400*(Zobs(6,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Foreign Real Interest Rate')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,7)
h=area(matm_gy(:,:),20);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_gy(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*exp(Zobs(7,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Government Spending-to-Output Ratio')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,8)
h=area(matm_by(:,:),0);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_by(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(8,:)/4), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Public Debt-to-Output Ratio')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,9)
h=area(matm_revcy(:,:),3);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_revcy(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(9,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Consumption Tax Revenue-to-Output Ratio')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,10)
h=area(matm_revwy(:,:),9.4);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_revwy(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(10,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Labour Income Tax Revenue-to-Output Ratio')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,11)
h=area(matm_revKy(:,:),2);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_revKy(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
hold on
plot(100*(Zobs(11,:)), 'Color','k', 'LineWidth', 0.8)
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Capital Income Tax Revenue-to-Output Ratio')
ylabel('\%')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')

subplot(4,3,12)
h=area(matm_k(:,:),20);
set(h,'LineStyle','none')
set(h(1),'FaceColor',[1 1 1])
for i =2:19
    set(h(i),'FaceColor',cmap(i,:))
end
hold on
plot(mat_quant_k(1:end,median_forecast),'Color',[0.4 0.4 0.4],'Linewidth',1.2);
xticks([1 20 40 60 80 100 120 140 160 180])
xticklabels({1983,1988,1993,1998,2003,2008,2013,2018,2023,2028})
xlim([1 180])
title('Capital Accumulation')
ax = gca;
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
set(gca, 'Box', 'off', 'TickDir', 'in', 'TickLabelInterpreter','latex')


set(gcf,'Color','w')
set(findall(gcf,'-property','FontSize'),'FontSize',12);
print -depsc '-r600' Figure_5.eps